Deuterated Molecules

Mon. Not. R. Astron. Soc.
330, 660–674 (2002)
A model of the chemistry in cometary comae: deuterated molecules
S. D. RodgersP and S. B. Charnley

Space Science Division, MS 245-3, NASA Ames Research Center, Moffett Field, CA 94035, USA
Accepted 2001 October 29. Received 2001 October 5; in original form 2001 February 12
A B S T R AC T
Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

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.
Key words: astrochemistry – comets: general – comets: individual: Hale – Bopp – comets:
individual: Hyakutake.
Ehrenfreund & Charnley 2000; Irvine et al. 2000; Bockelée-

1 INTRODUCTION
Morvan et al. 2000). The consensus appears to be that, whilst
Comets are thought to be the most pristine objects in our Solar comets retain a strong signature of their ultimate interstellar origins,
system. As such, their chemical composition reflects both the significant processing must have occurred in the protosolar nebula.
conditions in the outer solar nebula some 4:5 £ 109 yr ago, and the If this is correct, one of the most important effects may be on the
nature of the natal interstellar cloud from which the Solar system deuterium content of the ices. In dark molecular clouds, deuterated
was formed. Assuming that the chemistry in the natal cloud was the molecules are observed with abundances relative to their H-bearing
same as in dark clouds in the local interstellar medium, then analogues of up to 10 per cent (e.g. Howe et al. 1994; Roueff et al.
observations of such clouds tell us the initial chemical composition 2000), over 1000 times the cosmic D/H ratio of ,1025. This is
of the protosolar nebula. Thus a detailed comparison of cometary because replacing an H with a D typically lowers the ground-state
and interstellar ices can be used to constrain models of Solar energy level of a molecule by several hundred kelvin. At low
system formation, since any differences must be due to processes temperatures the ions Hþ þ þ
3 , C2 H2 and CH3 undergo D/H exchange
occurring as the nebula condensed [e.g. many models predict an reactions with HD, becoming highly fractionated (Smith, Adams &
accretion shock sufficiently strong to evaporate ice mantles from Alge 1982). The enhanced deuterium levels in these ions are then
infalling grains (Chick & Cassen 1997; Lunine et al. 1991)]. transferred to other molecules by gas-phase chemical reactions.
Until recently, such comparisons were not possible, since very Models of deuterium chemistry show that the predicted
few data were available. However, the last few years have seen an fractionation is extremely sensitive to the gas temperature (Millar,
explosion in our knowledge of both cometary and interstellar ices. Bennett & Herbst 1989; Rodgers & Millar 1996). Thus, if cometary
The recent apparitions of two very bright comets – Hyakutake in volatiles were indeed evaporated from infalling grains and
1996 and Hale– Bopp in 1997 – led to the detection of many processed before refreezing, it may be expected that the molecular
molecules never before seen in cometary comae (e.g. Mumma et al. D/H ratios would be reduced in the warm gas relative to the
1996; Irvine et al. 1996; Biver et al. 1997, 1999a). At the same original values set in the cold ð10 – 20 KÞ interstellar phase.
time, the launch of the ISO satellite allowed the detailed For example, in the inner region of the protosolar nebula, the
investigation of many interstellar solid state absorption lines (e.g. densities and temperatures were sufficiently large that the
Ehrenfreund 1999; Ehrenfreund & Schutte 2000). Together with deuterium exchange reaction
theoretical and observational progress in the understanding of ‘hot
HDO þ H2 O H2 O þ HD ð1Þ
cores’ (e.g. Charnley 1997; Millar, Macdonald & Gibb 1997b;
Rodgers & Charnley 2001a), this has given us a good reduced the HDO fractionation to the appropriate high-temperature
understanding of the nature of interstellar ices. These new equilibrium value (Geiss & Reeves 1981). Viscous mixing could
measurements have permitted comparisons to be made between have potentially transported material from the inner region to the
interstellar and cometary ices (see e.g. Altwegg et al. 1999; cometary formation zone, thereby diluting the cometary HDO
fractionation (Drouart et al. 1999; Mousis et al. 2000). Yung et al.
P
E-mail: rodgers@dusty.arc.nasa.gov (1988) showed that in the surface layer of the disc, photodestruction
q 2002 RAS
Deuterium chemistry in comets 661
of H2O and HDO followed by the conversion of OH to OD by the give a brief overview. In addition to calculating the chemistry, our
exchange reaction with atomic D, followed by the re-formation of model also follows the dynamical evolution of the gas after it
water, could also alter the HDO/H2O ratio. Again, if turbulent sublimates from the nuclear surface. This is done because previous
mixing were efficient, this process would affect the entire disc. An hydrodynamical models of the coma have shown that the gas
alternative approach was adopted by Aikawa & Herbst (1999, temperatures show strong variations with cometocentric distance
2001), who developed a chemical kinetic model of the deuterium (e.g. Marconi & Mendis 1983, 1986; Crovisier 1984; Huebner
chemistry occurring as matter moves through the disc toward the 1985; Körösmezey et al. 1987; Schmidt et al. 1988; Crifo 1991).
central protostar. They showed that the chemistry in the disc is very The initial expansion is adiabatic so the gas cools, but further out
similar to that in dark clouds, and that ion– molecule reactions in photodissociation and ionization reactions heat the gas. Eventually,
the outer disc may have modified the initial D/H ratios. the temperatures of the different components (ions, electrons and
It is also possible that the observed cometary D/H ratios may be neutrals) can decouple. The electron temperature, in particular, can
modified by gas-phase chemical reactions taking place in the coma. reach extremely high values, .104 K (Eberhardt & Krankowsky
Early models of coma chemistry showed that reactions can occur 1995).
rapidly in the inner coma, where the number densities can be as Many previous chemical models have ignored these effects, and

high as 1012 cm23 (Giguere & Huebner 1978; Huebner & Giguere have simply assumed constant outflow velocities and temperatures
1980; Biermann et al. 1981). The most important reactions were throughout the coma (e.g. Giguere & Huebner 1978; Irvine et al.
found to be proton transfer reactions. These can potentially cycle 1998). Since the gas temperature can have important consequences
deuterium between the different coma molecules, altering the for the chemistry, especially with regard to deuterium fractionation
initial D/H ratios released from the nuclear ice. Thus it is necessary reactions, we have calculated the gas velocity and temperatures in
to construct accurate models of cometary deuterium chemistry, in the coma. Also, it turns out that chemical reactions are an
order that gas-phase coma observations can be safely extrapolated important heat source for the coma gas. As was mentioned in the
to give the nuclear D/H ratios. Introduction, proton transfer reactions – which are typically
To date, only two deuterated molecules have been detected in exothermic by ,1 eV – occur rapidly in the inner coma, and can
cometary comae: HDO and DCN. The former was observed in heat the ion fluid well above what would otherwise be predicted in
comets Halley, Hyakutake and Hale – Bopp (Eberhardt et al. 1995; their absence (see Section 3.1 and Körösmezey et al. 1987).
Bockelée-Morvan et al. 1998; Meier et al. 1998a), and in each Therefore, for any model of the coma to be fully self-consistent, it
comet the HDO/H2O ratio was <6 £ 1024 (this is sometimes must simultaneously solve both the chemical and hydrodynamical
expressed as a D/H ratio of 3 £ 1024 , since the fact that water continuity equations.
contains two hydrogens increases the HDO/H2O over the D/H ratio Having said this, although we include all the major micro-
by a stoichiometric factor of 2). DCN was seen for the first time in physical processes occurring in the coma, and our results are in
Hale– Bopp, with a DCN/HCN ratio of 0.002 (Meier et al. 1998b). good agreement with other hydrodynamical models of comets (see
The fact that the DCN/HCN ratio in interstellar hot cores is also Section 3.1), it must be stressed that the primary aim of our model
,0.002 (Hatchell, Millar & Rodgers 1998) has led some authors to is to follow the coma chemistry. A full multifluid hydrodynamic
conclude that cometary ices must consist of relatively unprocessed treatment requires an iterative calculation (e.g. Körösmezey et al.
interstellar material. However, it is not clear that the DCN and 1987), and the inclusion of magnetic fields, which necessitates
HCN in hot cores are in fact injected directly from ices. Rather, abandoning the assumption of spherical symmetry (e.g. Wegmann
HCN may be formed in the hot gas from reactions initiated by et al. 1987). Both these requirements are not feasible for a chemical
evaporated ammonia or N2 (Rodgers & Charnley 2001a), and the reaction scheme as large as the one in our model. So, although we
DCN/HCN ratio will depend on the relative amounts of N2, NH3 have tried to ensure that our dynamical calculations are as accurate
and NH2D in the gas. Kessler, Qi & Blake (2000) have detected as possible, our model is not intended to be a state of the art
DCN in protoplanetary discs, and derive DCN/HCN ratios of physical model, but rather the hydrodynamical effects are included
,0.01, much larger than the values seen in hot cores and comets, so as to derive realistic number density and temperature profiles in
but consistent with the values seen in cold interstellar gas (Wootten the coma.
1987; Hatchell et al. 1998). Thus, this is further evidence for The chemical reaction scheme contains 181 species and over
processing of D/H ratios in the protosolar nebula. 3500 reactions. It is based on an original scheme developed for
Although only two deuterated cometary molecules have so far modelling interstellar hot cores (Charnley, Tielens & Millar 1992),
been observed, the large number detected in the interstellar with additional reactions taken from the UMIST data base (Millar,
medium (over 15 to date) suggests it is likely that many more will Farquhar & Willacy 1997a). Photodissociation rates appropriate
be detected in the future. Also, the Stardust and Rosetta comet for the solar UV field were taken from Huebner, Keady & Lyon
rendezvous space missions will allow cometary D/H ratios to be (1992). Some additional reactions (principally electron impact
determined via in situ measurements (Wright & Pillinger 1998), reactions) were taken from Schmidt et al. (1988). Also, some
and the FIRST/Herschel mission will allow the HDO/H2O ratio to proton transfer reactions were added that are exothermic but are
be determined in a large sample of comets (Bockelée-Morvan & not included in the UMIST rate file; such reactions have a
Crovisier 2001). We have therefore developed a model of the crucial role in determining the chemical production of molecules
deuterium chemistry occurring in the comae of comets. such as HNC (e.g. Rodgers & Charnley 1998, 2001c) and large
Comparison of our results with future observations will help to O-bearing organics (Rodgers & Charnley 2001a,b; Charnley et al.
elucidate many key questions regarding where, when and how 2001).
comets formed. 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
2 THE MODEL
ratios in the ice were assumed, and the dependence of the results on
A full description of our model appears in Appendix A; here we these values is discussed in Section 3.3. The initial temperature of
q 2002 RAS, MNRAS 330, 660–674

662 S. D. Rodgers and S. B. Charnley
Table 1. Initial abundances assumed for preserve the structure of the molecule. Thus HCND+ recombines to
parent species, relative to water. give HCN or DNC but not DCN or HNC (Millar 1992).
Nevertheless, this leaves many reactions for which the branching
Species Abundance (per cent) Ref. ratios will depend on the preferences of the person or the software
H2 O 100 – generating the deuterated rate file. For example, Markwick,
CO 20 1,2 Charnley & Millar (2001) assume that the reaction of CD+ with
CO2 6 3 CH3OH can form either CH2 DOHþ +
2 or CH3OHD , whereas in our
CH3OH 2 1,2
H2CO† 1 1,2
model only the second route is allowed. Although this is not a
CH4 0.7 4 particularly important reaction, because there are so many
C2H6 0.4 4 reactions in a typical rate file where the final position of the D is
C2H2 0.1 3 ambiguous, there are many possible deuterated versions of a
NH3 1 5 specific rate file. Hence the model results may well depend
HCN 0.1 1,2,6
N2 0.04 7 crucially upon the assumptions that were made when deuterating
the rate file. In the long term, these uncertainties can only be

Notes: † H2CO has an extended dust source resolved through laboratory measurements of deuterated reactions.
in most comets; see e.g. references 1,2 and 8. However, by identifying the key reactions involved in a specific
References: (1) Biver et al. (1999a); (2) Lis
process, and running the model for different values of the
et al. (1997); (3) Bockelée-Morvan et al.
(2000); (4) Mumma et al. (1996); (5) Bird branching ratios in these reactions, one can gain an understanding
et al. (1997); (6) Irvine et al. (1996); (7) of the relative importance of the original assumptions.
Krankowsky (1991); (8) Meier et al. (1993). Photodissociation cross-sections and branching ratios for a few
deuterated molecules have been studied both theoretically and in
the gas was given by (Biver et al. 1997) the laboratory. Zhang & Imre (1988) showed that the cross-section
for photodissociation of HDO is similar to that of H2O, but that
T 0 ¼ 100r 21:2
h ðKÞ; ð2Þ OD þ H production is favoured over OH þ D. The exact value of
where rh is the heliocentric distance in au. The initial H2O number the OD/OH branching ratio will depend on the spectrum of the
density is obtained from incident radiation; where both cross-sections peak ðl <
159 – 172 nmÞ the value is about 2.5 (cf. Bar et al. 1991). At
Q0 ðH2 OÞ longer wavelengths, or if the HDO is initially in a vibrationally
n0 ðH2 OÞ ¼ ; ð3Þ
4pR20 v0 excited state, the OD/OH ratio can be even higher (Zhang, Imre &
Frederick 1989; Vander Wal, Scott & Crim 1990; Vander Wal et al.
where Q0(H2O) is the total water production rate (molecule s21),
1991; Plusquellic, Votava & Nesbitt 1998). Experimental studies of
R0 is the nuclear radius and v0 is the initial gas velocity, assumed to
the photodissociation of deuterated acetylene, methane and
be 1.1 times the sound speed at temperature T0. Different values for
ammonia also show a predominance of H as opposed to D ejection
Q0 and R0, appropriate for different comets and/or heliocentric
(Wang, Hsu & Liu 1997; Wang & Liu 1998; Reid, Loomis & Leone
distances, were used as inputs.
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.
2.1 Deuterium chemistry
Brownsword et al. 1997). Finally, note that Shin et al. (1996)
The original chemical reaction scheme (which comprised some studied the photodissociation of deuterated acetaldehyde ions and
1200 reactions) was extended to include all analogous singly showed that some D/H ‘scrambling’ can occur. Our assumption
deuterated reactions using the techniques described by Rodgers & that deuterated reactions always preserve the integrity of functional
Millar (1996). Essentially, the rules used to generate the deuterated groups is therefore likely to be somewhat naive.
reaction set were as follows. Given the scant data that exist on the photodissociation rates of
deuterated molecules, we have assumed that these reactions occur
(i) The total rate coefficients for deuterated reactions were
at the same rates as their H-bearing analogues. Two different
assumed to be the same as for their H-bearing analogues.
assumptions regarding the branching ratios were made. In the first
(ii) Where more than one hydrogen is involved in a reaction, and
model, these ratios were assumed to be statistical, as for all other
their positions on the product species are not symmetrical, the
reactions involving deuterated species. In a second model, we
branching ratios for the D finishing up in each position are
assumed that cleavage of the X – H bond is three times more likely
proportional to the statistical probabilities (except for some
than for the X – D bond, where X represents O, C or N. There is also
photodissociation and recombination reactions – see below).
evidence for preferential ejection of H atoms in electron
Note that the use of the word ‘involved’ in rule 2 implies that dissociative recombination reactions (Gellene & Porter 1984;
generating a deuterated rate file is necessarily a subjective process. Larsson et al. 1996; Jensen et al. 1999, 2000). For those ions where
Only if the reaction mechanism is known in detail is it possible to the branching ratios have been measured we use the experimental
say which hydrogens take part in the reaction and which are values, but we have continued to assume statistical ejection
passive spectators, otherwise one must make informed guesses. probabilities for all other dissociative recombination reactions.
As in Rodgers & Millar (1996), functional groups present on In contrast to the overall reaction cross-sections, the energetics
both products and reactants were assumed not to be involved in the of several deuterated photodissociation reactions have been well
reaction. For example, photodissociation of CH3OD is assumed to studied in the laboratory. In particular, the H2O/D2O system has
produce CH3 and OD, whereas CH2DOH leads to CH2D and OH, been investigated extensively (e.g. Zanganeh et al. 2000; Harich
in agreement with experimental results (Harich et al. 1999). et al. 2001). These studies demonstrated that the kinetic and
Similarly, dissociative recombination reactions are assumed to internal energies of the emitted D and OD are almost identical to
q 2002 RAS, MNRAS 330, 660–674

those of the H and OH produced in the non-deuterated reaction. It heated mainly by photodissociation reactions. Within ,3000 km,
is important to know the energy distributions of the ejected atoms, all the energy from such reactions is deposited locally in the coma,
since they may be important in driving a supra-thermal chemistry so the heating rate varies as R 22. Beyond this distance, fast H
in the coma (Rodgers & Charnley 1998; see Section A6.2). In atoms and molecules carry off increasing amounts of energy
addition, these studies have also shown that the branching ratios for without being thermalized, leading to a sharper fall-off in the
the three-body channels are similar for D2O and H2O dissociation heating rate.
(Harich et al. 2001). The ions are heated mainly by exothermic proton transfer
Finally, several D/H exchange reactions were included; these are reactions. In the inner coma, each photoionization of H2O to H2O+
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:
H3 Oþ þ HDO $ H2 DOþ þ H2 O: ð4Þ

Reaction (4) was assumed to possess a rate coefficient of
2 £ 1029 cm3 s21 , and to be exothermic to the right by 32 K
(Eberhardt et al. 1995).
3 R E S U LT S
3.1 Hydrodynamics
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 s21 ; 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
Figure 2. Relative importance of different heating and cooling mechanisms.

Solid lines represent heating and dashed lines represent cooling. The total
Figure 1. Radial profiles of gas temperatures and velocity, for a Hyakutake- heating rate for each fluid is shown by the bold lines. Note that
like comet at 1 au. photodissociation reactions are included as part of the ‘Chemistry’ value.
q 2002 RAS, MNRAS 330, 660–674

ultimately results in NHþ 4 (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. Körö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.
3.2 Chemistry
Fig. 3 shows the abundances of a number of important species.
Rather than plotting number densities, the graphs show the value
4pnvR 2, 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 22 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
,1025 s21 and the outflow velocity of <1 km s21. 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; Körösmezey et al. 1987)], but different ions
Figure 3. Molecular fluxes, 4pnvR 2, versus cometocentric distance, R.
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
3.3 Deuterated molecules
an ion collides with a molecule with a higher proton affinity; a
similar chemistry occurs in interstellar hot cores (Rodgers & Because H2O ice is the largest reservoir of hydrogen atoms in the
Charnley 2001a). Thus in the inner coma NHþ 4 accounts for the nucleus, the D/H ratio in the coma as a whole is determined by the
bulk of the ionization, since ammonia has the largest proton affinity HDO/H2O ratio. Therefore we initially ran the model with a
of the parent species. As the coma density decreases, the time-scale HDO/H2O ratio equal to the observed value of 6 £ 1024 , and no
for proton transfer reactions to ammonia becomes longer than the deuterium in the rest of the ice, in order to see how effectively the
dynamical time-scale, so the abundance of NHþ 4 flattens off, and D in HDO could be transferred into other molecules in the coma.
protonated methanol becomes the major ion (CH3OH is twice as We find that, although a small amount of mixing does occur,
abundant as NH3). Eventually, H3O+ takes over as the dominant principally through deuteron transfer reactions, the D/H ratios in
ion, but by this point most of the chemistry in the coma has molecules other than water remain extremely low.
effectively ceased. Thereafter, we ran the model numerous times, in each case
q 2002 RAS, MNRAS 330, 660–674

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
Figure 4. Deuterium fractionation in coma molecules. The values for HDO
D/H ratios remain constant throughout the coma. Therefore we
and DCN correspond to the observed ratios. Values for other parent
conclude that the D/H ratios observed in cometary parent
molecules are arbitrarily chosen. Many model runs were performed using
different input values, and this graph illustrates just one such scenario. For molecules are the same as in the nucleus, and that the D/H ratios
different initial values, the output will obviously differ from that shown in daughter species derive from the ratios in their parent(s). A
here, but we find that the qualitative behaviour of the D/H ratios is always similar conclusion was reached by Rodgers & Millar (1996)
the same (see text). regarding D/H ratios in interstellar hot cores.
3.4 Possible formation mechanisms for DNC and HNC

varying the initial D/H ratios in the parent molecules. Fig. 4 plots
the deuterium fractionation for a number of important coma We have shown that the molecular D/H ratios in daughter
molecules, for one such set of initial values. It is clear from the molecules reflect the D/H ratios in their parents, so it may be
figure that parent molecules such as HDO, C2HD and CH3OD possible to use observed D/H ratios to identify the principal
preserve the initial D/H ratios present in the nuclear ice. It also production mechanisms for cometary daughter molecules. A prime
apparent that the fractionation of daughter species (e.g. DC3N) example is HNC, the origin of which is uncertain despite its
remains constant in the coma. This is because the D/H ratios in detection in comets Hyakutake, Hale– Bopp and Lee (Irvine et al.
these species are determined by the D/H ratios in their parents. For 1996, 1998; Biver et al. 1999b). Irvine et al. (1998) proposed that
example, cyanoacetylene is formed from the reaction between HNC is formed from ion– molecule chemistry via protonation of
C2H2 and CN, so DC3 N=HC3 N ¼ 0:5 £ C2 HD=C2 H2 , where the the parent species HCN followed by dissociative recombination.
factor of 0.5 is due to the fact that the D and H each have an equal However, in an earlier paper (Rodgers & Charnley 1998) we
probability of being incorporated into the cyanoacetylene showed that such reactions cannot form sufficient HNC, and
molecule. suggested a mechanism whereby HCN is isomerized to HNC by
Fig. 4 shows that the only species with changing D/H ratios in the energetic H atoms in the coma. Since then, however, the large
the coma are ions, and atomic D and H. The H2 DOþ /H3 Oþ ratio is HNC/HCN ratio seen in comet Lee appears to rule this mechanism
determined by the deuterium exchange reaction (4), so it varies out also, and it appears that HNC must be a photodissociation
with the gas temperature. The NH3 Dþ /NHþ 4 ratio is determined by product of some unknown parent, possibly large organic polymers
the D/H ratios in ions which transfer protons or deuterons to NH3. (Rodgers & Charnley 2001c). In addition, in our earlier papers we
In the inner coma, most of the proton transfer comes from assumed that the mean kinetic energy of the H atoms produced in
methanol, so NH3 Dþ /NHþ 4 is half the CH3 OD=CH3 OH ratio. the photodissociation of water is 1.75 eV, roughly half the total
Further out, H2DO+ and HDO+ are the main source of protons, so excess energy available to the system (Crovisier 1989). As we
the NH3 Dþ /NHþ 4 ratio falls toward the values in these species. discuss in Section A6.2, recent laboratory work suggests that as
Note that the high NH3 Dþ /NHþ 4 ratio causes a slight increase in the much as 75 per cent of the available energy ends up in internal
NH2 D=NH3 ratio as these ions recombine. However, this increase is energy of OH, which means that fast H atoms may have less energy
extremely small, and it is likely that the CH3OD fractionation in than we originally assumed. Nevertheless, it is clear that
Fig. 4 is an overestimate (we chose this high value to illustrate the suprathermal reactions are potentially much more important than
maximum possible effects of a large CH3 OD=CH3 OH ratio on ion– molecule chemistry, so we consider the possibility that such
other species), and so we may conclude that chemical mixing does reactions may lead to DNC.
not significantly alter the initial D/H ratios in parent molecules. Fig. 5 shows the predicted DNC/HNC ratios for different
The atomic D/H ratio varies throughout the coma because these formation mechanisms for these species. If HNC and DNC are
atoms are produced by different mechanisms in different regions of formed via ion – molecule chemistry, DNC results from deuteron
the coma. In the very inner coma they are produced mainly via transfer to HCN to form HCND+, followed by proton transfer to
dissociative recombination of NH3D+, but further out when the NH3. Thus the DNC/HNC ratio reflects the Dþ /Hþ ratio in species
coma becomes optically thin in the UV they come from that can transfer protons or deuterons to HCN, primarily H2DO+
photodissociation of HDO. At large R, the atomic D/H ratio and H3O+. Line (a) shows the predicted DNC/HNC ratio in this
gradually increases because the heavier D atoms thermalize more case; it is clear from comparison with Fig. 4 that, as expected, the
quickly than H atoms, so fewer escape from the coma. value is one-third of the H2 DOþ /H3 Oþ ratio.
q 2002 RAS, MNRAS 330, 660–674

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,
Figure 5. DNC/HNC ratios in the coma for different production routes to
although it could possibly explain the enhanced DCN/HCN ratio, it
HNC. (a) DNC formed via ion–molecule chemistry; (b) DNC formed via
cannot explain why the HDO fractionation is also larger in the dust
reactions of HCN and fast D atoms, for statistical D production in
photodissociation reactions; (c) as (b), but assuming that H atoms are three jets.
times more likely to be ejected in photodissociation reactions; (d) DNC An alternative explanation is that the DCN and HCN in these
formed via reactions of DCN with fast H atoms in addition to HCN with fast regions are coming from the photodestruction of large organic
D atoms. 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
For HNC formed via reactions of HCN with fast supra-thermal of HCN, in addition to the nucleus. Indeed, observations indicate
hydrogen atoms, the DNC fractionation depends on the precise an extended source of HCN around 15 – 20 per cent as strong as the
mechanism of this reaction. Specifically, it depends on whether the nuclear production rate (Wright et al. 1998; Blake et al. 1999; Veal
collision energy is sufficient to drive isomerization between the et al. 2000). In this case, the large DCN fractionation in the dust
two deuterated isomers of the excited intermediary, DCNH* and jets would trace not molecular DCN in the nucleus, but the D/H
HCND*. If not, then DNC can only be formed from collisions of ratio in its refractory parent. Nevertheless, this still fails to explain
fast D atoms with HCN. Lines (b) and (c) in Fig. 5 show the results the large HDO=H2 O ratio in the jets.
in this case; line (b) is for statistical photodissociation branching Fortunately, because H2O is the primary constituent of the
ratios, whereas (c) is for preferential ejection of H atoms. In both nuclear ice it is also the principal reservoir of cometary deuterium,
cases the predicted DNC/HNC ratio is low, since the fast D/H ratio so even if D/H exchange reactions do occur in a liquid layer they
is equal to 0.5 or 0.25 times the HDO/H2O ratio. However, if fast H will not be able to alter the HDO=H2 O ratio substantially. Therefore
atoms are capable of reacting with DCN to form HCND*, then it is we can be confident that the HDO fractionation of 6 £ 1024 seen in
clear that the large DCN fractionation of 0.002 can lead to comets Halley, Hyakutake and Hale– Bopp represents the overall
increased DNC fractionation (line d). nuclear value. As discussed in the Introduction, it is instructive to
Herbst & Talbi (1998) calculated the potential energy surfaces of compare cometary D/H ratios with those in interstellar hot cores,
the HDCN system, and showed that the activation energy to form since the values in the latter reflect the values in interstellar ices.
the HDCN complex from DCN þ H or HCN þ D is <0.35 eV. HDO has been detected in a number of hot cores (Jacq et al. 1990;
Although they did not calculate the energy required to form the Gensheimer, Mauersberger & Wilson 1996), with inferred
HCND and DCNH complexes, Talbi, Ellinger & Herbst (1996) HDO=H2 O ratios of a few times 1024, similar to the cometary
showed that the activation energy to form HCNH from HCN or value.
HNC þ H is 0.18 or 0.84 eV. If the energy barrier between the As discussed in the Introduction, DCN has also been seen in hot
deuterated isotopomers is of a similar magnitude, then HCND* and cores with a fractionation of ,0.001 (Hatchell et al. 1998), but it is
DCNH* formed by fast H atom impacts will have sufficient energy not certain that hot core HCN and DCN are actually coming from
to move freely between the two structures. Therefore line (d) will grains. A better comparison of cometary and interstellar ices may
be appropriate and we predict a large DNC/HNC ratio for cometary be obtained if the coma D/H ratios in the prototypical hot core
DNC and HNC formed via fast particle reactions. species methanol and ammonia can be obtained. In particular, the
Owing to the fact that HNC may also be coming from dust, it is CH2 DOH=CH3 OD and NH3 : NH2 D : NHD2 ratios can be used to
uncertain whether or not the DNC/HNC ratio could reliably be constrain the formation pathways to methanol and ammonia
used to constrain its source. However, it may be possible to rule out (Charnley, Tielens & Rodgers 1997; Rodgers & Charnley 2001d).
ion– molecule reactions and/or fast H reactions if this ratio can be Finally, it is possible to compare the cometary HDO=H2 O ratio
observed. It will also be instructive to compare the DNC/HNC ratio with that in the Earth’s oceans, which is measured to be 3 £ 1024
with the overall D/H ratios in refractory cometary organics which (De Witt, van der Straaten & Mook 1980). This is half the cometary
will be measured by the Stardust and Rosetta probes. 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
4 DISCUSSION
formed from out-gassing of H2O trapped in the planetesimals from
We have shown that the HDO=H2 O and DCN/HCN ratios in which the Earth formed, since the mean D/H ratio in meteorites is
q 2002 RAS, MNRAS 330, 660–674

similar to the terrestrial HDO=H2 O ratio (e.g. Lecluse & Robert deuterated molecules have been seen are all Oort Cloud comets,
1992; Bockelée-Morvan et al. 1998). On the other hand, Delsemme and it will be interesting to see if the fractionation is the same in
(2000) argued that the comets that brought volatiles to the Earth Kuiper Belt comets. If Oort Cloud and Kuiper Belt comets consist
would originate primarily from the inner, warmer regions of the of chemically distinct populations, it is likely that the molecular
protosolar nebula, and so would have reduced D/H ratios; such D/H ratios will be significantly different, since they are so
‘deuterium-poor’ comets should account for around 4 per cent of temperature-sensitive. A’Hearn et al. (1995) showed that Kuiper
the Oort Cloud population. Tentative evidence for this scenario was Belt comets are much more likely to be depleted in C2 and C3, and
provided by the anomalous chemical composition of comet suggested that this is related to a temperature gradient in the
LINEAR S4, which was observed to be depleted in a number of protosolar nebula which inhibited carbon-chain formation. More
molecules relative to other Oort Cloud comets (Mumma et al. recently, Mumma et al. (2000, 2001) suggested that the anomalous
2001; Bockelée-Morvan et al. 2001). Mumma et al. suggested that chemical composition of comets LINEAR S4 and Giacobini–
this could be explained if this comet originated in the proto-Jovian Zinner may be related to their place of origin.
region of the protosolar nebula. In addition, the chemical Finally, it is worth noting the parallels between deuterium
composition of comet Giacobini– Zinner also showed deviations fractionation and ortho– para ratios. Both values can be used to

from ‘typical’ cometary values (Weaver et al. 1999; Mumma et al. infer the formation temperature of the cometary molecules, and
2000). This comet is thought to be a Kuiper Belt object, and so it both can be potentially altered in the coma by similar mechanisms.
appears that the chemical composition of comets depends on their Just like D/H ratios, ortho– para ratios are altered via proton
place of origin in the protosolar nebula. If this is the case, it is likely transfer reactions, and are also affected by ion – molecule spin-
that substantial variations in D/H ratios will exist amongst comets. 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
5 CONCLUSIONS
that the ortho– para ratios observed in the coma are equal to the
We have modelled the deuterium chemistry occurring in the comae values in the nuclear ice. In order to quantify the degree of
of comets, following the sublimation of D-rich nuclear ices. modification that can occur in the coma, we are currently
Because deuterium fractionation is extremely sensitive to the gas modifying our chemical model to calculate ortho– para ratios
temperature, we have also modelled the hydrodynamics of the explicitly.
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 AC K N O W L E D G M E N T S
mixing actually occurs. Thus observations of molecular D/H ratios Theoretical astrochemistry at NASA Ames is supported by
in cometary comae can be used to infer the ratios in the nuclear ice. NASA’s Origins of Solar Systems and Exobiology Programs
For daughter species, formed in the coma via chemical reactions, through NASA Ames Interchange NCC2-1162. SDR is supported
the D/H ratios reflect that in their parents. Hence deuterium by a National Research Council postdoctoral research
fractionation can be used as a tracer of molecular formation associateship.
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 REFERENCES
it can be formed from isomerization of DCN. A’Hearn M. F., Millis R. L., Schleicher D. G., Osip D. J., Birch P. V., 1995,
Both deuterated species so far detected in comets, HDO and Icarus, 118, 223
DCN, have fractionations which are similar to those observed in Aikawa Y., Herbst E., 1999, ApJ, 526, 314
interstellar hot cores. This has led some authors to conclude that Aikawa Y., Herbst E., 2001, A&A, 371, 1107
comets must consist of relatively pristine interstellar material. Allen C. W., 1973, Astrophysical Quantities. Athlone Press, London
However, enhanced molecular D/H ratios are a general feature of Altwegg K. et al., 1993, A&A, 279, 260
low-temperature chemistry, so the cometary ratios may simply Altwegg K., Ehrenfreund P., Geiss J., Huebner W. F., Levasseur-Regourd
reflect the conditions in the outer protosolar nebula where comets A.-C., 1999, Space Sci. Rev., 90, 373
were formed. Thus the similarity of the observed ratios is not Bar I., Cohen Y., David D., Arusi-Parpar T., Rosenwaks S., Valentini J. J.,
1991, J. Phys. Chem., 95, 3341
necessarily an indicator that comets are made of pristine
Biermann L., Giguere P. T., Huebner W. F., 1982, A&A, 108, 221
interstellar grains (Bockelée-Morvan et al. 1998, 2000; Irvine Bird M. K., Huchtmeier W. K., Gensheimer P., Wilson T. L., Janardhan P.,
et al. 2000). Lemme C., 1997, A&A, 325, L5
Models of the deuterium fractionation in the protosolar nebula Biver N. et al., 1997, Sci, 275, 1915
have in fact shown that D/H ratios can be altered from their original Biver N. et al., 1999a, AJ, 118, 1850
values. These models fall into two categories: those that calculate Biver N. et al., 1999b, IAU Circ. 7203
the D/H ratios via diffusion from the inner nebula where reaction Blake G. A., Qi C., Hogerheijde M. R., Gurwell M. A., Muhleman D. O.,
(1) reduces the HDO=H2 O ratio (e.g. Drouart et al. 1999; Mousis 1999, Nat, 398, 213
et al. 2000); and those that follow the ion – molecule chemistry Bockelée-Morvan D., Crovisier J., 2001, in Pilbratt G. L., Cernicharo J.,
occurring in the nebula (e.g. Aikawa & Herbst 1999, 2001). Future Heras A. M., Prusti T., Harris R., eds, The Promise of First, ESA SP-
460, in press
theoretical work should combine these viewpoints, in order that a
Bockelée-Morvan D. et al., 1998, Icarus, 133, 147
fully self-consistent description can be obtained. Bockelée-Morvan D. et al., 2000, A&A, 353, 1101
It is also important to obtain more detections of other deuterated Bockelée-Morvan D. et al., 2001, Sci, 292, 1339
molecules in comets. Measuring the D/H ratios in species such as Bouchez A. H., Brown M. E., Spinrad H., Misch A., 1998, Icarus, 137, 62
ammonia and methanol will be extremely useful in comparing Brownsword R. A., Laurent T., Hillenkamp M., Vatsa R. K., Volpp H.-R.,
cometary and interstellar ices. Also, the three comets in which 1997, J. Chem. Phys., 106, 9563
q 2002 RAS, MNRAS 330, 660–674

Charnley S. B., 1997, ApJ, 481, 396 Irvine W. M., Bergin E. A., Dickens J. E., Jewitt D., Lovell A. J., Matthews
Charnley S. B., Tielens A. G. G. M., Millar T. J., 1992, ApJ, 399, L71 H. E., Schloerb F. P., Senay M., 1998, Nat, 393, 547
Charnley S. B., Tielens A. G. G. M., Rodgers S. D., 1997, ApJ, 482, L203 Irvine W. M., Schloerb F. P., Crovisier J., Fegley B., Mumma M. J., 2000, in
Charnley S. B., Rodgers S. D., Kuan Y.-J., Huang H.-C., 2001, Adv. Space Mannings V., Boss A. P., Russel S. S., eds, Protostars and Planets IV.
Res., in press Univ. Arizona Press, Tucson, p. 1159
Chesnavich W. J., Su T., Bowers M. T., 1978, in Ausloos P., ed., Kinetics of Israelevich P. L., Ershkovich A. I., Neubauer F. M., 1998, A&A, 329, 765
Ion–Molecule Reactions. NATO, La Baule, p. 152 Jacq T., Walmsley C. M., Henkel C., Baudry A., Mauersberger R., Jewell
Chick K. M., Cassen P., 1997, ApJ, 477, 398 P. R., 1990, A&A, 228, 447
Chyba C. F., 1990, Nat, 343, 129 Jensen M. J., Bilodeau R. C., Heber O., Pedersen H. B., Safvan C. P., Urbain
Cook P. A., Langford S. R., Ashfold M. N. R., Dixon R. N., 2000, J. Chem. X., Zajfman D., Andersen L. H., 1999, Phys. Rev. A, 60, 2970
Phys., 113, 994 Jensen M. J., Bilodeau R. C., Safvan C. P., Seirsen K., Andersen L. H.,
Cravens T. E., Körösmezey A., 1986, Planet. Space Sci., 34, 961 Pedersen H. B., Heber O., 2000, ApJ, 543, 764
Crifo J. F., 1991, in Newburn R. L., Neugebauer M., eds, Comets in the Keenan J., Chao J., Kaye J., 1983, Gas Tables. John Wiley, New York
Post-Halley Era. Kluwer, Dordrecht, p. 937 Keller H. U. et al., 1987, A&A, 187, 807
Crifo J. F., Rodionov A. V., 1997a, Icarus, 127, 319 Kessler J. E., Qi C., Blake G. A., 2000, AAS DPS Meeting No. 32, 65.30

Crifo J. F., Rodionov A. V., 1997b, Icarus, 129, 72 Körösmezey A. et al., 1987, J. Geophys. Res., 92, 7331
Crifo J. F., Rodionov A. V., 2000, Icarus, 148, 464 Krankowsky D., 1991, in Newburn R. L., Neugebauer M. et al., eds, Comets
Crovisier J., 1984, A&A, 130, 361 in the Post-Halley Era. Kluwer, Dordrecht, p. 855
Crovisier J., 1989, A&A, 213, 459 Lai L.-H., Che D.-C., Liu K., 1996, J. Phys. Chem., 100, 6376
Dalgarno A., Black J. H., Weisheit J. C., 1973, Astrophys. Lett., 14, 77 Larsson M. et al., 1996, A&A, 309, L1
Delsemme A. H., 2000, Icarus, 146, 313 Lecluse C., Robert F., 1992, Meteoritics, 27, 248
De Witt J. P., van der Straaten C. M., Mook W. G., 1980, Geostandard. Lis D. C. et al., 1997, Icarus, 130, 355
News., 4, 33 Lisse C. M. et al., 1999, Icarus, 140, 189
Dixon R. N., Hwang D. W., Yang X. F., Harich S., Lin J. J., Yang X., 1999, Liu L.-Z., 1999, ApJ, 520, 399
Sci, 285, 1249 Lunine J. I., Engel S., Rizk B., Horanyi M., 1991, Icarus, 94, 333
Draine B. T., 1980, ApJ, 241, 1021 Marconi M. L., Mendis D. A., 1983, ApJ, 273, 381
Draine B. T., 1986, MNRAS, 220, 133 Marconi M. L., Mendis D. A., 1986, Geophys. Res. Lett., 13, 405
Drouart A., Dubrulle B., Gautier D., Robert F., 1999, Icarus, 140, 129 Markwick A. J., Charnley S. B., Millar T. J., 2001, A&A, 376, 1054
Eberhardt P., Krankowsky D., 1995, A&A, 295, 795 Meier R., Eberhardt P., Krankowsky D., Hodges R. R., 1993, A&A, 277,
Eberhardt P., Reber M., Krankowsky D., Hodges R. R., 1995, A&A, 302, 677
301 Meier R., Owen T. C., Matthews H. E., Jewitt D. C., Bockelée-Morvan D.,
Ehrenfreund P., 1999, Space Sci. Rev., 90, 233 Biver N., Crovisier J., Gautier D., 1998a, Sci, 279, 842
Ehrenfreund P., Charnley S. B., 2000, ARA&A, 38, 427 Meier R. et al., 1998b, Sci, 279, 1707
Ehrenfreund P., Schutte W. A., 2000, in Minh Y. C., van Dishoeck E. F., eds, Millar T. J., 1992, in Singh P. D., ed., Astrochemistry of Cosmic
IAU Symp. 197, Astrochemistry: From Molecular Clouds to Planetary Phenomena. Kluwer, Dordrecht, p. 211
Systems. Astron. Soc. Pac., San Francisco, p. 135 Millar T. J., Bennett A., Herbst E., 1989, ApJ, 340, 906
Flower D. R., Pineau des Forêts G., Hartquist T., 1985, MNRAS, 268, 775 Millar T. J., Farquhar P. R. A., Willacy K., 1997a, A&AS, 121, 139
Gail H.-P., Sedlmayr E., 1985, A&A, 148, 183 Millar T. J., Macdonald G. H., Gibb A. G., 1997b, A&A, 325, 1163
Geiss J., Reeves H., 1981, A&A, 93, 189 Mordaunt D. H., Lambert I. R., Morley G. P., Ashfold M. N. R., Dixon R. N.,
Gellene G. I., Porter R. F., 1984, J. Phys. Chem., 88, 6680 Western C. M., Schnieder L., Welge K. H., 1993, J. Chem. Phys., 98,
Gensheimer P. D., Mauersberger R., Wilson T. L., 1996, A&A, 314, 281 2054
Giguere P. T., Huebner W. F., 1978, ApJ, 223, 638 Mordaunt D. H., Ashfold N. R., Dixon R. N., 1996, J. Chem. Phys., 104,
Harich S. A., Lin J. J., Lee Y. T., Yang X., 1999, J. Phys. Chem. A, 103, 6460
10 324 Mousis O., Gautier D., Bockelée-Morvan D., Robert F., Dubrulle B.,
Harich S. A., Hwang D. W. H., Yang X., Lin J. J., Yang X., Dixon R. N., Drouart A., 2000, Icarus, 148, 513
2000, J. Chem. Phys., 113, 10 073 Mumma M. J., DiSanti M. A., Dello Russo N., Fomenkova M., Magee-
Harich S. A., Yang X., Hwang D. W. H., Lin J. J., Yang X., Dixon R. N., Sauer K., Kaminski C. D., Xie D. X., 1996, Sci, 272, 1310
2001, J. Chem. Phys., 114, 7830 Mumma M. J., DiSanti M. A., Dello Russo N., Magee-Sauer K., Rettig
Hatchell J., Millar T. J., Rodgers S. D., 1998, A&A, 332, 695 T. W., 2000, ApJ, 531, L155
Herbst E., Talbi D., 1998, A&A, 338, 1080 Mumma M. J. et al., 2001, Sci, 292, 1334
Hindmarsh A. C., 1983, in Stepleman R. S., ed., Scientific Computing. Neubauer F. M. et al., 1986, Nat, 321, 352
North-Holland, Amsterdam, p. 55 Pack J. L., Voshall R. E., Phelps A. V., 1962, Phys. Rev., 127, 2084
Howe D. A., Millar T. J., Schilke P., Walmsley C. M., 1994, MNRAS, 267, Plusquellic D. F., Votava O., Nesbitt D. J., 1998, J. Chem. Phys., 109, 6631
59 Reid J. P., Loomis R. A., Leone S. R., 2000a, J. Chem. Phys., 112, 3181
Huebner W. F., 1985, in Diercksen G. H. F., Huebner W. F., Langhoff P. W., Reid J. P., Loomis R. A., Leone S. R., 2000b, J. Phys. Chem. A, 104, 10 139
eds, Molecular Astrophysics. Reidel, Dordrecht, p. 311 Rodgers S. D., Charnley S. B., 1998, ApJ, 501, L227
Huebner W. F., Giguere P. T., 1980, ApJ, 238, 753 Rodgers S. D., Charnley S. B., 2001a, ApJ, 546, 324
Huebner W. F., Keady J. J., 1984, A&A, 135, 177 Rodgers S. D., Charnley S. B., 2001b, MNRAS, 320, L61
Huebner W. F., Keady J. J., Lyon S. P., 1992, Ap&SS, 195, 1 Rodgers S. D., Charnley S. B., 2001c, MNRAS, 323, 84
Hunter E. P., Lias S. G., 1998, in Mallard W. G., Linstrom P. J., eds, NIST Rodgers S. D., Charnley S. B., 2001d, ApJ, 553, 613
Chemistry WebBook, NIST Standard Reference Database No. 69. Rodgers S. D., Millar T. J., 1996, MNRAS, 280, 1046
NIST, Gaithersburg, MD (http://webbook.nist.gov), p. 20 899 Roueff E., Tiné S., Coudert L. H., Pineau des Forêts G., Falgarone E., Gerin
Hwang D. H., Yang X. F., Yang X., 1999a, J. Chem. Phys., 110, 4119 M., 2000, A&A, 354, L63
Hwang D. H., Yang X. F., Harich S., Lin J. J., Yang X., 1999b, J. Chem. Sarmecanic J., Fomenkova M., Jones B., Lavezzi T., 1997, ApJ, 483, L69
Phys., 110, 4123 Schmidt H. U., Wegmann R., Huebner W. F., Boice D. C., 1988, Comput.
Ip W.-H., 1983, ApJ, 264, 726 Phys. Commun., 49, 17
Irvine W. M. et al., 1996, Nat, 383, 418 Shimizu M., 1976, Ap&SS, 40, 149
q 2002 RAS, MNRAS 330, 660–674

Shin S. K., Kim B., Haldeman J. G., Han S.-J., 1996, J. Phys. Chem., 100, number densities become sufficiently low effectively to halt the
8280 chemistry) is ,3 d.
Smith D., Adams N. G., Alge E., 1982, ApJ, 263, 123 The model also neglects the effects of gravity and magnetic
Su T., Bowers M. T., 1973, J. Chem. Phys., 58, 3027
fields. The first of these assumptions is reasonable, since the
Talbi D., Ellinger Y., Herbst E., 1996, A&A, 314, 688
thermal energy density of the gas at the nuclear surface is over
Tielens A. G. G. M., 1983, ApJ, 271, 702
Vander Wal R. L., Scott J. L., Crim F. F., 1990, J. Chem. Phys., 92, 803 1000 times the gravitational potential energy density. The neglect
Vander Wal R. L., Scott J. L., Crim F. F., Weide K., Schinke R., 1991, of magnetic fields is a more serious shortcoming of the model, as
J. Chem. Phys., 94, 3548 they may have significant effects on the plasma flow in actual
Veal J. M. et al., 2000, AJ, 119, 1498 comets. Their omission is justified purely on grounds of simplicity
Wang J.-H., Liu K., 1998, J. Chem. Phys., 109, 7105 – the primary aim of this paper is to follow the coma chemistry,
Wang J.-H., Hsu Y.-T., Liu K., 1997, J. Phys. Chem. A, 101, 6593 which necessitates certain approximations when calculating the
Weaver H. A. et al., 1999, Icarus, 142, 482 hydrodynamics. Note that neglect of magnetic field effects will
Wegmann R., Schmidt H. U., Huebner W. F., Boice D. C., 1987, A&A, 187, have little effect on the accuracy of our modelling of the inner
339 coma, as theoretical models predict that the magnetic field should

Wootten A., 1987, in Vardya M. S., Tarafdar S. P., eds, Astrochemistry.
be very small in this region (Wegmann et al. 1987; Liu 1999). In
Reidel, Dordrecht, p. 311
Wright I. P., Pillinger C. T., 1998, Adv. Space Res., 21, 1537 comet Halley, such a diamagnetic cavity was indeed observed; the
Wright M. C. H. et al., 1998, AJ, 116, 3018 magnetic field strength was less than 0.1 nT within 5000 km of the
Yung Y. L., Friedl R. R., Pinto J. P., Bayes K. D., Wen J.-S., 1988, Icarus, nucleus (Neubauer et al. 1986).
74, 121
Zanganeh A. H., Fillion J. H., Ruiz J., Castillejo M., Lemaire J. L., A1 Fluid equations
Shafizadeh N., Rostas F., 2000, J. Chem. Phys., 112, 5660
Zhang J., Imre D. G., 1988, Chem. Phys. Lett., 149, 233 In the absence of gravitational and magnetic fields, the steady,
Zhang J., Imre D. G., Frederick J. H., 1989, J. Phys. Chem., 93, 1840 spherically symmetric conservation equations for particle number,
mass, momentum and energy for each fluid are
APPENDIX A: DESCRIPTION OF THE 1 d

ðR 2 nvÞ ¼ N; ðA1Þ
MODEL R 2 dR
1 d
The model separates the coma into several distinct constituents: ðR 2 rvÞ ¼ M; ðA2Þ
neutral, ion and electron fluids, and fast H, H2, D and HD photo- R 2 dR
products. The latter are important since they affect both the 1 d d
ðR 2 rv 2 Þ þ ðnkB TÞ ¼ F; ðA3Þ
hydrodynamics, as they remove energy which would otherwise R 2 dR dR
heat the gas, and the chemistry, as they can drive reactions with 2
1 d 2 v g kB T
large activation energy barriers. As in most previous coma R rv þ ¼ Q; ðA4Þ
R 2 dR 2 g21 m
chemistry models, it is assumed that the outflow is steady with
spherical symmetry. Owing to the existence of a sonic point in the where R is the cometocentric distance, n is the number density, r is
plasma when a full multifluid treatment is attempted in the steady the mass density, v is the velocity, kB is Boltzmann’s constant, T is
approximation, all three fluids are assumed to flow outward with the temperature, g is the ratio of specific heats, and m is the mean
the same velocity (Marconi & Mendis 1986; Körösmezey et al. molecular mass. N, M, F and Q are source terms representing the
1987); this issue is discussed in Appendix B. However, the net generation rate per unit volume of particles, mass, momentum
neutrals, ions and electrons have distinct temperatures. The model and energy respectively. In principle, the above equations can be
calculates the abundances of 181 chemical species linked by over rearranged to give the radial derivatives of v, r and T for each fluid,
3500 reactions. which can then be solved to produce a profile of the coma. For
The nucleus of comet Halley was far from symmetrical (Keller steady flows the resulting equations contain a singularity at the
et al. 1987), and comets are known to experience episodic sonic point where the velocity gradient, dv/dR, becomes infinite.
outbursts. Hence, to model the gas temperatures and densities For the neutral fluid this is not a problem as the flow is everywhere
adequately around rotating nuclei of arbitrary shape represents an supersonic [in fact, when the dynamical effects of dust particles are
extremely large computational challenge. Currently, the most considered, it turns out that the flow in the very innermost region,
advanced models of complex asymmetric circumnuclear processes R # 50 m, is subsonic (Marconi & Mendis 1983); but this is
are being developed by Crifo & Rodionov (1997a,b, 2000), but we sufficiently thin a region that it is reasonable to neglect the gas–
are still a long way from a full understanding of how the properties dust interaction]. However, the large temperatures attained by the
of the nucleus (e.g. shape, size, rotation period, inclination of axis, electrons in the outer coma cause the plasma flow to become
etc.) determine the gas flow near the surface. Luckily, it appears subsonic at large R (see Appendix B), so we have had to simplify
that these initial inhomogeneities are smoothed out as the gas our model slightly in order to remove the resulting singularity from
expands, and beyond around 1000 km the gaseous coma displays our system of equations.
spherical symmetry (Crifo 1991). For a spherically symmetric To circumvent the problem of the plasma sonic transition, we
coma, the flow will be steady as long as there are no major diurnal have followed Marconi & Mendis (1986) and assumed that all
variations in the gas production rate. Of course, the gas production three fluids flow outward with the same velocity v. In reality, the
rate will increase as the comet approaches the Sun, but the time- ion– neutral collisional coupling becomes weaker with increasing
scale for significant changes in rh is of the order of weeks. In cometocentric distance and so, as vi and vn decouple before the
comparison, the time taken for the coma gas to travel a distance of outer coma is reached, this leads to ne and Te being overestimated
3 £ 105 km (roughly the distance at which most of the H2O has there. However, in a model similar to ours, i.e. using a single fluid
been photodissociated, and over 10 times the distance at which the velocity and ignoring dust, Marconi & Mendis (1986) were able to
q 2002 RAS, MNRAS 330, 660–674

reproduce the observed electron density and temperature structure conditions at the nuclear surface are known, then all the physical
of comet Giacobini – Zinner to within a factor of about 2. and chemical properties of the outflowing gas as a function of R can
Considering all the uncertainties in the chemistry (e.g. rates, be calculated by numerical integration of equations (A6), (A7)
branching ratios, etc.), we therefore consider this approximation to (one for each fluid) and (A5) (one for each species). This
be an acceptable one to make. integration was performed using the LSODE package, an updated
In our model, a common velocity for all three fluids ensures that version of GEAR (Hindmarsh 1983). Therefore the problem comes
our equations remain self-consistent, and that the fundamental down to calculating the mass, momentum and energy source terms
conservation principles enshrined in equations (A1) – (A4) are not for each fluid arising from all the relevant physical processes
violated. It also has the benefit that we are still able to calculate occurring in the coma, and the chemical source terms generated by
individual temperature profiles for the different fluids. Mathemat- reactions occurring in the coma.
ically, this assumption can be incorporated into the above equations
by setting F ¼ F 0 þ f and Q ¼ Q0 þ fv where f represents the drag
force acting on each fluid, and F0 and Q0 account for all external A2 Source terms
momentum and energy sources. Here f does not equal the actual The chemical source term for a species A, NA, is found by summing

physical drag force that is present, but rather it must be artificially the creation and destruction of A over all reactions, i.e.
assumed to have the correct value for each fluid to prevent their X
velocities from diverging. However, note that for the coupling NA ¼ N^ A R; ðA10Þ
All reactions
between the ion and electron fluids, f will accurately reflect the
Coulomb interaction, since this force is indeed sufficiently strong where N̂A is the net number of molecules of A created in each
to ensure that the two fluids always have the same velocity (in this reaction, and R is the rate (reactions cm23 s21) of the reaction. The
case f ¼ ni eE, where e is the electronic charge and E is the electric rate coefficients used in the calculation of R are discussed in
field strength in the coma, cf. Körösmezey et al. 1987). Section A3.1.
With this alteration, and the assumption that g is constant, the Note that equation (A6) requires only the total mass source term,
above equations can be manipulated to give the following Mt, and not the individual terms for each fluid. Therefore, since
expressions for the radial derivatives of n, v and T: chemical reactions and scattering processes must conserve mass,
dnA N A nA dv 2nA they can be neglected. The only process capable of altering the
¼ 2 2 ; ðA5Þ mass of the fluids as a whole is the escape of fast particles from the
dR v v dR R
coma. Thus
dv 1
¼ M t ¼ 2M f ; ðA11Þ
dR rt v 2 2 ½gnkB Tt
f
where M is the mass removed by such particles.
2v
£ F t v 2 ½ðg 2 1ÞGt 2 M t v þ ½gnkB Tt ;
2
ðA6Þ Similarly, only the total momentum source term need be
R
calculated, so although inter-fluid scattering and chemical

dT x ðgx 2 1ÞT x Gx 2v dv Nx reactions will transfer momentum between fluids, such processes
¼ 2 2 2 ; ðA7Þ can be ignored when calculating Ft. Again, escape of fast particles
dR v nx kB T x R dR ðgx 2 1Þnx
is the only mechanism capable of removing momentum from the
where the subscripts A, x and t respectively imply a particular fluids as a whole, so
chemical species, a particular fluid and summation over all three
fluids. The energy source term Q has been replaced by the thermal F t ¼ 2F f : ðA12Þ
energy source term G, defined by The situation regarding the thermal energy source terms is not so
simple. Equations (A7) require Gx for each fluid, so all processes
2 Mv :
1 2
G ; Q 2 Fv þ ðA8Þ
that transfer energy between the three fluids must be considered.
Although it is also possible to derive an expression for drx /dR, it is These include chemical reactions, elastic collisions, loss of fast
easier to calculate the mass density from the sum of the densities of particles, and inelastic collisions resulting in rovibrational
individual species: excitation of water molecules which subsequently emit a photon.
X The thermal energy source terms for each fluid are calculated from
rx ¼ mA nA : ðA9Þ
Gn ¼ Gchem n þ Gn 2 Gn ;
þ Gen 2 Grad inel f
A[x n ðA13Þ
The value of g for water has been measured over a wide range of Gi ¼ Gchem
i e ;
2 Gion ðA14Þ
temperatures (Keenan, Chao & Kaye 1983), and it is found to vary
Ge ¼ Gchem
e þ Gion
e 2 Ginel
e 2 Gen : ðA15Þ
from 1.333 for T # 100 K to 1.252 at T ¼ 1000 K. This variation is
chem
sufficiently small that we assumed a constant value for gn of 1.3 The G terms are due to chemistry, and are found by
throughout the coma. An identical value for gi was used, based on summing Ĝx multiplied by the reaction rate R over all reactions, as
the fact that the dominant ions (H3O+, NHþ þ
4 and CH3 OH2 Þ are non- in equation (A10), where Ĝx represents the mean thermal energy
linear, and so like H2O will possess three rotational degrees of source term for fluid x in each reaction, and will depend on the type
freedom at low temperatures (i.e. g ¼ 4=3Þ, with successively more of the reaction, the masses of the reactants and products, the
vibrational degrees of freedom becoming available as the velocities and temperatures of the reactants, and the exo/endo-
temperature increases. The electrons were assumed to have the thermicity of the reaction. In general, calculating Ĝx for all but the
classical value of ge ¼ 5=3. simplest reactions is a non-trivial task. However, with certain
If the chemical and hydrodynamic source terms N, M, F and G simplifying assumptions general formulae can be derived; these are
can be calculated for every species and fluid, and the boundary discussed in detail in Section A3.2.
q 2002 RAS, MNRAS 330, 660–674

Gen and Gion
e refer to the heating and cooling caused by elastic where Tin is an effective collision temperature defined by (Flower,
collisions of electrons with water molecules and ions. Expressions Pineau des Forêts & Hartquist 1985):
for these quantities are described below in Section A4. Elastic
mn T i þ mi T n
scattering between neutrals and ions is calculated as part of G chem T in ¼ ; ðA18Þ
mn þ mi
(see Section A4).
G inel represents the energy lost by electrons and gained by the where mi and mn are the respective masses of the ion and the
neutral fluid because of inelastic collisions of electrons with water. neutral. The constant b in equation (A17) depends on the dipole
The neutral fluid will also lose energy as a result of inelastic moment and polarizability of the neutral, and a ‘locking constant’
water – water collisions. In the inner coma, the photons emitted by which measures the degree of alignment between the dipole and
excited H2O molecules are likely to be reabsorbed by another H2O the ion. For a value of the locking constant of 0.255, which is
molecule before they escape from the coma, and the excitation within the range of values predicted by numerous theories (see e.g.
energy will eventually be returned to the gas as heat following Chesnavich, Su & Bowers 1978), the constant b can be calculated
super-elastic collisions. Only the fraction of photons that are able via
to escape from the coma effectively remove energy from the gas; mD

this is represented by the term Grad
n . The numerical calculation of
b¼ pffiffiffiffiffi ; ðA19Þ
mD þ aP
these terms is described is Section A5.
where mD is the dipole moment in debye and aP is the polarizability
in Å3. This low-temperature rate enhancement is calculated for all
A3 Chemistry
molecules with dipole moments larger than 1 debye.
A3.1 Rate coefficients
All photo-rates are appropriate for a ‘quiet’ Sun, and are scaled by
a factor r 22
h exp½2tuv ðRÞ, where rh is the heliocentric distance of A3.2 Hydrodynamic source terms
the comet in au, and tuv is the UV optical depth in the coma.
Strictly speaking, determining tuv requires an iterative calculation, In addition to creating and destroying molecules, chemical
since the value at radius R depends on the column density between reactions can add or remove mass, momentum and energy to/from
R and infinity. However, such a calculation would be prohibitively each fluid. As was discussed in Section A2, only the total (over
costly in terms of computer time, so tuv is calculated assuming a all three fluids) mass and momentum source terms enter the
1/R 2 density distribution. In this case, the value is given by equation for dv/dR, so the contributions of chemical reactions
ð1 can be ignored, since they act only to transfer these quantities
tl ðRÞ ¼ sl nðR0 Þ dR0 < sl nðRÞR: ðA16Þ between fluids. However, it is necessary to know the thermal
R energy source terms for each fluid, and these are generally complex
We use a value for suv of 2 £ 10217 cm2 . A more accurate functions of the reaction type and properties of the reactants and
calculation of the photo-rates could be achieved by splitting the products.
solar flux into a number of wavelength bins, and calculating t for The principal obstacle to deriving the mean energy source terms
each bin, as has been done by previous authors [e.g. Giguere & resulting from chemical reactions comes from the dependence of
Huebner (1978) used five bins; Schmidt et al. (1988) used the cross-section on the relative velocity of the reactants. In
120 bins]. Again, however, this would greatly increase the run time general, one must calculate the source terms for each pair of
of the model, because the rates for every photo-reaction must be reactant velocities, (vA, vB), and then average over all such pairs,
recalculated at every value of R by integrating the cross-section weighted according to vAB £ sðvAB Þ, where vAB ¼ jvA 2 vB j.
over the UV field. Since the values of t at different UV This requires a double integral over the Maxwellian velocity
wavelengths vary by at most a factor of 2 (cf. fig. 2 of Giguere & distributions of A and B. Such integrals do not always have an
Huebner 1978), and since the values only climb above unity in the analytic solution, so it is not possible to write down general
very innermost coma, it was decided to use a single value for tuv. formulae for source terms in terms of s(vAB). Fortunately, the most
The largest error caused by this assumption will be in the flux of important reactions in cometary comae are ion – molecule
high-energy photons, which have low absorption cross-sections reactions, which typically have cross-sections proportional to
and so will be more abundant in the inner coma than our value of v21
AB . In this case, the double integral is not necessary, since the
tuv would suggest. Since these photons are only important in probability of two molecules A and B reacting is independent of
determining the rates of photodissociative ionization (PDI) their relative velocity. This means that the energy change of each
reactions, this limitation is surmounted by assuming that PDI fluid involved in the reaction depends only the mean properties of
reactions occur at the unattenuated rate throughout the coma. A and B averaged over the fluid to which they belong.
The most important chemical reactions occurring in the inner Draine (1986) used this fact to derive source terms for a
coma are exothermic proton transfer reactions, which take place at variety of reaction types, and we have extended this analysis to
approximately the collisional rate. If the neutral reactant has a cover all of the reaction types included in our model. The fact
permanent dipole moment, then these rates will depend on the that in our model the neutrals and plasma share a common
temperature in the coma. The rate coefficient for ion – neutral velocity serves to simplify the expressions a great deal. Table A1
reactions involving polar molecules can be calculated using the lists the mean thermal energy source terms, Ĝ, for the neutrals,
average dipole orientation theory of Su & Bowers (1973). The rate ions and electrons for the 12 generic reaction types that we
coefficients generated by this theory can be parametrized in the consider. Note that for many reactions the energy change DE is
form (Eberhardt & Krankowsky 1995) unknown. However, the most important reactions in cometary
rffiffiffiffiffiffiffi comae are photo-reactions and proton transfer reactions. Average
300 excess energies of the former appropriate for the solar UV field
kðT in Þ ¼ kð300Þ 1 þ b 21 ; ðA17Þ
T in were calculated by Huebner et al. (1992), and the exothermicities
q 2002 RAS, MNRAS 330, 660–674

Table A1. Thermal energy source terms per reaction, Ĝ, for each of the 12 generic reaction types contained in the model.
Reaction type Example Ĝn Ĝi Ĝe
1. Radiative association of neutrals A þ B ! C þ hn 2Qn 0 0

2. Radiative association of ions & neutrals A þ Bþ ! Cþ þ hn 2Qn mA
mC ðQn 2 Qi Þ 0
3. Neutral–neutral AþB!CþD DE 0 0
mA mD þmB mC mA mD þmB mC
4.† Ion–neutral A þ B þ ! C þ Dþ m2
ðQi 2 Qn Þ þ mmDT DE m2
ðQn 2 Qi Þ þ mmCT DE 0
5.† Radiative recombination Aþ þ e2 ! A þ hn Qi T 2Qi T 2Qe
6.† Dissociative recombination Aþ þ e2 ! B þ C Qi þ Qe þ DE 2Qi 2Qe
7. Photoionization A þ hn ! Aþ þ e2 2Qn Qn DE
8.† Photodissociation of neutrals A þ hn ! B þ C DE 0 0
9.† Photodissociation of ions Aþ þ hn ! B þ Cþ mC mB
mA DE þ mA Qi
mB
mA ðDE 2 Qi Þ 0
10. Photodissociative ionisation (PDI) A þ hn ! B þ Cþ þ e2 mC
2 mA Qn mC
m A Qn DE
11. Electron impact ionization (EII) A þ e2 ! Aþ þ e2 þ e2 2Qn Qn 2DE
12. Electron impact dissociation (EID) A þ e2 ! B þ C þ e2 0 0 2DE

Notes: †The expressions for these reactions were derived by Draine (1986). mT is the total mass of the reactants. DE is the average exo/endo-thermicity of each
reaction. Qn ; 32 kB T n , Qi ; 32 kB T i , Qe ; kB T e .
of the latter were calculated from the proton affinity data of

A5 Inelastic collisions
Hunter & Lias (1998).
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
A4 Elastic scattering gas and radiates it out of the coma. There are two main sources of
excited H2O: collisions between two water molecules, and
As there are three fluids, three types of elastic collisions are
collisions between an electron and a water molecule. The first of
possible. Ion – neutral elastic scattering will be dominated by
these removes energy from the neutral fluid, whilst the latter
collisions involving water molecules, although there will be a
removes it from the electrons.
contribution from other parent molecules. The most abundant ions
A semi-empirical formula for the energy loss resulting from
in the coma that do not react chemically with H2O, CO or CO2 are
water – water collisions was given by Shimizu (1976). When this is
H3O+, CH3 OHþ + þ
2 , HCNH and NH4 , so each time one of these
adapted to include the effects of radiation trapping in the inner
species collides with an H2O, CO or CO2 molecule it will be
coma (Schmidt et al. 1988), the neutral cooling rate is equal to
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 DE ¼ 0Þ, these 12 8:5 £ 10219 T 2n n2w
Grad
n ¼ expð2tir Þ; ðA22Þ
collisions are included in the reaction rate file, and the energy nw þ ð2:7 £ 107 T n Þ
source terms resulting from them are calculated as for all other
ion– molecule reactions. where tir is the optical depth in the infrared, given by
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 tir ðRÞ ¼ sir nðRÞR; ðA23Þ
velocity distribution this gives a rate coefficient of 2:63 £
1025 T 20:5
e cm3 s21 : The mean energy change of the particles in here, sir is the mean infrared absorption cross-section per
each collision can be obtained from the formula for elastic ion– molecule, assumed to be 4 £ 10215 cm2 (Schmidt et al. 1988).
neutral collisions (see previous paragraph), but replacing mi by me, Crovisier (1984) argued that equation (A22) overestimates the
and recognizing that the mean kinetic energy of a colliding electron cooling rate by a factor of 3 – 10, depending on the temperature.
is kBTe, not 32 kB T e (Draine 1986). This leads to the following Fortunately, the optical depth correction ensures that this cooling
expression for the neutral heating caused by electron – water mechanism is unimportant in the coma (see Section 3.1), so the
scattering [cf. equations (22) and (24) of Körösmezey et al. uncertainties in this expression will have negligible effects on our
(1987)]: solutions.
The energy loss rate of the electrons resulting from inelastic
Gen ¼ 1:1 £ 10225 nw ne T 21=2
e ð2T e 2 3T n Þ erg cm23 s21 ; ðA20Þ
collisions with water was calculated by Cravens & Körösmezey
where nw is the number density of H2O. (1986), and the expressions derived by these authors are used in our
The heat transfer between ions and electrons owing to Coulomb model. As discussed in the previous paragraph, radiation trapping
interactions was calculated by Draine (1980), who derived in the inner coma means that not all of this energy will be radiated
0 sffiffiffiffiffi1 away and a proportion will be transferred back to the neutrals via
n2
T 3e A super-elastic water – water collisions. Again, approximating this
Gion
e ¼ 1:37 £ 10 T e ðT i 2 T e Þ ln@1:24 £ 104
242 i 21:5
; effect by an optical depth term, we have
mi ni
ðA21Þ Ginel inel

n ¼ Ge ½1 2 expð2tir Þ; ðA24Þ
where mi is the mean molecular mass of the ions, and all quantities
are in cgs units. where Ginel
e is the energy lost by the electrons.
q 2002 RAS, MNRAS 330, 660–674

(1984). Fig. A1 shows this value as a function of R/Rc , where Rc is
A6 Non-thermal photo-products
the ‘critical radius’ where the mean free path of the fast particles is
A6.1 Transition to free molecular flow equal to R. Note that we plot the escape probability, Pesc, whereas
Huebner & Keady (1984) plotted current density j1. The two values
As the density of the coma falls off, the use of fluid dynamics to
are related via Pesc ¼ d=dRðR 2 j1 /Rc Þ. Other authors have calculated
calculate the coma properties becomes less and less realistic. In
Pesc using different methods, and although the results differ from
particular, when the mean free path of a molecule becomes
model to model, the qualitative agreement between them is very
comparable to the size of the coma, it is unreasonable to assume
good (Crifo 1991). Huebner & Keady (1984) also calculated the
that photodissociation products share their excess energy with
total flux of fast particles as a function of R/Rc . Dividing this by the
other particles. This means that fluid equations can significantly
velocity of these particles gives the number densities, which are
overestimate the temperatures in the outer coma (Ip 1983). A
necessary to calculate the rates of chemical reactions involving
detailed description of the transition region of the coma would
such particles.
require calculating the non-Maxwellian velocity distribution of
each species. However, this would make the model far too
A6.2 Energies and mean free paths of fast particles

complex, so a rough approximation has been made, whereby heavy
species such as OH and NH2 produced in photodissociation The vast majority of H and H2 molecules are formed via
reactions are assumed to be in thermal equilibrium with the rest of photodissociation of water:
the neutrals, but light particles (H, H2, D and HD) produced via the
H2 O þ hn ! OH þ H; ðA25Þ
same reactions are treated separately from their thermal
counterparts. H2 O þ hn ! O þ H2 : ðA26Þ
Atomic and molecular hydrogen are selected for this special
Early hydrodynamical models assumed that all of the excess
treatment because (i) H and H2 are formed in greater numbers by
energy of photodissociation reactions – of the order of several eV
photodissociation reactions than other species, (ii) most of the
– ended up as kinetic energy of the emitted particles. However,
excess kinetic energy of such reactions goes to the lighter particle,
laboratory measurements have shown that a substantial fraction of
and (iii) these species require more collisions to be thermalized
the energy is often accounted for by ro-vibrational excitation of the
than do heavier molecules. Thus much of the energy deposited in
emitted radicals. For example, recent work on the photodissocia-
the coma by photons ends up as kinetic energy of H and H2
tion of H2O at Ly-a has shown that around 75 per cent of the excess
molecules. This energy should only be considered to be heating the
energy typically ends up as rotational excitation of the resulting
neutrals as a whole as and when these fast particles are
OH, which can possess rotational quantum numbers as high as
thermalized. The reason why D and HD are also considered is
J ¼ 45 (Dixon et al. 1999; Hwang et al. 1999b; Harich et al. 2000).
that these energetic particles may have an effect on the coma
Similar results are also seen for longer wavelengths (Hwang, Yang
chemistry. Specifically, they can drive reactions that are
& Yang 1999a; Zanganeh et al. 2000), and for many other
endothermic or possess energy barriers too large for their thermal
molecules such as NH3, CH4, C2H2, CH3OH and HCN (Mordaunt
counterparts to overcome.
et al. 1993; Mordaunt, Ashfold & Dixon 1996; Lai, Che & Liu
The net fraction of fast particles that are immediately
1996; Harich et al. 1999; Cook et al. 2000).
thermalized (or, conversely, the fraction that escape from the
Crovisier (1989) was the first to point out that this will affect the
coma without colliding) was calculated by Huebner & Keady
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
Figure A1. Fraction of fast particles that are thermalized locally (i.e. calculated is the number of particles that escape with at least half of
1 2 Pesc Þ as a function of R/Rc , derived from the first-flight current density their initial energy. Although this figure of one-half is arbitrary, it
calculations of Huebner & Keady (1984). neither drastically over- nor under-estimates the number of fast
q 2002 RAS, MNRAS 330, 660–674

particles that leave the coma with a substantial fraction of the It is clear from equation (B2) that dvi /dR becomes infinite as
energy imparted to them by reactions (A25) and (A26). vi ! ci , and thus it is not possible to numerically integrate the
The amount of energy lost by a fast H or H2 in a collision with a conservation equations across the sonic point. The value of dv/dR
H2O molecule can be estimated by using the same expression as for can remain finite at the sonic point, but only if c ! 0 as v ! c. In
elastic ion– neutral collisions (see Table A1). It can be shown that this case the sonic transition is smooth, and the value of dv/dR
the mean number of collisions necessary for fast H and H2 particles can be obtained by using L’Hôpital’s rule, so standard numerical
to lose half their energy, c¯, is equal to 6 and 3 respectively. Thus a integration can be used to solve the flow equations. Previous
pffiffi particle will travel an average random walk distance equal to
fast work on transonic spherical outflows have focused on sonic points
cLcoll before it loses halfpffiits
ffi energy, where Lcoll is the collisional in the neutral fluid caused by gas – dust interactions. Such situations
mean free path, equal to ð 2nscoll Þ21 . Treating the H2O molecules arise in the very inner coma of dusty comets (Marconi & Mendis
as hard spheres of diameter 3.5Å (Allen 1973) implies that 1983), and in circumstellar outflows (Tielens 1983; Gail &
scoll < 10215 cm2 . With the further assumptions that nvR 2 is Sedlmayr 1985). In these cases, the sonic point arises very close
constant, and the fact that the outflow velocity in the outer coma is to the inner boundary, and these authors were able to make the
<2.5 times the initial velocity, the critical radius Rc is given by a priori assumption that the transition is smooth, and then work

back from the sonic point to determine the necessary initial
n0 R2
Rc ¼ 5 £ 10216 pffiffi0 ; ðA27Þ conditions for this to occur. Unfortunately, because the sonic
c transition in cometary plasmas occurs far from the nucleus, this
where n0 and R0 are the initial gas number density and nuclear technique is not appropriate. Also, for all sets of realistic initial
radius, respectively. For typical cometary gas production rates, Rc conditions, we find that c at the sonic point cannot be made to
equals several thousand kilometres (cf. Huebner 1985). approach zero, since it is dominated by the electron heating term.
The effects on the hydrodynamic source terms resulting from the Indeed, the large temperature increase of the electrons is
escape of fast particles are obtained by multiplying the number of the underlying cause of the sonic transition, since ci increases
such particles being created, the fraction that escape, and the mass, with Te (equation B1).
momentum and energy carried by every such particle. To calculate The physical explanation for the infinite value of dvi /dR is
the latter it is assumed that all fast Hs travel with an rms velocity of that the plasma undergoes a shock, which results in a
12.4 km s21, and all fast H2s travel at 17.9 km s21, corresponding discontinuity in the velocity profile. If such a shock is indeed
to the excess energies of reactions (A25) and (A26). 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
APPENDIX B: SONIC POINT IN THE PLASMA number of comets (Altwegg et al. 1993; Eberhardt & Krankowsky
As discussed in Section A1, the Coulomb force ensures that the 1995; Bouchez et al. 1998). However, the processes that cause
ions and electrons always have the same velocity. So, if we drop these regions are not well understood, but most probably involve
our assumption that the plasma and neutral fluids share a common the interaction of cometary plasma with the solar wind and the
outflow velocity, we can still use equation (A6) to calculate dvi /dR, interplanetary magnetic field (e.g. Israelevich, Ershkovich &
but must now take the summation to be over the ion and electron Neubauer 1998).
fluids only. Since the plasma sound speed is given by Of course, solving the corresponding system of partial
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi differential equations removes the singularity. However, as these
gi Pi þ ge Pe kB ðgi T i þ ge T e Þ models are extremely computationally intensive, and our primary
ci ¼ ¼ ; ðB1Þ focus is in the chemical aspects of the innermost coma, we have
ri mi
simply chosen to adapt the flow equations so as to remove the sonic
the expression for dvi /dR can be rewritten in the form point. As noted in Section A.1, previous models (Marconi &
dvi c Mendis 1986) that made this approximation have been able to
¼ ; ðB2Þ produce models of cometary comae in good agreement with the
dR v2i 2 c2i
observations.
where c is given by

1 2vi ri c2i
c ; ðF i þ F e Þvi 2 ðgi 2 1ÞGi 2 ðge 2 1ÞGe 2 M i v2i þ
ri R
ðB3Þ This paper has been typeset from a TEX/LATEX file prepared by the author.
q 2002 RAS, MNRAS 330, 660–674

Deuterated Molecules

Uploaded by

Copyright:

Available Formats

Deuterated Molecules

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Deuterated Molecules

Uploaded by

Copyright:

Available Formats

Mon. Not. R. Astron. Soc.

330, 660–674 (2002)

A model of the chemistry in cometary comae: deuterated molecules

S. D. RodgersP and S. B. Charnley

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

Ehrenfreund & Charnley 2000; Irvine et al. 2000; Bockelée-

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

Figure 2. Relative importance of different heating and cooling mechanisms.

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

3.4 Possible formation mechanisms for DNC and HNC

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

APPENDIX A: DESCRIPTION OF THE 1 d

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Reaction type Example Ĝn Ĝi Ĝe

1. Radiative association of neutrals A þ B ! C þ hn 2Qn 0 0

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

of the latter were calculated from the proton affinity data of

ðA21Þ Ginel inel

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

Downloaded from https://academic.oup.com/mnras/article/330/3/660/1049200 by guest on 11 July 2021

q 2002 RAS, MNRAS 330, 660–674

You might also like