We have developed a combined chemical and hydrodynamical model of cometary comae. The model is described and used to follow the deuterium chemistry occurring in the coma after the sublimation of D-rich nuclear ices. It is shown that the molecular D/H ratios in parent species remain unchanged in the coma gas, confirming that observations of such molecules can be used to infer the D/H ratio in the nuclear ices. The D/H ratios in daughter molecules reflect that of their parents, and so can potentially be used to elucidate the formation mechanisms of these species. In particular, observations of DNC may be able to resolve the uncertainty surrounding the origin of HNC in cometary comae.
Comets are thought to be the most pristine objects in our Solar system. As such, their chemical composition reflects both the conditions in the outer solar nebula some 4.5×109 yr ago, and the nature of the natal interstellar cloud from which the Solar system was formed. Assuming that the chemistry in the natal cloud was the same as in dark clouds in the local interstellar medium, then observations of such clouds tell us the initial chemical composition of the protosolar nebula. Thus a detailed comparison of cometary and interstellar ices can be used to constrain models of Solar system formation, since any differences must be due to processes occurring as the nebula condensed [e.g. many models predict an accretion shock sufficiently strong to evaporate ice mantles from infalling grains (Chick & Cassen 1997; Lunine et al. 1991)].
Until recently, such comparisons were not possible, since very few data were available. However, the last few years have seen an explosion in our knowledge of both cometary and interstellar ices. The recent apparitions of two very bright comets – Hyakutake in 1996 and Hale–Bopp in 1997 – led to the detection of many molecules never before seen in cometary comae (e.g. Mumma et al. 1996; Irvine et al. 1996; Biver et al. 1997, 1999a). At the same time, the launch of the ISO satellite allowed the detailed investigation of many interstellar solid state absorption lines (e.g. Ehrenfreund 1999; Ehrenfreund & Schutte 2000). Together with theoretical and observational progress in the understanding of ‘hot cores’ (e.g. Charnley 1997; Millar, Macdonald & Gibb 1997b; Rodgers & Charnley 2001a), this has given us a good understanding of the nature of interstellar ices. These new measurements have permitted comparisons to be made between interstellar and cometary ices (see e.g. Altwegg et al. 1999; Ehrenfreund & Charnley 2000; Irvine et al. 2000; Bockelée-Morvan et al. 2000). The consensus appears to be that, whilst comets retain a strong signature of their ultimate interstellar origins, significant processing must have occurred in the protosolar nebula.
If this is correct, one of the most important effects may be on the deuterium content of the ices. In dark molecular clouds, deuterated molecules are observed with abundances relative to their H-bearing analogues of up to 10 per cent (e.g. Howe et al. 1994; Roueff et al. 2000), over 1000 times the cosmic D/H ratio of ∼10−5. This is because replacing an H with a D typically lowers the ground-state energy level of a molecule by several hundred kelvin. At low temperatures the ions and undergo D/H exchange reactions with HD, becoming highly fractionated (Smith, Adams & Alge 1982). The enhanced deuterium levels in these ions are then transferred to other molecules by gas-phase chemical reactions. Models of deuterium chemistry show that the predicted fractionation is extremely sensitive to the gas temperature (Millar, Bennett & Herbst 1989; Rodgers & Millar 1996). Thus, if cometary volatiles were indeed evaporated from infalling grains and processed before refreezing, it may be expected that the molecular D/H ratios would be reduced in the warm gas relative to the original values set in the cold (10–20 K) interstellar phase.Geiss & Reeves 1981). Viscous mixing could have potentially transported material from the inner region to the cometary formation zone, thereby diluting the cometary HDO fractionation (Drouart et al. 1999; Mousis et al. 2000). Yung et al. (1988) showed that in the surface layer of the disc, photodestruction of H2O and HDO followed by the conversion of OH to OD by the exchange reaction with atomic D, followed by the re-formation of water, could also alter the HDO/H2O ratio. Again, if turbulent mixing were efficient, this process would affect the entire disc. An alternative approach was adopted by Aikawa & Herbst (1999, 2001), who developed a chemical kinetic model of the deuterium chemistry occurring as matter moves through the disc toward the central protostar. They showed that the chemistry in the disc is very similar to that in dark clouds, and that ion–molecule reactions in the outer disc may have modified the initial D/H ratios.
It is also possible that the observed cometary D/H ratios may be modified by gas-phase chemical reactions taking place in the coma. Early models of coma chemistry showed that reactions can occur rapidly in the inner coma, where the number densities can be as high as 1012 cm−3 (Giguere & Huebner 1978; Huebner & Giguere 1980; Biermann et al. 1981). The most important reactions were found to be proton transfer reactions. These can potentially cycle deuterium between the different coma molecules, altering the initial D/H ratios released from the nuclear ice. Thus it is necessary to construct accurate models of cometary deuterium chemistry, in order that gas-phase coma observations can be safely extrapolated to give the nuclear D/H ratios.
To date, only two deuterated molecules have been detected in cometary comae: HDO and DCN. The former was observed in comets Halley, Hyakutake and Hale–Bopp (Eberhardt et al. 1995; Bockelée-Morvan et al. 1998; Meier et al. 1998a), and in each comet the HDO/H2O ratio was ≈6×10−4 (this is sometimes expressed as a D/H ratio of 3×10−4, since the fact that water contains two hydrogens increases the HDO/H2O over the D/H ratio by a stoichiometric factor of 2). DCN was seen for the first time in Hale–Bopp, with a DCN/HCN ratio of 0.002 (Meier et al. 1998b). The fact that the DCN/HCN ratio in interstellar hot cores is also ∼0.002 (Hatchell, Millar & Rodgers 1998) has led some authors to conclude that cometary ices must consist of relatively unprocessed interstellar material. However, it is not clear that the DCN and HCN in hot cores are in fact injected directly from ices. Rather, HCN may be formed in the hot gas from reactions initiated by evaporated ammonia or N2 (Rodgers & Charnley 2001a), and the DCN/HCN ratio will depend on the relative amounts of N2, NH3 and NH2D in the gas. Kessler, Qi & Blake (2000) have detected DCN in protoplanetary discs, and derive DCN/HCN ratios of ∼0.01, much larger than the values seen in hot cores and comets, but consistent with the values seen in cold interstellar gas (Wootten 1987; Hatchell et al. 1998). Thus, this is further evidence for processing of D/H ratios in the protosolar nebula.
Although only two deuterated cometary molecules have so far been observed, the large number detected in the interstellar medium (over 15 to date) suggests it is likely that many more will be detected in the future. Also, the Stardust and Rosetta comet rendezvous space missions will allow cometary D/H ratios to be determined via in situ measurements (Wright & Pillinger 1998), and the FIRST/Herschel mission will allow the HDO/H2O ratio to be determined in a large sample of comets (Bockelée-Morvan & Crovisier 2001). We have therefore developed a model of the deuterium chemistry occurring in the comae of comets. Comparison of our results with future observations will help to elucidate many key questions regarding where, when and how comets formed.
A full description of our model appears in Appendix A; here we give a brief overview. In addition to calculating the chemistry, our model also follows the dynamical evolution of the gas after it sublimates from the nuclear surface. This is done because previous hydrodynamical models of the coma have shown that the gas temperatures show strong variations with cometocentric distance (e.g. Marconi & Mendis 1983, 1986; Crovisier 1984; Huebner 1985; Körösmezey et al. 1987; Schmidt et al. 1988; Crifo 1991). The initial expansion is adiabatic so the gas cools, but further out photodissociation and ionization reactions heat the gas. Eventually, the temperatures of the different components (ions, electrons and neutrals) can decouple. The electron temperature, in particular, can reach extremely high values, <104 K (Eberhardt & Krankowsky 1995).
Many previous chemical models have ignored these effects, and have simply assumed constant outflow velocities and temperatures throughout the coma (e.g. Giguere & Huebner 1978; Irvine et al. 1998). Since the gas temperature can have important consequences for the chemistry, especially with regard to deuterium fractionation reactions, we have calculated the gas velocity and temperatures in the coma. Also, it turns out that chemical reactions are an important heat source for the coma gas. As was mentioned in the Introduction, proton transfer reactions – which are typically exothermic by ∼1 eV – occur rapidly in the inner coma, and can heat the ion fluid well above what would otherwise be predicted in their absence (see Section 3.1 and Körösmezey et al. 1987). Therefore, for any model of the coma to be fully self-consistent, it must simultaneously solve both the chemical and hydrodynamical continuity equations.
Having said this, although we include all the major micro-physical processes occurring in the coma, and our results are in good agreement with other hydrodynamical models of comets (see Section 3.1), it must be stressed that the primary aim of our model is to follow the coma chemistry. A full multifluid hydrodynamic treatment requires an iterative calculation (e.g. Körösmezey et al. 1987), and the inclusion of magnetic fields, which necessitates abandoning the assumption of spherical symmetry (e.g. Wegmann et al. 1987). Both these requirements are not feasible for a chemical reaction scheme as large as the one in our model. So, although we have tried to ensure that our dynamical calculations are as accurate as possible, our model is not intended to be a state of the art physical model, but rather the hydrodynamical effects are included so as to derive realistic number density and temperature profiles in the coma.
The chemical reaction scheme contains 181 species and over 3500 reactions. It is based on an original scheme developed for modelling interstellar hot cores (Charnley, Tielens & Millar 1992), with additional reactions taken from the UMIST data base (Millar, Farquhar & Willacy 1997a). Photodissociation rates appropriate for the solar UV field were taken from Huebner, Keady & Lyon (1992). Some additional reactions (principally electron impact reactions) were taken from Schmidt et al. (1988). Also, some proton transfer reactions were added that are exothermic but are not included in the UMIST rate file; such reactions have a crucial role in determining the chemical production of molecules such as HNC (e.g. Rodgers & Charnley 1998, 2001c) and large O-bearing organics (Rodgers & Charnley 2001a,b; Charnley et al. 2001).
The initial chemical composition of the nuclear ice was based where possible on observations of Hyakutake and Hale–Bopp, and is shown in Table 1. Various values for the initial molecular D/H ratios in the ice were assumed, and the dependence of the results on these values is discussed in Section 3.3. The initial temperature of the gas was given by (Biver et al. 1997)
The original chemical reaction scheme (which comprised some 1200 reactions) was extended to include all analogous singly deuterated reactions using the techniques described by Rodgers & Millar (1996). Essentially, the rules used to generate the deuterated reaction set were as follows.
The total rate coefficients for deuterated reactions were assumed to be the same as for their H-bearing analogues.
Where more than one hydrogen is involved in a reaction, and their positions on the product species are not symmetrical, the branching ratios for the D finishing up in each position are proportional to the statistical probabilities (except for some photodissociation and recombination reactions – see below).
Note that the use of the word ‘involved’ in rule 2 implies that generating a deuterated rate file is necessarily a subjective process. Only if the reaction mechanism is known in detail is it possible to say which hydrogens take part in the reaction and which are passive spectators, otherwise one must make informed guesses.
As in Rodgers & Millar (1996), functional groups present on both products and reactants were assumed not to be involved in the reaction. For example, photodissociation of CH3OD is assumed to produce CH3 and OD, whereas CH2DOH leads to CH2D and OH, in agreement with experimental results (Harich et al. 1999). Similarly, dissociative recombination reactions are assumed to preserve the structure of the molecule. Thus HCND+ recombines to give HCN or DNC but not DCN or HNC (Millar 1992). Nevertheless, this leaves many reactions for which the branching ratios will depend on the preferences of the person or the software generating the deuterated rate file. For example, Markwick, Charnley & Millar (2001) assume that the reaction of CD+ with CH3OH can form either or CH3OHD+, whereas in our model only the second route is allowed. Although this is not a particularly important reaction, because there are so many reactions in a typical rate file where the final position of the D is ambiguous, there are many possible deuterated versions of a specific rate file. Hence the model results may well depend crucially upon the assumptions that were made when deuterating the rate file. In the long term, these uncertainties can only be resolved through laboratory measurements of deuterated reactions. However, by identifying the key reactions involved in a specific process, and running the model for different values of the branching ratios in these reactions, one can gain an understanding of the relative importance of the original assumptions.
Photodissociation cross-sections and branching ratios for a few deuterated molecules have been studied both theoretically and in the laboratory. Zhang & Imre (1988) showed that the cross-section for photodissociation of HDO is similar to that of H2O, but that OD+H production is favoured over OH+D. The exact value of the OD/OH branching ratio will depend on the spectrum of the incident radiation; where both cross-sections peak (λ≈159–172 nm) the value is about 2.5 (cf. Bar et al. 1991). At longer wavelengths, or if the HDO is initially in a vibrationally excited state, the OD/OH ratio can be even higher (Zhang, Imre & Frederick 1989; Vander Wal, Scott & Crim 1990; Vander Wal et al. 1991; Plusquellic, Votava & Nesbitt 1998). Experimental studies of the photodissociation of deuterated acetylene, methane and ammonia also show a predominance of H as opposed to D ejection (Wang, Hsu & Liu 1997; Wang & Liu 1998; Reid, Loomis & Leone 2000a,b), in this case by factors of ∼2–4. However, in most cases, the total photodissociation cross-sections for the deuterated molecules are similar to their non-deuterated counterparts (cf. Brownsword et al. 1997). Finally, note that Shin et al. (1996) studied the photodissociation of deuterated acetaldehyde ions and showed that some D/H ‘scrambling’ can occur. Our assumption that deuterated reactions always preserve the integrity of functional groups is therefore likely to be somewhat naive.
Given the scant data that exist on the photodissociation rates of deuterated molecules, we have assumed that these reactions occur at the same rates as their H-bearing analogues. Two different assumptions regarding the branching ratios were made. In the first model, these ratios were assumed to be statistical, as for all other reactions involving deuterated species. In a second model, we assumed that cleavage of the X–H bond is three times more likely than for theX–D bond, where X represents O, C or N. There is also evidence for preferential ejection of H atoms in electron dissociative recombination reactions (Gellene & Porter 1984; Larsson et al. 1996; Jensen et al. 1999, 2000). For those ions where the branching ratios have been measured we use the experimental values, but we have continued to assume statistical ejection probabilities for all other dissociative recombination reactions.
In contrast to the overall reaction cross-sections, the energetics of several deuterated photodissociation reactions have been well studied in the laboratory. In particular, the H2O/D2O system has been investigated extensively (e.g. Zanganeh et al. 2000; Harich et al. 2001). These studies demonstrated that the kinetic and internal energies of the emitted D and OD are almost identical to those of the H and OH produced in the non-deuterated reaction. It is important to know the energy distributions of the ejected atoms, since they may be important in driving a supra-thermal chemistry in the coma (Rodgers & Charnley 1998; see Section A6.2). In addition, these studies have also shown that the branching ratios for the three-body channels are similar for D2O and H2O dissociation (Harich et al. 2001).
Finally, several D/H exchange reactions were included; these are listed in table 2 of Rodgers & Millar (1996). One alteration that we have made is to remove the reaction HCN+D↔DCN+H, since recent calculations by Herbst & Talbi (1998) suggest that it has an activation energy of ∼5000 K. We have also added the following reaction:
Fig. 1 shows the temperature and velocity profiles produced by our model for a Hyakutake-like comet at 1 au (Q0=1.7×1029 molecule s−1,Mumma et al. 1996; R0=2.2 km, Sarmecanic et al. 1997; Lisse et al. 1999). The results are in excellent agreement with previous models of the coma (Crifo 1991). We reproduce the initial acceleration of the gas, and the temperature drop in the inner coma. Beyond ∼100 km solar photons deposit large amounts of energy into the gas and the temperatures rise.
Fig. 2 shows the heating and cooling mechanisms in the coma, also for a Hyakutake-like comet at 1 au. Except for the very inner coma where the optical depth is greater than 1, the neutrals are heated mainly by photodissociation reactions. Within ∼3000 km, all the energy from such reactions is deposited locally in the coma, so the heating rate varies as R−2. Beyond this distance, fast H atoms and molecules carry off increasing amounts of energy without being thermalized, leading to a sharper fall-off in the heating rate.
The ions are heated mainly by exothermic proton transfer reactions. In the inner coma, each photoionization of H2O to H2O+ ultimately results in (either by direct reaction with NH3 or via a succession of steps), which releases 2.7 eV of energy into the gas. Around half this energy goes to the neutral fluid, but the low fractional ionization means that the mean energy per particle is much larger for the ions. This explains why the ion temperature never drops as low as the neutrals and electrons in the inner coma, and also shows the importance of including chemical reactions in hydrodynamic models (cf. Körösmezey et al. 1987).
The electron temperature is very strongly coupled to the neutrals in the inner coma, owing to inelastic scattering with H2O molecules. Almost all the energy imparted to the electrons when they are created in photoionization reactions is immediately transferred to the neutrals. Further out, the electrons can no longer lose their energy as the number densities decrease and inelastic scattering becomes inefficient, so the electron temperature becomes extremely large.
Fig. 3 shows the abundances of a number of important species. Rather than plotting number densities, the graphs show the value 4πnvR2, which corresponds to the total number of molecules per second passing through a spherical shell at radius R. Such plots make it easier to analyse the chemistry, since this removes the effects of the R−2 geometric dilution, and any stretching or compression arising from acceleration or deceleration of the gas. Thus, if a curve in Fig. 3 is rising (falling), then that species must be being created (destroyed) by chemical reactions.
The top panel shows the abundances of ‘parent’ species, present initially in the nuclear ice mixture. Most such species survive for more than 105 km, consistent with their photodestruction rates of ∼10−5 s−1 and the outflow velocity of ≈1 km s−1. The exception is ammonia, which has a particularly rapid photodissociation rate. The middle panel shows ‘daughter’ species, most of which are produced via photodissociation reactions. The slight kinks in the curves at ∼10 km are due to the fact that within this radius the coma is optically thick to UV radiation (cf. Huebner & Giguere 1980). Note that the abundance of HC3N flattens off after 1000 km. This is because HC3N is created by the reaction between C2H2 and CN, and beyond 1000 km the number densities become sufficiently low to ensure that the time-scale for this reaction exceeds the dynamical time-scale.
The bottom panel in Fig. 3 shows the abundances of ions in the coma. The overall ionization level increases steadily throughout the coma, and is roughly proportional to R [i.e. n(e)∼1/R (cf. Wegmann et al. 1987; Körösmezey et al. 1987)], but different ions dominate in different regions of the coma. This is because the ion–molecule chemistry in comets consists mainly of a ‘protonation cascade’, whereby successive proton transfer reactions occur when an ion collides with a molecule with a higher proton affinity; a similar chemistry occurs in interstellar hot cores (Rodgers & Charnley 2001a). Thus in the inner coma accounts for the bulk of the ionization, since ammonia has the largest proton affinity of the parent species. As the coma density decreases, the time-scale for proton transfer reactions to ammonia becomes longer than the dynamical time-scale, so the abundance of flattens off, and protonated methanol becomes the major ion (CH3OH is twice as abundant as NH3). Eventually, H3O+ takes over as the dominant ion, but by this point most of the chemistry in the coma has effectively ceased.
Because H2O ice is the largest reservoir of hydrogen atoms in the nucleus, the D/H ratio in the coma as a whole is determined by the HDO/H2O ratio. Therefore we initially ran the model with a HDO/H2O ratio equal to the observed value of 6×10−4, and no deuterium in the rest of the ice, in order to see how effectively the D in HDO could be transferred into other molecules in the coma. We find that, although a small amount of mixing does occur, principally through deuteron transfer reactions, the D/H ratios in molecules other than water remain extremely low.
Thereafter, we ran the model numerous times, in each case varying the initial D/H ratios in the parent molecules. Fig. 4 plots the deuterium fractionation for a number of important coma molecules, for one such set of initial values. It is clear from the figure that parent molecules such as HDO, C2HD and CH3OD preserve the initial D/H ratios present in the nuclear ice. It also apparent that the fractionation of daughter species (e.g. DC3N) remains constant in the coma. This is because the D/H ratios in these species are determined by the D/H ratios in their parents. For example, cyanoacetylene is formed from the reaction between C2H2 and CN, so DC3N/HC3N=0.5×C2HD/C2H2, where the factor of 0.5 is due to the fact that the D and H each have an equal probability of being incorporated into the cyanoacetylene molecule.
Fig. 4 shows that the only species with changing D/H ratios in the coma are ions, and atomic D and H. The H2DO+/H3O+ ratio is determined by the deuterium exchange reaction (4), so it varies with the gas temperature. The ratio is determined by the D/H ratios in ions which transfer protons or deuterons to NH3. In the inner coma, most of the proton transfer comes from methanol, so is half the CH3OD/CH3OH ratio. Further out, H2DO+ and HDO+ are the main source of protons, so the ratio falls toward the values in these species. Note that the high ratio causes a slight increase in the NH2D/NH3 ratio as these ions recombine. However, this increase is extremely small, and it is likely that the CH3OD fractionation in Fig. 4 is an overestimate (we chose this high value to illustrate the maximum possible effects of a large CH3OD/CH3OH ratio on other species), and so we may conclude that chemical mixing does not significantly alter the initial D/H ratios in parent molecules.
The atomic D/H ratio varies throughout the coma because these atoms are produced by different mechanisms in different regions of the coma. In the very inner coma they are produced mainly via dissociative recombination of NH3D+, but further out when the coma becomes optically thin in the UV they come from photodissociation of HDO. At large R, the atomic D/H ratio gradually increases because the heavier D atoms thermalize more quickly than H atoms, so fewer escape from the coma.
The results shown in Fig. 4 assume statistical branching ratios for photodissociation reactions. As was discussed in Section 2.1, we also ran the model with altered photo-rates, assuming preferential ejection of H rather than D from deuterated molecules. Because the total rate coefficients are the same, the abundances of parent species remain as shown in Fig. 4. For daughter species, the qualitative behaviour is the same – the fractionation is constant at some value determined by the parent – but the D/H ratios are slightly different. For example, when we assume that the O–H bond is three times more likely to break than the O–D bond, the amount of OD in the coma is increased by 50 per cent, but the amount of atomic D is a factor of 2 smaller.
In addition to running the model with different assumed values for the initial molecular D/H ratios and photodissociation branching ratios, we also used different values for the parameters R0, rh and Q0. On every occasion, we found that the behaviour depicted in Fig. 4 remains qualitatively the same; the molecular D/H ratios remain constant throughout the coma. Therefore we conclude that the D/H ratios observed in cometary parent molecules are the same as in the nucleus, and that the D/H ratios in daughter species derive from the ratios in their parent(s). A similar conclusion was reached by Rodgers & Millar (1996) regarding D/H ratios in interstellar hot cores.
Possible formation mechanisms for DNC and HNC
We have shown that the molecular D/H ratios in daughter molecules reflect the D/H ratios in their parents, so it may be possible to use observed D/H ratios to identify the principal production mechanisms for cometary daughter molecules. A prime example is HNC, the origin of which is uncertain despite its detection in comets Hyakutake, Hale–Bopp and Lee (Irvine et al. 1996, 1998; Biver et al. 1999b). Irvine et al. (1998) proposed that HNC is formed from ion–molecule chemistry via protonation of the parent species HCN followed by dissociative recombination. However, in an earlier paper (Rodgers & Charnley 1998) we showed that such reactions cannot form sufficient HNC, and suggested a mechanism whereby HCN is isomerized to HNC by the energetic H atoms in the coma. Since then, however, the large HNC/HCN ratio seen in comet Lee appears to rule this mechanism out also, and it appears that HNC must be a photodissociation product of some unknown parent, possibly large organic polymers (Rodgers & Charnley 2001c). In addition, in our earlier papers we assumed that the mean kinetic energy of the H atoms produced in the photodissociation of water is 1.75 eV, roughly half the total excess energy available to the system (Crovisier 1989). As we discuss in Section A6.2, recent laboratory work suggests that as much as 75 per cent of the available energy ends up in internal energy of OH, which means that fast H atoms may have less energy than we originally assumed. Nevertheless, it is clear that suprathermal reactions are potentially much more important than ion–molecule chemistry, so we consider the possibility that such reactions may lead to DNC.
Fig. 5 shows the predicted DNC/HNC ratios for different formation mechanisms for these species. If HNC and DNC are formed via ion–molecule chemistry, DNC results from deuteron transfer to HCN to form HCND+, followed by proton transfer to NH3. Thus the DNC/HNC ratio reflects the D+/H+ ratio in species that can transfer protons or deuterons to HCN, primarily H2DO+ and H3O+. Line (a) shows the predicted DNC/HNC ratio in this case; it is clear from comparison with Fig. 4 that, as expected, the value is one-third of the H2DO+/H3O+ ratio.
For HNC formed via reactions of HCN with fast supra-thermal hydrogen atoms, the DNC fractionation depends on the precise mechanism of this reaction. Specifically, it depends on whether the collision energy is sufficient to drive isomerization between the two deuterated isomers of the excited intermediary, DCNH* and HCND*. If not, then DNC can only be formed from collisions of fast D atoms with HCN. Lines (b) and (c) in Fig. 5 show the results in this case; line (b) is for statistical photodissociation branching ratios, whereas (c) is for preferential ejection of H atoms. In both cases the predicted DNC/HNC ratio is low, since the fast D/H ratio is equal to 0.5 or 0.25 times the HDO/H2O ratio. However, if fast H atoms are capable of reacting with DCN to form HCND*, then it is clear that the large DCN fractionation of 0.002 can lead to increased DNC fractionation (line d).
Herbst & Talbi (1998) calculated the potential energy surfaces of the HDCN system, and showed that the activation energy to form the HDCN complex from DCN+H or HCN+D is ≈0.35 eV. Although they did not calculate the energy required to form the HCND and DCNH complexes, Talbi, Ellinger & Herbst (1996) showed that the activation energy to form HCNH from HCN or HNC+H is 0.18 or 0.84 eV. If the energy barrier between the deuterated isotopomers is of a similar magnitude, then HCND* and DCNH* formed by fast H atom impacts will have sufficient energy to move freely between the two structures. Therefore line (d) will be appropriate and we predict a large DNC/HNC ratio for cometary DNC and HNC formed via fast particle reactions.
Owing to the fact that HNC may also be coming from dust, it is uncertain whether or not the DNC/HNC ratio could reliably be used to constrain its source. However, it may be possible to rule out ion–molecule reactions and/or fast H reactions if this ratio can be observed. It will also be instructive to compare the DNC/HNC ratio with the overall D/H ratios in refractory cometary organics which will be measured by the Stardust and Rosetta probes.
We have shown that the HDO/H2O and DCN/HCN ratios in cometary comae remain constant, and so represent the values in the molecules sublimating from the nucleus. So it appears that observations of D/H ratios in coma molecules can be directly translated into nuclear values. However, Blake et al. (1999) mapped the HDO and DCN distributions in Hale–Bopp, and showed that both species peaked in abundance ≈3000 km sunward of the nucleus, at the same position as dust jets were observed. They derived D/H ratios at this point that were 5–10 times larger than in the coma as a whole. In order to explain the increased DCN/HCN ratio, Blake et al. (1999) proposed that these molecules were subliming directly from grains, whereas DCN molecules in the nucleus passed through a thin layer of warm liquid water just below the surface and underwent D/H exchange reactions with this water. Thus the DCN/HCN ratio in the dust jets would correspond to the actual value in the nuclear ice, whereas the ‘nuclear’ value would be diluted via reactions with the subsurface liquid water layer. This mechanism has not been verified experimentally and, although it could possibly explain the enhanced DCN/HCN ratio, it cannot explain why the HDO fractionation is also larger in the dust jets.
An alternative explanation is that the DCN and HCN in these regions are coming from the photodestruction of large organic molecules or dust particles. It is probable that large CHON particles are a source of cometary HNC (Rodgers & Charnley 2001c), so it is likely that such particles will also provide a source of HCN, in addition to the nucleus. Indeed, observations indicate an extended source of HCN around 15–20 per cent as strong as the nuclear production rate (Wright et al. 1998; Blake et al. 1999; Veal et al. 2000). In this case, the large DCN fractionation in the dust jets would trace not molecular DCN in the nucleus, but the D/H ratio in its refractory parent. Nevertheless, this still fails to explain the large HDO/H2O ratio in the jets.
Fortunately, because H2O is the primary constituent of the nuclear ice it is also the principal reservoir of cometary deuterium, so even if D/H exchange reactions do occur in a liquid layer they will not be able to alter the HDO/H2O ratio substantially. Therefore we can be confident that the HDO fractionation of 6×10−4 seen in comets Halley, Hyakutake and Hale–Bopp represents the overall nuclear value. As discussed in the Introduction, it is instructive to compare cometary D/H ratios with those in interstellar hot cores, since the values in the latter reflect the values in interstellar ices. HDO has been detected in a number of hot cores (Jacq et al. 1990; Gensheimer, Mauersberger & Wilson 1996), with inferred HDO/H2O ratios of a few times 10−4, similar to the cometary value.
As discussed in the Introduction, DCN has also been seen in hot cores with a fractionation of ∼0.001 (Hatchell et al. 1998), but it is not certain that hot core HCN and DCN are actually coming from grains. A better comparison of cometary and interstellar ices may be obtained if the coma D/H ratios in the prototypical hot core species methanol and ammonia can be obtained. In particular, the CH2DOH/CH3OD and NH3:NH2D:NHD2 ratios can be used to constrain the formation pathways to methanol and ammonia (Charnley, Tielens & Rodgers 1997; Rodgers & Charnley 2001d).
Finally, it is possible to compare the cometary HDO/H2O ratio with that in the Earth's oceans, which is measured to be 3×10−4 (De Witt, van der Straaten & Mook 1980). This is half the cometary value, which casts doubt upon the hypothesis that cometary impacts brought the bulk of terrestrial water to the Earth (e.g. Chyba 1990). Rather, it is more consistent with the oceans being formed from out-gassing of H2O trapped in the planetesimals from which the Earth formed, since the mean D/H ratio in meteorites is similar to the terrestrial HDO/H2O ratio (e.g. Lecluse & Robert 1992; Bockelée-Morvan et al. 1998). On the other hand, Delsemme (2000) argued that the comets that brought volatiles to the Earth would originate primarily from the inner, warmer regions of the protosolar nebula, and so would have reduced D/H ratios; such ‘deuterium-poor’ comets should account for around 4 per cent of the Oort Cloud population. Tentative evidence for this scenario was provided by the anomalous chemical composition of comet LINEAR S4, which was observed to be depleted in a number of molecules relative to other Oort Cloud comets (Mumma et al. 2001; Bockelée-Morvan et al. 2001). Mumma et al. suggested that this could be explained if this comet originated in the proto-Jovian region of the protosolar nebula. In addition, the chemical composition of comet Giacobini–Zinner also showed deviations from ‘typical’ cometary values (Weaver et al. 1999; Mumma et al. 2000). This comet is thought to be a Kuiper Belt object, and so it appears that the chemical composition of comets depends on their place of origin in the protosolar nebula. If this is the case, it is likely that substantial variations in D/H ratios will exist amongst comets.
We have modelled the deuterium chemistry occurring in the comae of comets, following the sublimation of D-rich nuclear ices. Because deuterium fractionation is extremely sensitive to the gas temperature, we have also modelled the hydrodynamics of the outflowing gas. We find that, although chemical reactions can potentially alter the initial D/H ratios via deuteron transfer reactions, the relatively low ionization ensures that very little mixing actually occurs. Thus observations of molecular D/H ratios in cometary comae can be used to infer the ratios in the nuclear ice.
For daughter species, formed in the coma via chemical reactions, the D/H ratios reflect that in their parents. Hence deuterium fractionation can be used as a tracer of molecular formation pathways. If DNC can be detected in comets, this may be extremely useful in determining the origin of cometary HNC, as we have shown that the DNC/HNC ratio will depend on whether or not it can be formed from isomerization of DCN.
Both deuterated species so far detected in comets, HDO and DCN, have fractionations which are similar to those observed in interstellar hot cores. This has led some authors to conclude that comets must consist of relatively pristine interstellar material. However, enhanced molecular D/H ratios are a general feature of low-temperature chemistry, so the cometary ratios may simply reflect the conditions in the outer protosolar nebula where comets were formed. Thus the similarity of the observed ratios is not necessarily an indicator that comets are made of pristine interstellar grains (Bockelée-Morvan et al. 1998, 2000; Irvine et al. 2000).
Models of the deuterium fractionation in the protosolar nebula have in fact shown that D/H ratios can be altered from their original values. These models fall into two categories: those that calculate the D/H ratios via diffusion from the inner nebula where reaction (1) reduces the HDO/H2O ratio (e.g. Drouart et al. 1999; Mousis et al. 2000); and those that follow the ion–molecule chemistry occurring in the nebula (e.g. Aikawa & Herbst 1999, 2001). Future theoretical work should combine these viewpoints, in order that a fully self-consistent description can be obtained.
It is also important to obtain more detections of other deuterated molecules in comets. Measuring the D/H ratios in species such as ammonia and methanol will be extremely useful in comparing cometary and interstellar ices. Also, the three comets in which deuterated molecules have been seen are all Oort Cloud comets, and it will be interesting to see if the fractionation is the same in Kuiper Belt comets. If Oort Cloud and Kuiper Belt comets consist of chemically distinct populations, it is likely that the molecular D/H ratios will be significantly different, since they are so temperature-sensitive. A'Hearn et al. (1995) showed that Kuiper Belt comets are much more likely to be depleted in C2 and C3, and suggested that this is related to a temperature gradient in the protosolar nebula which inhibited carbon-chain formation. More recently, Mumma et al. (2000, 2001) suggested that the anomalous chemical composition of comets LINEAR S4 and Giacobini–Zinner may be related to their place of origin.
Finally, it is worth noting the parallels between deuterium fractionation and ortho–para ratios. Both values can be used to infer the formation temperature of the cometary molecules, and both can be potentially altered in the coma by similar mechanisms. Just like D/H ratios, ortho–para ratios are altered via proton transfer reactions, and are also affected by ion–molecule spin-exchange reactions (Dalgarno, Black & Weisheit 1973). Our results presented in this paper – that chemical reactions are unable to alter significantly the nuclear D/H ratios in the coma – suggest that the ortho–para ratios observed in the coma are equal to the values in the nuclear ice. In order to quantify the degree of modification that can occur in the coma, we are currently modifying our chemical model to calculate ortho–para ratios explicitly.
Theoretical astrochemistry at NASA Ames is supported by NASA's Origins of Solar Systems and Exobiology Programs through NASA Ames Interchange NCC2-1162. SDR is supported by a National Research Council postdoctoral research associateship
Theoretical astrochemistry at NASA Ames is supported by NASA's Origins of Solar Systems and Exobiology Programs through NASA Ames Interchange NCC2-1162. SDR is supported by a National Research Council postdoctoral research associateship.
Appendix A: Description of The Model
The model separates the coma into several distinct constituents: neutral, ion and electron fluids, and fast H, H2, D and HD photo-products. The latter are important since they affect both the hydrodynamics, as they remove energy which would otherwise heat the gas, and the chemistry, as they can drive reactions with large activation energy barriers. As in most previous coma chemistry models, it is assumed that the outflow is steady with spherical symmetry. Owing to the existence of a sonic point in the plasma when a full multifluid treatment is attempted in the steady approximation, all three fluids are assumed to flow outward with the same velocity (Marconi & Mendis 1986; Körösmezey et al. 1987); this issue is discussed in Appendix B. However, the neutrals, ions and electrons have distinct temperatures. The model calculates the abundances of 181 chemical species linked by over 3500 reactions.
The nucleus of comet Halley was far from symmetrical (Keller et al. 1987), and comets are known to experience episodic outbursts. Hence, to model the gas temperatures and densities adequately around rotating nuclei of arbitrary shape represents an extremely large computational challenge. Currently, the most advanced models of complex asymmetric circumnuclear processes are being developed by Crifo & Rodionov (1997a,b, 2000), but we are still a long way from a full understanding of how the properties of the nucleus (e.g. shape, size, rotation period, inclination of axis, etc.) determine the gas flow near the surface. Luckily, it appears that these initial inhomogeneities are smoothed out as the gas expands, and beyond around 1000 km the gaseous coma displays spherical symmetry (Crifo 1991). For a spherically symmetric coma, the flow will be steady as long as there are no major diurnal variations in the gas production rate. Of course, the gas production rate will increase as the comet approaches the Sun, but the time-scale for significant changes in rh is of the order of weeks. In comparison, the time taken for the coma gas to travel a distance of 3×105 km (roughly the distance at which most of the H2O has been photodissociated, and over 10 times the distance at which the number densities become sufficiently low effectively to halt the chemistry) is ∼3 d.
The model also neglects the effects of gravity and magnetic fields. The first of these assumptions is reasonable, since the thermal energy density of the gas at the nuclear surface is over 1000 times the gravitational potential energy density. The neglect of magnetic fields is a more serious shortcoming of the model, as they may have significant effects on the plasma flow in actual comets. Their omission is justified purely on grounds of simplicity – the primary aim of this paper is to follow the coma chemistry, which necessitates certain approximations when calculating the hydrodynamics. Note that neglect of magnetic field effects will have little effect on the accuracy of our modelling of the inner coma, as theoretical models predict that the magnetic field should be very small in this region (Wegmann et al. 1987; Liu 1999). In comet Halley, such a diamagnetic cavity was indeed observed; the magnetic field strength was less than 0.1 nT within 5000 km of the nucleus (Neubauer et al. 1986).
A1 Fluid equationsMarconi & Mendis 1983); but this is sufficiently thin a region that it is reasonable to neglect the gas–dust interaction]. However, the large temperatures attained by the electrons in the outer coma cause the plasma flow to become subsonic at large R (see Appendix B), so we have had to simplify our model slightly in order to remove the resulting singularity from our system of equations.
To circumvent the problem of the plasma sonic transition, we have followed Marconi & Mendis (1986) and assumed that all three fluids flow outward with the same velocity v. In reality, the ion–neutral collisional coupling becomes weaker with increasing cometocentric distance and so, as vi and vn decouple before the outer coma is reached, this leads to ne and Te being overestimated there. However, in a model similar to ours, i.e. using a single fluid velocity and ignoring dust, Marconi & Mendis (1986) were able to reproduce the observed electron density and temperature structure of comet Giacobini–Zinner to within a factor of about 2. Considering all the uncertainties in the chemistry (e.g. rates, branching ratios, etc.), we therefore consider this approximation to be an acceptable one to make.
In our model, a common velocity for all three fluids ensures that our equations remain self-consistent, and that the fundamental conservation principles enshrined in equations (A1)–(A4) are not violated. It also has the benefit that we are still able to calculate individual temperature profiles for the different fluids. Mathematically, this assumption can be incorporated into the above equations by setting F=F′+ƒ and Q=Q′+fv where ƒ represents the drag force acting on each fluid, and F′ and Q′ account for all external momentum and energy sources. Here ƒ does not equal the actual physical drag force that is present, but rather it must be artificially assumed to have the correct value for each fluid to prevent their velocities from diverging. However, note that for the coupling between the ion and electron fluids, ƒ will accurately reflect the Coulomb interaction, since this force is indeed sufficiently strong to ensure that the two fluids always have the same velocity (in this case ƒ=nieE, where e is the electronic charge and E is the electric field strength in the coma, cf. Körösmezey et al. 1987).
The value of γ for water has been measured over a wide range of temperatures (Keenan, Chao & Kaye 1983), and it is found to vary from 1.333 for T≤100 K to 1.252 at T=1000 K. This variation is sufficiently small that we assumed a constant value for γn of 1.3 throughout the coma. An identical value for γi was used, based on the fact that the dominant ions (H3O+, and are non-linear, and so like H2O will possess three rotational degrees of freedom at low temperatures (i.e. γ=4/3), with successively more vibrational degrees of freedom becoming available as the temperature increases. The electrons were assumed to have the classical value of γe=5/3.
If the chemical and hydrodynamic source terms N, M, F and G can be calculated for every species and fluid, and the boundary conditions at the nuclear surface are known, then all the physical and chemical properties of the outflowing gas as a function of R can be calculated by numerical integration of equations (A6), (A7) (one for each fluid) and (A5) (one for each species). This integration was performed using the lsode package, an updated version of gear (Hindmarsh 1983). Therefore the problem comes down to calculating the mass, momentum and energy source terms for each fluid arising from all the relevant physical processes occurring in the coma, and the chemical source terms generated by reactions occurring in the coma.
A2 Source termsSection A3.1.
Note that equation (A6) requires only the total mass source term, Mt, and not the individual terms for each fluid. Therefore, since chemical reactions and scattering processes must conserve mass, they can be neglected. The only process capable of altering the mass of the fluids as a whole is the escape of fast particles from the coma. Thus
Similarly, only the total momentum source term need be calculated, so although inter-fluid scattering and chemical reactions will transfer momentum between fluids, such processes can be ignored when calculating Ft. Again, escape of fast particles is the only mechanism capable of removing momentum from the fluids as a whole, so
The situation regarding the thermal energy source terms is not so simple. Equations (A7) require Gx for each fluid, so all processes that transfer energy between the three fluids must be considered. These include chemical reactions, elastic collisions, loss of fast particles, and inelastic collisions resulting in rovibrational excitation of water molecules which subsequently emit a photon. The thermal energy source terms for each fluid are calculated from
The Gchem terms are due to chemistry, and are found by summing G^x multiplied by the reaction rate ℛ over all reactions, as in equation (A10), where G^x represents the mean thermal energy source term for fluid x in each reaction, and will depend on the type of the reaction, the masses of the reactants and products, the velocities and temperatures of the reactants, and the exo/endo-thermicity of the reaction. In general, calculating G^x for all but the simplest reactions is a non-trivial task. However, with certain simplifying assumptions general formulae can be derived; these are discussed in detail in Section A3.2.
and refer to the heating and cooling caused by elastic collisions of electrons with water molecules and ions. Expressions for these quantities are described below in Section A4. Elastic scattering between neutrals and ions is calculated as part of Gchem (see Section A4).
Ginel represents the energy lost by electrons and gained by the neutral fluid because of inelastic collisions of electrons with water. The neutral fluid will also lose energy as a result of inelastic water–water collisions. In the inner coma, the photons emitted by excited H2O molecules are likely to be reabsorbed by another H2O molecule before they escape from the coma, and the excitation energy will eventually be returned to the gas as heat following super-elastic collisions. Only the fraction of photons that are able to escape from the coma effectively remove energy from the gas; this is represented by the term . The numerical calculation of these terms is described is Section A5.
A3.1 Rate coefficients
All photo-rates are appropriate for a ‘quiet’ Sun, and are scaled by a factor , where rh is the heliocentric distance of the comet in au, and τuv is the UV optical depth in the coma. Strictly speaking, determining τuv requires an iterative calculation, since the value at radius R depends on the column density between R and infinity. However, such a calculation would be prohibitively costly in terms of computer time, so τuv is calculated assuming a 1/R2 density distribution. In this case, the value is given by
We use a value for σuv of 2×10−17 cm2. A more accurate calculation of the photo-rates could be achieved by splitting the solar flux into a number of wavelength bins, and calculating τ for each bin, as has been done by previous authors [e.g. Giguere & Huebner (1978) used five bins; Schmidt et al. (1988) used 120 bins]. Again, however, this would greatly increase the run time of the model, because the rates for every photo-reaction must be recalculated at every value of R by integrating the cross-section over the UV field. Since the values of τ at different UV wavelengths vary by at most a factor of 2 (cf. fig. 2 of Giguere & Huebner 1978), and since the values only climb above unity in the very innermost coma, it was decided to use a single value for τuv. The largest error caused by this assumption will be in the flux of high-energy photons, which have low absorption cross-sections and so will be more abundant in the inner coma than our value of τuv would suggest. Since these photons are only important in determining the rates of photodissociative ionization (PDI) reactions, this limitation is surmounted by assuming that PDI reactions occur at the unattenuated rate throughout the coma.
The most important chemical reactions occurring in the inner coma are exothermic proton transfer reactions, which take place at approximately the collisional rate. If the neutral reactant has a permanent dipole moment, then these rates will depend on the temperature in the coma. The rate coefficient for ion–neutral reactions involving polar molecules can be calculated using the average dipole orientation theory of Su & Bowers (1973). The rate coefficients generated by this theory can be parametrized in the form (Eberhardt & Krankowsky 1995)Flower, Pineau des Forêts & Hartquist 1985): equation (A17) depends on the dipole moment and polarizability of the neutral, and a ‘locking constant’ which measures the degree of alignment between the dipole and the ion. For a value of the locking constant of 0.255, which is within the range of values predicted by numerous theories (see e.g. Chesnavich, Su & Bowers 1978), the constant β can be calculated via
A3.2 Hydrodynamic source terms
In addition to creating and destroying molecules, chemical reactions can add or remove mass, momentum and energy to/from each fluid. As was discussed in Section A2, only the total (over all three fluids) mass and momentum source terms enter the equation for dv/dR, so the contributions of chemical reactions can be ignored, since they act only to transfer these quantities between fluids. However, it is necessary to know the thermal energy source terms for each fluid, and these are generally complex functions of the reaction type and properties of the reactants and products.
The principal obstacle to deriving the mean energy source terms resulting from chemical reactions comes from the dependence of the cross-section on the relative velocity of the reactants. In general, one must calculate the source terms for each pair of reactant velocities, (vA, vB), and then average over all such pairs, weighted according to vAB×σ(vAB), where vAB=|vA−vB|. This requires a double integral over the Maxwellian velocity distributions of A and B. Such integrals do not always have an analytic solution, so it is not possible to write down general formulae for source terms in terms of σ(vAB). Fortunately, the most important reactions in cometary comae are ion–molecule reactions, which typically have cross-sections proportional to . In this case, the double integral is not necessary, since the probability of two molecules A and B reacting is independent of their relative velocity. This means that the energy change of each fluid involved in the reaction depends only the mean properties of A and B averaged over the fluid to which they belong.
Draine (1986) used this fact to derive source terms for a variety of reaction types, and we have extended this analysis to cover all of the reaction types included in our model. The fact that in our model the neutrals and plasma share a common velocity serves to simplify the expressions a great deal. Table A1 lists the mean thermal energy source terms, G^, for the neutrals, ions and electrons for the 12 generic reaction types that we consider. Note that for many reactions the energy change ΔE is unknown. However, the most important reactions in cometary comae are photo-reactions and proton transfer reactions. Average excess energies of the former appropriate for the solar UV field were calculated by Huebner et al. (1992), and the exothermicities of the latter were calculated from the proton affinity data of Hunter & Lias (1998).
A4 Elastic scattering
As there are three fluids, three types of elastic collisions are possible. Ion–neutral elastic scattering will be dominated by collisions involving water molecules, although there will be a contribution from other parent molecules. The most abundant ions in the coma that do not react chemically with H2O, CO or CO2 are H3O+, , HCNH+ and , so each time one of these species collides with an H2O, CO or CO2 molecule it will be elastically scattered. Since such elastic collisions are only a special case of the more general process of ion–neutral reactions, (i.e. type 4 in Table A1 with mA=mC, mB=mD and ΔE=0), these 12 collisions are included in the reaction rate file, and the energy source terms resulting from them are calculated as for all other ion–molecule reactions.
Electron–neutral scattering will also be dominated by collisions with water molecules. Pack, Voshall & Phelps (1962) measured the cross-section for e–H2O collisions; averaged over a Maxwellian velocity distribution this gives a rate coefficient of The mean energy change of the particles in each collision can be obtained from the formula for elastic ion–neutral collisions (see previous paragraph), but replacing mi by me, and recognizing that the mean kinetic energy of a colliding electron is kBTe, not (Draine 1986). This leads to the following expression for the neutral heating caused by electron–water scattering [cf. equations (22) and (24) of Körösmezey et al. (1987)]:
The heat transfer between ions and electrons owing to Coulomb interactions was calculated by Draine (1980), who derived
A5 Inelastic collisions
Water molecules can be collisionally excited into higher rovibrational states, which then decay with the emission of a photon. This process effectively removes kinetic energy from the gas and radiates it out of the coma. There are two main sources of excited H2O: collisions between two water molecules, and collisions between an electron and a water molecule. The first of these removes energy from the neutral fluid, whilst the latter removes it from the electrons.
A semi-empirical formula for the energy loss resulting from water–water collisions was given by Shimizu (1976). When this is adapted to include the effects of radiation trapping in the inner coma (Schmidt et al. 1988), the neutral cooling rate is equal toSchmidt et al. 1988). Crovisier (1984) argued that equation (A22) overestimates the cooling rate by a factor of 3–10, depending on the temperature. Fortunately, the optical depth correction ensures that this cooling mechanism is unimportant in the coma (see Section 3.1), so the uncertainties in this expression will have negligible effects on our solutions.
The energy loss rate of the electrons resulting from inelastic collisions with water was calculated by Cravens & Körösmezey (1986), and the expressions derived by these authors are used in our model. As discussed in the previous paragraph, radiation trapping in the inner coma means that not all of this energy will be radiated away and a proportion will be transferred back to the neutrals via super-elastic water–water collisions. Again, approximating this effect by an optical depth term, we have
A6 Non-thermal photo-products
A6.1 Transition to free molecular flow
As the density of the coma falls off, the use of fluid dynamics to calculate the coma properties becomes less and less realistic. In particular, when the mean free path of a molecule becomes comparable to the size of the coma, it is unreasonable to assume that photodissociation products share their excess energy with other particles. This means that fluid equations can significantly overestimate the temperatures in the outer coma (Ip 1983). A detailed description of the transition region of the coma would require calculating the non-Maxwellian velocity distribution of each species. However, this would make the model far too complex, so a rough approximation has been made, whereby heavy species such as OH and NH2 produced in photodissociation reactions are assumed to be in thermal equilibrium with the rest of the neutrals, but light particles (H, H2, D and HD) produced via the same reactions are treated separately from their thermal counterparts.
Atomic and molecular hydrogen are selected for this special treatment because (i) H and H2 are formed in greater numbers by photodissociation reactions than other species, (ii) most of the excess kinetic energy of such reactions goes to the lighter particle, and (iii) these species require more collisions to be thermalized than do heavier molecules. Thus much of the energy deposited in the coma by photons ends up as kinetic energy of H and H2 molecules. This energy should only be considered to be heating the neutrals as a whole as and when these fast particles are thermalized. The reason why D and HD are also considered is that these energetic particles may have an effect on the coma chemistry. Specifically, they can drive reactions that are endothermic or possess energy barriers too large for their thermal counterparts to overcome.
The net fraction of fast particles that are immediately thermalized (or, conversely, the fraction that escape from the coma without colliding) was calculated by Huebner & Keady (1984). Fig. A1 shows this value as a function of R/Rc, where Rc is the ‘critical radius’ where the mean free path of the fast particles is equal to R. Note that we plot the escape probability, Pesc, whereas Huebner & Keady (1984) plotted current density j1. The two values are related via Pesc=d/dR(R2j1/Rc). Other authors have calculated Pesc using different methods, and although the results differ from model to model, the qualitative agreement between them is very good (Crifo 1991). Huebner & Keady (1984) also calculated the total flux of fast particles as a function of R/Rc. Dividing this by the velocity of these particles gives the number densities, which are necessary to calculate the rates of chemical reactions involving such particles.
A6.2 Energies and mean free paths of fast particles
Early hydrodynamical models assumed that all of the excess energy of photodissociation reactions – of the order of several eV – ended up as kinetic energy of the emitted particles. However, laboratory measurements have shown that a substantial fraction of the energy is often accounted for by ro-vibrational excitation of the emitted radicals. For example, recent work on the photodissociation of H2O at Ly-α has shown that around 75 per cent of the excess energy typically ends up as rotational excitation of the resulting OH, which can possess rotational quantum numbers as high as J=45 (Dixon et al. 1999; Hwang et al. 1999b; Harich et al. 2000). Similar results are also seen for longer wavelengths (Hwang, Yang & Yang 1999a; Zanganeh et al. 2000), and for many other molecules such as NH3, CH4, C2H2, CH3OH and HCN (Mordaunt et al. 1993; Mordaunt, Ashfold & Dixon 1996; Lai, Che & Liu 1996; Harich et al. 1999; Cook et al. 2000).
Crovisier (1989) was the first to point out that this will affect the thermodynamics in the coma, since the excited OH radicals decay rapidly through radiative emission, and so this energy is lost from the coma. The mean total excess energy of reaction (A25) is 3.4 eV (Huebner et al. 1992), and Crovisier (1989) calculated that on average 1.85 eV of that is in the form of translational energy of the products. However, the recent laboratory work discussed above suggests that an even larger fraction may go into internal energy of OH. Thus in our model we assume that 75per cent of the excess energy is lost, and only 0.85 eV goes into kinetic energy. We find that this change has only a very small effect on the coma temperature: comparison of the present Fig. 1 with fig. 1 of Rodgers & Charnley (1998) (where we assumed an energy input of 1.85 eV per photodissociation event) shows that the coma is colder, but only slightly. For reaction (A26) we neglect excitation of the O atom, and assume that all the excess energy (3.84 eV: Huebner et al. 1992) is kinetic.
The energies of these fast H and H2 photoproducts are ≈50–100 times greater than the average thermal energy of a neutral molecule, so each fast particle must undergo many collisions before it is thermalized. Thus an accurate calculation of the energy removed from the coma by these particles requires knowledge of how many collisions each will undergo before it escapes. However, for the sake of simplicity, what is actually calculated is the number of particles that escape with at least half of their initial energy. Although this figure of one-half is arbitrary, it neither drastically over- nor under-estimates the number of fast particles that leave the coma with a substantial fraction of the energy imparted to them by reactions (A25) and (A26).
The amount of energy lost by a fast H or H2 in a collision with a H2O molecule can be estimated by using the same expression as for elastic ion–neutral collisions (see Table A1). It can be shown that the mean number of collisions necessary for fast H and H2 particles to lose half their energy, c¯, is equal to 6 and 3 respectively. Thus a fast particle will travel an average random walk distance equal to before it loses half its energy, where Λcoll is the collisional mean free path, equal to . Treating the H2O molecules as hard spheres of diameter 3.5Å (Allen 1973) implies that σcoll≈10−15 cm2. With the further assumptions that nvR2 is constant, and the fact that the outflow velocity in the outer coma is ≈2.5 times the initial velocity, the critical radius Rc is given byHuebner 1985).
The effects on the hydrodynamic source terms resulting from the escape of fast particles are obtained by multiplying the number of such particles being created, the fraction that escape, and the mass, momentum and energy carried by every such particle. To calculate the latter it is assumed that all fast Hs travel with an rms velocity of 12.4 km s−1, and all fast H2s travel at 17.9 km s−1, corresponding to the excess energies of reactions (A25) and (A26).
Appendix B: Sonic Point in The Plasma
As discussed in Section A1, the Coulomb force ensures that the ions and electrons always have the same velocity. So, if we drop our assumption that the plasma and neutral fluids share a common outflow velocity, we can still use equation (A6) to calculate dvi/dR, but must now take the summation to be over the ion and electron fluids only. Since the plasma sound speed is given by
It is clear from equation (B2) that dvi/dR becomes infinite as vi→ci, and thus it is not possible to numerically integrate the conservation equations across the sonic point. The value of dv/dR can remain finite at the sonic point, but only if ψ→0 as v→c. In this case the sonic transition is smooth, and the value of dv/dR can be obtained by using L'Hôpital's rule, so standard numerical integration can be used to solve the flow equations. Previous work on transonic spherical outflows have focused on sonic points in the neutral fluid caused by gas–dust interactions. Such situations arise in the very inner coma of dusty comets (Marconi & Mendis 1983), and in circumstellar outflows (Tielens 1983; Gail & Sedlmayr 1985). In these cases, the sonic point arises very close to the inner boundary, and these authors were able to make the a priori assumption that the transition is smooth, and then work back from the sonic point to determine the necessary initial conditions for this to occur. Unfortunately, because the sonic transition in cometary plasmas occurs far from the nucleus, this technique is not appropriate. Also, for all sets of realistic initial conditions, we find that ψ at the sonic point cannot be made to approach zero, since it is dominated by the electron heating term. Indeed, the large temperature increase of the electrons is the underlying cause of the sonic transition, since ci increases with Te (equation B1).
The physical explanation for the infinite value of dvi/dR is that the plasma undergoes a shock, which results in a discontinuity in the velocity profile. If such a shock is indeed present in cometary atmospheres then one would expect to see a region of hot, dense, slow-moving ions at around 104 km from the nucleus. In fact, so called ‘ion pile-up’ regions have been seen in a number of comets (Altwegg et al. 1993; Eberhardt & Krankowsky 1995; Bouchez et al. 1998). However, the processes that cause these regions are not well understood, but most probably involve the interaction of cometary plasma with the solar wind and the interplanetary magnetic field (e.g. Israelevich, Ershkovich & Neubauer 1998).
Of course, solving the corresponding system of partial differential equations removes the singularity. However, as these models are extremely computationally intensive, and our primary focus is in the chemical aspects of the innermost coma, we have simply chosen to adapt the flow equations so as to remove the sonic point. As noted in Section A.1, previous models (Marconi & Mendis 1986) that made this approximation have been able to produce models of cometary comae in good agreement with the observations.
This paper has been typeset from a tex/latex file prepared by the author.