Contents

• Complete depletion in prestellar cores
• 1 Introduction
• 2 Model
• 3 Results
• 4 Discussion
• 5 Conclusions
• References
• 6 Online Material
• Appendix A: Chemical reactions and rate coefficients
• Appendix B: Ice Mantle Thickness

A&A 418, 1035-1043 (2004)
DOI: 10.1051/0004-6361:20035718

Complete depletion in prestellar cores

C. M. Walmsley¹ - D. R. Flower² - G. Pineau des Forêts^3,4

1 - INAF, Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, 50125 Firenze, Italy
2 - Physics Department, The University, Durham DH1 3LE, UK
3 - IAS, Université de Paris-Sud, 92405 Orsay Cedex, France
4 - LUTH, Observatoire de Paris, 92195 Meudon Cedex, France

Received 20 November 2003 / Accepted 9 February 2004

Abstract
We have carried out calculations of ionization equilibrium and deuterium fractionation for conditions appropriate to a completely depleted, low mass pre-protostellar core, where heavy elements such as C, N, and O have vanished from the gas phase and are incorporated in ice mantles frozen on dust grain surfaces. We put particular emphasis on the interpretation of recent observations of H₂D⁺ towards the centre of the prestellar core L 1544 (Caselli et al. 2003) and also compute the ambipolar diffusion timescale. We consider explicitly the ortho and para forms of H₂, H₃⁺, and H₂D⁺. Our results show that the ionization degree under such conditions depends sensitively on the grain size distribution or, more precisely, on the mean grain surface area per hydrogen nucleus. Depending upon this parameter and upon density, the major ion may be H⁺, H₃⁺, or D₃⁺. We show that the abundance of ortho-H₂D⁺ observed towards L 1544 can be explained satisfactorily in terms of a complete depletion model and that this species is, as a consequence, an important tracer of the kinematics of prestellar cores.

Key words: submillimeter - molecular processes - radio lines: ISM

1 Introduction

Determining the physical structure of pre-protostellar cores is one of the keys to understanding the development of protostars. In the simplest view of the development of such cores, they contract slowly towards some "pivotal state'', subsequent to which dynamical collapse sets in. Clearly, one should attempt to determine the parameters of this pivotal state, since these parameters dictate the subsequent evolution. It seems probable that the pivotal state is marked by high column density and low temperature. Likely candidates can be identified using measurements of dust emission and absorption (e.g. André et al. 2000), which indicate temperatures (of both gas and dust) of about 10 K or below, molecular hydrogen column densities in the range 10²² to 10²³ cm^-2, and sizes of several thousand AU.

An additional result from recent observational studies of prestellar cores is that depletion of molecular species on to dust grain surfaces is also a marker of relatively evolved cores (i.e. close to the pivotal state) with high central densities and column densities (e.g. Caselli et al. 2002; Tafalla et al. 2002). Tafalla et al. show that CO becomes depleted by at least a factor of 10, relative to its canonical abundance, above densities of $5\times 10^4$ cm^-3. CS shows similar behaviour, whereas nitrogen containing species such as NH₃ and N₂H⁺ have abundances which are either constant or (in the case of ammonia) increase somewhat in the high density gas surrounding the dust emission peak. Their interpretation of these observations is that molecular nitrogen (N₂), which is the source of the observed NH₃ and N₂H⁺, is sufficiently volatile to remain in the gas phase at densities of order 10⁵cm^-3 (essentially because of spot heating of the grain surfaces by cosmic rays), whereas CO (the source of carbon for most C-containing species) is not.

What happens at still higher densities? It seems likely that, at sufficiently high densities, all species containing heavy elements will condense out. The difference in sublimation energies for CO and N₂ is not great (e.g. Bergin & Langer 1997) and thus one might expect N₂ and other N-containing species to disappear from the gas phase at densities only somewhat higher (or dust temperature slightly lower) than found for CO. Whilst it is true that the laboratory experiments may not accurately simulate interstellar ice surfaces, it seems probable that N₂ disappears from the gas phase at densities above 10⁶ cm^-3.

In fact, the densities inferred from mm dust continuum emission in several objects are of order 10⁶ cm^-3 and there is some evidence for "holes'' in the N₂H⁺ distribution at densities above $3\times 10^5$ cm^-3. For example, Bergin et al. (2002) find a flattening of the N₂H⁺ column density distribution in B68 which could be interpreted as being due to an absence of N₂ towards the dust peak. Belloche (2002) finds a "N₂H⁺ hole'' around the peak of the Class 0 source IRAM04191, which is most easily interpreted in terms of N₂ condensing out at the highest densities in the prestellar core. Last, but not least, Caselli et al. (2003, in what follows CvTCB) have detected compact emission (2000 AU in size) in the 1₁₀-1₁₁ H₂D⁺ line towards the dust emission peak of L 1544. This latter discovery is significant in view of the fact that CvTCB infer an abundance of order 10^-9 relative to H₂. Estimates of the ionization degree in L 1544 (Caselli et al. 2002) are of order $2\times 10^{-9}$. Taken together, these results imply that H₂D⁺(and hence probably H₃⁺) is a major ion. It seems likely that this can be the case only if ions containing heavy elements, such as N₂H⁺ and HCO⁺  are absent, which in turn can be true only if the heavy element content of the gas has condensed on to grain surfaces. If this is the case, one must rely on species lacking heavy elements (such as H₂D⁺) to trace the kinematics of the high density nucleus of the pre-protostellar core. Also, it becomes relevant to establish the ionization degree and other physical parameters of a region where species containing heavy elements have condensed out.

In this paper, we study some of the consequences of complete depletion in an isolated, low mass prestellar core. In particular, we compute the ionization degree and its dependence upon physical parameters such as the gas density and the grain size. We consider in detail deuterium fractionation in these conditions and try to explain the observed H₂D⁺ (or, more precisely, the ortho-H₂D⁺) abundance in objects such as L 1544. In Sect. 2, we discuss the model assumptions and, in Sect. 3, we summarize the results. Then, in Sect. 4, we consider some of the implications of our calculations and in Sect. 5, we briefly summarize our conclusions. Some preliminary results of our study have been presented elsewhere (Walmsley et al. 2004).

2 Model

Our model derives from the work of Flower et al. (2003) and Flower & Pineau des Forêts (2003). The time-dependent chemical rate equations are solved numerically for a situation in which heavy elements have condensed out on to dust grain surfaces. As in the case of the primordial gas, the only significant gas-phase species are compounds containing the elements H, D and He. An important role is played by the grains which influence the charge distribution (ions, charged grains, and electrons) and enable the formation of H₂ and HD. Apart from these latter two processes, we neglect surface chemistry in the present study. We discuss below our assumptions concerning reaction rate coefficients.

First, we consider briefly what is implied by the statement that "no heavy elements are present in the gas phase''. For the purpose of the present study, this means in essence that ions such as H₃⁺, produced as a consequence of cosmic ray ionization of molecular hydrogen, are destroyed by dissociative recombination or recombination on grain surfaces, rather than by transferring a proton to species such as CO and N₂ which have larger proton affinities than H₂. Assuming a rate coefficient k_e for dissociative recombination of H₃⁺ with electrons of $4\times 10^{-7}$ cm³ s^-1(McCall et al. 2002, extrapolated to 10 K with a T^-0.52dependence: see Appendix A) and a rate coefficient for proton transfer reactions between H₃⁺ and species X such as CO or N₂ of the order of the Langevin value, $\alpha _{\rm L} = 10^{-9}$ cm³ s^-1, the condition for "complete depletion'' is [X] $< k_{e}[\rm {e}]/\alpha_{L}$, where the square brackets denote a fractional abundance; this implies that, for example, the CO abundance [CO] in the "complete depletion limit'' should be less than 400[e], or roughly $4\times 10^{-6}$ in cores such as L 1544, where $\rm [e] \approx 10^{-8}$. Thus, the calculations reported here will be relevant when species such as CO have abundances much less than 10^-6.

In what follows, we report the steady-state conditions predicted by our calculations. The timescales to reach chemical equilibrium are somewhat shorter than the dynamical timescales in cores such as L 1544. The free-fall time is approximately $3\times 10^4$ y for a density of 10⁶ cm^-3; but, when magnetic fields are present, the collapse time is about an order of magnitude larger (see Ciolek & Basu 2000; Aikawa et al. 2003). As shown below, the equilibrium timescale for the entire chemistry is roughly 10⁵ y, although ionization equilibrium (with a characteristic time $1/[k_{e}n\rm (e)]$, which is of the order of a few years) is reached much more rapidly. We conclude that the steady-state assumption is a valid first approximation but that some of our results may need to be reconsidered in the light of their time dependence.

2.1 Grain processes

In general, we have followed the approach and adopted the rate coefficients relating to "large grains'' of Flower & Pineau des Forêts (2003, in what follows FPdesF; see their Appendix A). The following processes affecting the grain charge are taken into account: electron attachment; electron detachment by the H₂ fluorescence photons, generated by collisions with the secondary electrons produced by cosmic ray ionization of hydrogen; attachment of positive ions to neutral grains and recombination of positive ions with negatively charged grains. However, for the present applications, their assumption of a Mathis et al. (1977, hereafter MRN) power law distribution of grain sizes (based on the extinction curve derived for diffuse interstellar gas) seemed inappropriate. We have therefore assumed single-size grains, with a dust-to-gas mass density ratio (0.013) consistent with all elements heavier than helium having condensed out. Moreover, we adopt a grain material density of 2 gm cm^-3 for the core plus "dirty ice'' mantle. We treat the grain radius $a_{\rm g}$ as a parameter of the model.

Recombination of ions with electrons on negatively charged grain surfaces is important because, if the available grain surface area is large enough, this process can dominate the destruction of ions such as H₃⁺. Recombination on grains is certainly the main destruction mechanism for protons, which recombine only slowly (radiatively) with electrons in the gas phase, and hence our results depend on the surface area of negatively charged grains. Accordingly, we have considered carefully the distribution of grain charge, adopting the rate coefficients of FPdesF. We assumed that H₂ and HD form on grain surfaces, at rates which scale with the total grain surface area.

2.2 Chemical scheme and D-fractionation

The assumption of complete depletion of the heavy elements on to grains leads to a considerable reduction in the number of chemical species which have to be taken into account: H, H₂, H⁺, H₂⁺, H₃⁺ (and their deuterated forms), He and He⁺. At the low temperatures which prevail in protostellar cores, deuterium fractionation occurs, resulting in greatly enhanced abundances of H₂D⁺, D₂H⁺ and D₃⁺. In this situation, the rates of the associated forwards and reverse deuterium-substitution reactions depend not only on the differences in zero-point energies along the deuteration sequence (Ramanlal et al. 2003) but also on the ortho/para abundance ratios of H₂, H₂⁺, H₃⁺, and H₂D⁺. Accordingly, we have treated the ortho and para forms of these molecules as separate species in our calculations.

A key role is played by the ortho/para ratio of molecular hydrogen, which is expected to be much higher than in LTE (see the discussion of Le Bourlot 1991); this has consequences for the ortho/para ratios of species such as H₂D⁺ (see Gerlich et al. 2002, hereafter GHR; Pagani et al. 1992), which become enhanced roughly in proportion to that of H₂. Also important in this context are the spin selection rules for ortho-para conversion, discussed, for example, by Uy et al. (1997).

Caution has to be exercised when using the data of Uy et al. (1997). These authors considered conditions in which the ortho and para forms are statistically populated, i.e. in proportion to their nuclear spin statistical weights, (2I + 1). This situation obtains at temperatures and densities much higher than those considered here. In prestellar cores, only the low-lying states of the ortho and para species are significantly populated. Under these circumstances, the statistical weights of the rotational levels, (2J + 1), must be taken explicitly into account. When compiling the rate coefficients for reactions involving the ortho and para forms of H₂, H₂⁺, H₃⁺, and H₂D⁺, we made the simplifying assumption that only the lowest rotational level is populated in each case and treated the corresponding ortho and para forms as distinct species, which interconvert through chemical reactions. The complete set of reactions included in the model and their rate coefficients are given in Appendix A.

In the limit of complete depletion, deuterium fractionation occurs mainly in the reactions

$${\rm H}_{3}^{+} \ + {\rm HD} \ \rightarrow \ {\rm H}_{2}{\rm D}^{+} \ + {\rm H}_{2}$$ (1)

with rate coefficient k₁,

$${\rm H}_{2}{\rm D}^{+} \ + {\rm HD} \ \rightarrow \ {\rm D}_{2}{\rm H}^{+} \ + \ {\rm H}_{2}$$ (2)

with rate coefficient k₂, and

$${\rm D}_{2}{\rm H}^{+} \ + \ {\rm HD} \rightarrow \ {\rm D}_{3}^{+} \ + \ {\rm H}_{2}$$ (3)

with rate coefficient k₃. These reactions are exothermic by roughly 200 K and consequently enhance abundance ratios such as $n({\rm H}_{2}{\rm D}^{+})/n({\rm H}_{3}^{+})$. However, such enhancement can be limited by a variety of processes. If heavy elements are present in the gas phase, the high abundances of molecules with high proton affinities, such as CO, result in the fractionation being transferred to ions such as DCO⁺; this limits the degree of enhancement of deuterium in H₂D⁺. It follows that the depletion of heavy elements favours deuterium fractionation (see Dalgarno & Lepp 1984), and so a "completely depleted core'' is a limiting case where one might expect high degrees of enhancement.

Even in the absence of heavy elements, there are several processes which limit deuterium substitution in species such as H₃⁺. For example, it has been known for a long time that dissociative recombination with electrons, at a rate $k_{e}~n(\rm {e})$, competes with deuterium enrichment in H₂D⁺, for example (Watson 1978). It is also clear that recombination of H₂D⁺, D₂H⁺ and D₃⁺ on grain surfaces can have a similar effect, and so models with a large grain surface area tend to have moderate D-fractionation. The relative importance of these two processes is roughly the same for D-fractionation and for ion neutralization: in models where the D-fractionation is limited by recombination on grain surfaces, the degree of ionization is limited by the same process.

D-fractionation in an ion such as H₂D⁺ is determined also by the transfer of the deuterium to the higher members of the chain, D₂H⁺ and D₃⁺. Thus reactions 2 and 3 restrict the fractionation in H₂D⁺ and D₂H⁺, respectively. As we shall see below, the effect of these reactions is to limit ratios such as [H₂D⁺]/[H₃⁺] and [D₂H⁺]/[H₂D⁺] to values of the order of unity.

GHR pointed out that ortho-H₂ can restrict deuterium fractionation because the reaction between ortho-H₂ and ortho-H₂D⁺ is exothermic and believed to be fast. Depending on internal excitation, this process may also be important for D₂H⁺. We find that the reaction with ortho-H₂ is not usually dominant but is significant because the abundance of ortho-H₂ is high.

The ratio of the H₂D⁺ and H₃⁺ abundances is given approximately by:

$$[{\rm H_2D^+}]/[{\rm H}_3^+] = \frac{2 ~ k_1[{\rm D}]}{k_e[{\rm e}] + k_{ig}[{\rm g}^-] + 2k_2[{\rm D}] + 2k^-_1}$$ (4)

where [D] is the fractional abundance of D in the gas phase, k_ig is the rate coefficient for neutralization of H₂D⁺ on negatively charged grains of number density $n_{{\rm g}^{-}}=[{\rm g}^{-}]~ n_{\rm H}$, and k₁^- is the rate coefficient for the reverse of reaction 1. We see that deuterium fractionation is limited by the terms in the denominator of the right hand side of this equation. One may note also that, whilst [H₂D⁺]/[H₃⁺] and [D₂H⁺]/[H₂D⁺] are limited to k₁/k₂ and k₂/k₃, respectively, there is no such restriction for [D₃⁺]/[D₂H⁺] because D₃⁺ does not react with HD. Thus, in certain circumstances, D₃⁺ can be the principal ion.

2.3 Numerical method

The code described by Flower et al. (2003) has been used to calculate the steady-state abundances of the chemical species, X_i. The program integrates the coupled differential equations

$$\frac{\rm d}{{\rm d}t}~ n({\rm X}_{i}) = f[n({\rm X}_{i}),T,t]$$

for the number densities $n({\rm X}_{i})$, as functions of time, t, until steady state is attained. The kinetic temperature, T, was assumed constant. By adopting this approach, not only are the steady-state values of the number densities obtained for sufficiently large tbut also information is recovered on the timescales required for the various chemical reactions to reach equilibrium.

Figure 1: Sample result for our reference model ( $n(\rm {H}_2) = 10^6$ cm-3, T = 10 K, $a_{\rm g}$ = 0.1 $\mu$m, $\zeta = 3\times 10^{-17}$ s-1) showing the approach to steady-state conditions for ortho/para H2, D and H.

In Fig. 1, we show the results of a sample calculation in which, initially, the hydrogen and deuterium are in molecular form, with an elemental abundance ratio $n_{\rm D}/n_{\rm H} = 1.6 \times 10^{-5}$. The initial ortho/para H₂ ratio was $3.5 \times 10^{-7}$, the equilibrium value at the gas kinetic temperature T = 10 K (assumed constant). We find that the time required to reach steady state for all species is of the order of 10⁵ y (almost independent of number density). However, as noted above, the time to attain ionization equilibrium (a few years) is much shorter than that required to attain equilibrium between ortho- and para-H₂ (the longest chemical timescale). Atomic hydrogen and deuterium are intermediate, reflecting the timescale for the formation of H₂ and HD on grains.

The time characterizing chemical steady state is comparable with the free-fall time, $\tau _{\rm ff}$, at the lowest density ( $n_{\rm H} = 2\times 10^5$ cm^-3) and is an order of magnitude larger than $\tau _{\rm ff}$ at the highest density ($n_{\rm H}$ = $2\times 10^7$ cm^-3) considered. However, magnetic support to the cloud increases the time for core formation and, in practice, the timescale for chemical equilibrium is likely to be comparable with that for contraction and collapse of the protostellar core.

3 Results

3.1 Dependence on model parameters

We present our results with reference to a model with the following parameters: $a_{\rm g}=0.1$ $\mbox{$\mu$ m}$, $\rho_{\rm g}=2$ g cm^-3, T = 10 K, and $\zeta = 3\times 10^{-17}$ s^-1, where $a_{\rm g}$ is the grain radius, $\rho_{\rm g}$ is the mean density of the grain material (core and ice mantle), T is the gas temperature, and $\zeta$ is the rate of cosmic ray ionisation of H₂. To these parameters correspond a fractional abundance of grains $n_{\rm g}/n_{\rm H} = 3.5\times 10^{-12}$ and a grain surface area per hydrogen nucleus $n_{\rm g}~ \sigma _{\rm g}/n_{\rm H} = 1.1\times 10^{-21}$ cm². These numbers have been derived assuming that all heavy elements are incorporated into either the cores or the ice mantles of the grains (corresponding to a dust-to-gas mass density ratio of 0.013).

Figure 2: Abundances of major ions and electrons in our standard model ($a_{\rm g}$ = 0.1 $\mu$m, T = 10 K, $\zeta = 3\times 10^{-17}~\rm s^{-1}$) as functions of the density of molecular hydrogen.

Figure 3: Ortho/para ratios of H2 (multiplied by 104), H3+, and H2D+ as functions of density for our reference model (cf. Fig. 2). The approximate expression of Gerlich et al. (2002) for H2D+ is shown for comparison.

In Fig. 2 are plotted the fractional abundances of the principal ions and of the electrons as functions of $n({\rm H_2})\approx~n_{\rm H}/2$; the results of computations of ortho/para ratios for the same model are presented in Fig. 3 and will be discussed in Sect. 3.4 below. At the lowest density considered, n(H₂) = 10⁵ cm^-3, H⁺ is the major ion and the removal of molecular ions by dissociative recombination is significant. H⁺, on the other hand, is removed by recombination on grain surfaces. As n(H₂) increases, recombination on grains dominates dissociative recombination for the molecular ions. H⁺ becomes a relatively minor species at high density.

The balance between the various deuterated forms of H₃⁺ is determined approximately by Eq. (4) and its analogues for [D₂H⁺]/[H₂D⁺] and [D₃⁺]/[D₂H⁺]. At densities below 10⁶ cm^-3, the dissociative recombination term k_e [e] is dominant in the denominator of the right hand side of Eq. (4), whereas, at higher densities, the term due to the reaction with HD becomes more important and results in [H₂D⁺]/[H₃⁺] (and [D₂H⁺]/[H₂D⁺]) being close to unity. On the other hand, D₃⁺ cannot be destroyed by HD and becomes the dominant ion.

Figure 4: Abundances of major ions and electrons in a model with a smaller grain size of 0.025 $\mu$m.

Figure 4 is analogous to Fig. 2, for an assumed grain radius of 0.025 $\mu$m. In this case, recombination of ions on grain surfaces is the dominant neutralization process for all ions and H⁺ is never dominant. Moreover, for densities below 10⁷ cm^-3, H₃⁺ is more abundant than its deuterated forms. We note that the increased rate of grain neutralization causes the ion abundance to fall and its density dependence to become steeper.

3.2 Ionization degree for complete depletion

The determination of the ionization degree in the region of complete depletion is important for a number of reasons. The degree of ionization determines the coupling of the magnetic field to the gas and hence the timescale for ambipolar diffusion. The electron abundance plays a crucial role also in the chemistry and in deuterium fractionation.

Figure 5: The free electron abundance computed as a function of density for grain sizes between 0.025 and 0.2 $\mu$m. The result of McKee (1989) is shown for comparison (broken curve). Note the decrease in the level of ionization for small grain sizes (large grain surface areas).

In Fig. 5, the free electron abundance, n(e)/ $n(\rm {H}_2)$, computed as a function of density, is compared with the relationship of McKee (1989), $n({\rm e})/n({\rm H}_2) = 1.3\times 10^{-5}n({\rm H}_2)^{-0.5}$, which is often used in discussions of core evolution. It may be seen that our model generally predicts lower degrees of ionization. We find that n(e)/ $n(\rm {H}_2)$varies approximately as $n_{\rm H}^{-0.75}$between $n_{\rm H} =2\times 10^6$ and $2\times 10^7$cm^-3. As already noted, the degree of ionization of the gas is sensitive to the assumed grain size, owing to the importance of the recombination of ions with electrons on grain surfaces. The fraction of electrons attached to grains remains small for all the grain sizes considered here: the ratio of negatively charged grains to free electrons attains 0.20 for for $a_{\rm g} = 0.025$ $\mu$m and $n({\rm H}_{2}) = 10^7$ cm^-3. However, we note that this ratio becomes larger for smaller grain sizes, becoming $\gg$1 for $a_{\rm g} < 0.01$ $\mu$m. The fraction of grains which is neutral decreases from 0.68 for $a_{\rm g} = 0.025$ $\mu$m to 0.20 for $a_{\rm g} = 0.2$ $\mu$m and for $n({\rm H}_{2}) = 10^7$ cm^-3 (but this ratio is insensitive to the gas density). The fraction of positively charged grains is always negligible.

3.3 Ambipolar diffusion in depleted cores

If the core is sub-critical (i.e. if the magnetic field is capable of resisting collapse), the timescale for star formation is determined by that for ambipolar diffusion. Here we estimate the timescale for ambipolar diffusion under the conditions of complete heavy element depletion.

The ambipolar diffusion timescale, $\tau _{\rm AD}$, depends on the ion abundance, n_i/n_n, the rate coefficient for momentum transfer between the ions and the neutrals (essentially H₂), $\langle\sigma v\rangle_{in} \approx 2\times 10^{-9}$ cm³ s^-1, and the masses m_i and m_n of the charged and neutral collision partners. We have taken

$$\tau_{\rm AD} \approx \frac{2}{\pi G~ m_{n}^2} \sum _{i} \fra... ...{n}} \frac{m_{i}m_{n}}{m_{i}+m_{n}}\langle\sigma v\rangle_{in}$$ (5)

where G is the gravitational constant and the summation is over the ionic species. The rate coefficient for momentum transfer approaches its Langevin form

$$\langle\sigma v\rangle_{in}=2\pi e(\alpha /m_{in})^{0.5}$$

at low temperatures (cf. Flower 2000), and hence $\tau _{\rm AD}$ becomes proportional to the square root of the reduced mass, m_in = [m_im_n/(m_i+m_n)]; $\alpha = 7.7\times 10^{-25}$ cm³ is the polarizability of H₂. Thus, $\tau _{\rm AD}$ increases by 50 percent as the dominant ion changes from being H⁺ to being D₃⁺. In Fig. 6, we plot the ratio of the ambipolar diffusion timescale $\tau _{\rm AD}$ to the free-fall timescale $\tau _{\rm ff}$as a function of density. We see from Fig. 6 that the ambipolar diffusion timescale in the limit of complete depletion is of the same order of magnitude as that inferred for normal ionic abundances in dense cores, i.e. about an order of magnitude larger than the free-fall time (see, for example, Shu et al. 1987; Ciolek & Basu 2000). As one might expect, smaller grain sizes lead to shorter ambipolar diffusion timescales, owing to more rapid neutralization on grain surfaces.

Figure 6: Ratio of ambipolar diffusion time $\tau _{\rm AD}$ (Eq. (5)) to free fall time  $\tau _{\rm ff}$ as function of molecular hydrogen density for values of the grain radius between 0.025 and 0.2 $\mu$m.

The above discussion neglects the collisional coupling between the neutral gas and charged grains (cf. Nakano & Umebayashi 1980). For the conditions that we have considered, the coupling to the charged grains is dominant, particularly for small grain sizes and high gas densities. We estimate that its effect is to increase the value of $\tau _{\rm AD}$ by one to two orders of magnitude. We conclude that the collapse of a magnetically sub-critical, completely depleted core is likely to be impossible on reasonable timescales.

3.4 The abundances of ortho-H₂ , ortho-H ₂D ⁺ and ortho-H ₃ ⁺

Figure 3 shows that the ortho/para ratios considered here are relatively insensitive to density. For molecular hydrogen, we derive a steady-state ortho/para ratio of approximately $5\times 10^{-5}$for the conditions in the nucleus of L 1544, and this value is reflected in the H₂D⁺ ortho/para ratio, as shown by GHR. In fact, we find that the H₂D⁺ ortho/para ratio predicted by Eq. (7) of GHR is a good approximation to our numerical results.

The expression of Le Bourlot (1991) for ortho/para H₂ was derived for different conditions (lower density and high abundances of molecular ions such as HCO⁺) and is not a good approximation in the conditions that we have considered. In our situation, the ortho/para H₂ ratio is given approximately by $\zeta /[k_{\rm op}~n(\rm {XH}^{+})]$, where $n(\rm {XH}^+)$ represents the sum of the number densities of ions such as H⁺ and H₃⁺ which convert ortho- to para-H₂ through proton transfer reactions, with a rate coefficient for ortho to para conversion $k_{\rm op}$ which is of order 10^-10 cm^-3 s^-1.

The ortho/para H₂D⁺ ratio has a flat dependence on density, reflecting the behaviour of ortho/para H₂. The ortho forms of both H₂D⁺ and H₃⁺ are produced mainly in reactions of the para forms with ortho-H₂  and so their high ortho/para ratios are attributable to the relatively high ortho-H₂ abundance; Fig. 7 illustrates this point. Because the ortho/para H₂ ratio is not thermalized under the conditions which we consider, the ortho/para H₂D⁺ and H₃⁺ ratios are not thermalized either. It follows that measurements of the ortho/para H₃⁺ ratio should not be used to infer a temperature.

Figure 7: Ortho/para ratios of H3+ and H2D+ as functions of the molecular hydrogen ortho/para ratio. The conditions are those of our reference model: n(H2) = 106 cm-3, T = 10 K, $a_{\rm g}$ = 0.1 $\mu$m, $\zeta = 3\times 10^{-17}$ s-1.

The basic question that we wish to answer is whether the ortho-H₂D⁺ emission observed in L 1544 can be reproduced by a model which assumes complete heavy element depletion. The abundance of ortho-H₂D⁺  relative to H₂  observed towards the peak of L 1544 is $5\times 10^{-10}$ according to CvTCB. In Fig. 8 are shown the computed values of [ortho-H₂D⁺]/[H₂] as a function of density for various assumed grain sizes. There is a tendency for [ortho-H₂D⁺]/[H₂] to decrease with increasing density, owing to the fact that deuteration proceeds further down the chain (towards D₃⁺) at high density. Our results are within a factor of 2 of the observed peak ortho-H₂D⁺ abundance, for grain radii below 0.1 $\mu$m. Given the observational uncertainties, this level of agreement might be considered satisfactory. It is also significant that computations for a = 0.01 $\mu$m show that the computed ortho-H₂D⁺ abundance decreases for all values of n(H₂) relative to the results in Fig. 8 for a = 0.025 $\mu$m, and thus grain radii in the range 0.025-0.05 $\mu$m maximize the ortho-H₂D⁺ abundance. We note that the "observed'' abundance assumes a rate of collisional excitation of H₂D⁺, a temperature, and a density distribution, of which the latter two are controversial (see, for example, Galli et al. 2002). Higher angular resolution studies of the mm-submm continuum are needed for further progress.

One result from our study is that a relatively high rate of cosmic ray ionization is required in the nucleus of L 1544 to account for the observed ortho-H₂D⁺ emission; this implies the penetration of low energy (100 MeV) cosmic rays into the high density core. Cosmic ray fluxes higher than assumed here would enhance the computed ortho-H₂D⁺ abundance; but to a higher flux would correspond an increased rate of heat input to the gas and temperatures in excess of those observed in NH₃ and other tracers.

Figure 8: The computed ortho-H2D+ abundance from our reference model (T = 10 K and $\zeta = 3\times 10^{-17}$ s-1) as a function of density for several values of the assumed grain radius. The values observed by Caselli et al. (2003) are shown for comparison.

3.5 Atomic deuterium in the completely depleted core

An interesting by-product of our work is the prediction of the atomic D abundance under the conditions of complete depletion. Figure 9 shows the abundances of atomic deuterium and hydrogen as functions of density for various grain radii. Atomic H is produced mainly by cosmic ray ionization and destroyed by the formation of H₂ on grains. As a result, the fractional abundance of H varies inversely with the density and decreases with the grain radius, owing to the greater grain surface area available to form H₂. Atomic D, on the other hand, is produced by the recombination of deuterated ions and becomes relatively more abundant as deuterium fractionation increases. It is noteworthy that the abundance of atomic D found from our computations [>0.1 cm^-3 for n(H₂) = 10⁶ cm^-3] is of the order of that required to explain the high fractionation of deuterium in methanol and formaldehyde on grain surfaces (Caselli et al. 2002).

Figure 9: Atomic H (full lines) and atomic D (broken lines) densities for four values of the grain radius (increasing from bottom to top) for densities between 105 and 107 cm-3.

4 Discussion

4.1 Grain size distribution in depleted cores

The grain size distribution in prestellar cores is not known. It depends on the timescales for grain coagulation and mantle accretion, relative to the dynamical and chemical timescales. In view of this uncertainty, we compared the results of computations assuming a unique grain radius with the corresponding calculations assuming a MRN power-law size distribution; we retained the same value of the grain surface area per H-nucleus, $n_{\rm g}~\sigma_{\rm g}/n_{\rm H}$, in both cases. In calculations with a unique value of $a_{\rm g}$, $n_{\rm g}~\sigma_{\rm g}/n_{\rm H} = 1.1\times 10^{-21}~ (0.10/a_{\rm g}$($\mu$m)) cm². When a power law dependence $a_{\rm g}^{-r}$ is adopted, providing the exponent r > 3 (MRN derived r = 3.5), the lower limit to the radius essentially determines the grain surface area. When r = 3.5, the lower limit to $a_{\rm g}$ is 600 Å for $n_{\rm g}~ \sigma _{\rm g}/n_{\rm H} = 1.1\times 10^{-21}$ cm². We found that the results of calculations with the same value of grain surface area per H-nucleus were quantitatively similar for the cases of a unique grain radius and a MRN power-law size distribution.

As shown above, the detection of H₂D⁺ in L 1544 by CvTCB can be explained if $a_{\rm g} < 0.1$ $\mu$m  or, equivalently, if the grain surface area per H-nucleus is at least $1.1\times 10^{-21}$ cm², which is close to the value for grains in diffuse clouds. This conclusion is consistent with deductions which have been made from the mm-infrared extinction curve in prestellar regions: Bianchi et al. (2003) found that the ratio of mm to NIR dust absorption is marginally larger in prestellar environments than in diffuse clouds; it is also consistent with the simulations of Ossenkopf & Henning (1994).

Under conditions of complete heavy element depletion, the thickness of the grain mantle is expected to be roughly 40 percent of the grain radius, or 400 Å for 0.1 $\mu$m grains. Then, 45 percent of the grain mass is in the refractory core and 55 percent in the ice mantle. A MRN size distribution yields an ice mantle 260 Å thick, independent of the grain radius: see Appendix B.

4.2 Ortho to para ratios

Our calculations show that the ortho/para ratios of both H₃⁺ and H₂D⁺ are directly dependent on the H₂ ortho/para ratio. Thus, measurements of the former would yield information on the latter. We emphasize that, whilst steady state may be a reasonable approximation for most of the chemical species that we have discussed, this may not be the case for the H₂ ortho/para ratio. Thus, reliable measurements of the ortho/para ratios of either H₃⁺ or H₂D⁺could provide information on the evolution of the source.

Unfortunately, the 1₀₁-0₀₀ transition to the ground state of para-H₂D⁺ (at 1.37 THz) will be difficult to excite at the temperatures of prestellar cores and may be observable only in absorption towards a background continuum source. On the other hand, H₃⁺ has no allowed submm transitions and is best observed in the NIR at 3.5 $\mu$m  in absorption towards background sources (McCall et al. 1999). It has been assumed that the H₃⁺ which has been observed exists in regions whose composition is similar to that of "normal'' interstellar molecular clouds and, furthermore, that the ortho/para forms are in LTE. In our view, both of these assumptions may prove to be invalid. Observations of H₂D⁺ and the other members of the family may help to shed light on this issue.

4.3 Observations of completely depleted cores

Observations of pre-protostellar cores, in which heavy elements have frozen on to grains, are important because they provide us with information on the structure and kinematics of the core at the time of its initial gravitational collapse. It is important to observe the mm continuum emission of these objects with high angular resolution in order to derive the density and temperature distributions for comparison with theoretical determinations (e.g. Galli et al. 2002). Likewise, kinematical information may be obtained from higher angular resolution studies of the 1₁₀-1₁₁line of H₂D⁺. The present study shows that it is plausible to conclude that the ortho-H₂D⁺ emission detected by CvTCB towards L 1544 arises in completely depleted material. Indeed, these authors failed by a factor of 10 to account for their observations of L 1544 using estimates based on less extreme assumptions regarding the depletion of heavy elements from the gas phase.

Our results can be compared with those of Roberts et al. (2003) who carried out a time dependent calculation of the chemistry in a prestellar core with density $3\times 10^6$ cm^-3  including freeze-out. They found, as do we, that D₃⁺ can be the most abundant ion, though their results suggest that H₃⁺ becomes the major ion when all heavy species freeze out; this may be attributable to a different way of dealing with grain neutralization of ions.

Further progress in this area will depend on our capacity to sample transitions other than the the 1₁₀-1₁₁ line of H₂D⁺. H₃⁺ and D₃⁺ do not have permanent dipole moments, and hence their detection is feasible only in absorption in the near infrared. Our calculations suggest that searches for D₃⁺ should be directed towards regions where high heavy element depletion is believed to have occurred. Searches for emission in low-lying transitions of D₂H⁺ should also be made: the abundance of D₂H⁺ in conditions of complete depletion is predicted to be comparable to that of H₂D⁺. Under conditions of LTE at temperature T, the relative populations of the lowest ortho (0₀₀) and para (1₀₁) levels of D₂H⁺ would be

$$\frac{n(1_{01})}{n(0_{00})} = (9/6){\rm exp}(-50.2/T)$$

where the nuclear spin statistical weights (2I + 1) of the ortho (I = 0, 2) and the para (I = 1) levels are 6 and 3 and their rotational degeneracies (2J + 1) are 1 and 3, respectively, yielding total statistical weights of 6 and 9 for 0₀₀ and 1₀₁, respectively. However, under the conditions in prestellar cores, the para/ortho D₂H⁺ ratio (like the ortho/para ratios discussed above, and for related reasons) is probably much in excess of its LTE value; we estimate that the para and ortho forms of this species may have comparable abundances. Thus a deep search for the 692 GHz 1₁₀-1₀₁transition of D₂H⁺ seems worthwhile.

5 Conclusions

We summarize here our most important conclusions.

• We believe that our data strongly suggest that the detection of H₂D⁺ by CvTCB is due to complete freeze-out of the heavy elements in the core of L 1544. If so, only species such as H₃⁺ and its deuterated counterparts are useful tracers of the kinematics in the dense centre of L 1544 and other similar objects.
• Knowledge of the mean grain surface area is needed in order to understand phenomena such as deuterium fractionation; it seems likely that this remains true even if some of the heavy species are in the gas phase. Our results imply that, in L 1544, the mean grain surface area per H-atom is above 10^-21 cm², which is comparable to values found in the diffuse interstellar medium.
• The degree of ionization and the timescale for ambipolar diffusion, determined under conditions of complete heavy element depletion, are not significantly different from values reported elsewhere in the literature, evaluated assuming degrees of depletion more typical of dense interstellar clouds. However, we estimate that, under the conditions considered here, the effect of the coupling between the neutral gas and charged grains is to increase the ambipolar diffusion timescale by between one and two orders of magnitude; this finding will be the subject of further investigation.

Acknowledgements

It is a pleasure to thank Evelyne Roueff and Eric Herbst for informative and helpful discussions, Paola Caselli for her comments on the text, and the referee, Ted Bergin, for a constructive and thought-provoking report. CMW thanks the Ecole Normale Supérieure in Paris for hospitality during a stay when part of the work described here was carried out as well as the Max Planck Institute in Bonn for hospitality on diverse occasions.

References

 
• Aikawa, Y., Ohashi, N., & Herbst, E. 2003, ApJ, 593, 906 NASA ADS  
• André, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Univ. of Arizona)  
• Belloche, A. 2002, Ph.D. Thesis, Univ. Paris  
• Bergin, E. A., & Langer, W. D. 1997, ApJ, 486, 316 NASA ADS  
• Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101 NASA ADS  
• Bianchi, S., Goncalves, J., Albrecht, M., et al. 2003, A&A, 399, L43 NASA ADS  
• Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002a, ApJ, 565, 331 NASA ADS  
• Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002b, ApJ, 523, 344  
• Caselli, P., Stantcheva, T., Shalabiea, O., Shematovich, V. I., & Herbst, E. 2002c, Plan. Sp. Sci., 50, 1257 NASA ADS  
• Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37 (CvTCB) NASA ADS  
• Ciolek, G. E., & Basu, S. 2000, ApJ, 529, 925 NASA ADS  
• Dalgarno, A., & Lepp, S. 1984, ApJ, 287, L47 NASA ADS  
• Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 NASA ADS  
• Flower, D. R. 2000, MNRAS, 313, L19 NASA ADS  
• Flower, D. R., Le Bourlot, J., Pineau des Forêts, G., & Cabrit, S. 2003, MNRAS, 341, 1272 NASA ADS  
• Flower, D. R., & Pineau des Forêts, G. 2003, MNRAS, 343, 390 (FPdesF) NASA ADS  
• Galli, D., Walmsley, C. M., & Gonçalves, J. 2002, A&A, 394, 275 NASA ADS  
• Gerlich, D., Herbst, E., & Roueff, E. 2002, Plan. Sp. Sci., 50, 1275 (GHR) NASA ADS  
• Giles, K., Adams, N. G., & Smith, D. 1992, J. Phys. Chem., 96, 7645  
• Le Bourlot, J. 1991, A&A, 242, 235 NASA ADS  
• Le Bourlot, J., Pineau des Forêts, G., & Roueff, E. 1995, A&A, 297, 251 NASA ADS  
• Linsky, J. 2003, Space Sci. Rev., 106, 49 NASA ADS  
• Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 NASA ADS  
• McCall, B. J., Geballe, T. R., Hinkle, K. H., & Oka, T. 1999, ApJ, 522, 338 NASA ADS  
• McCall, B. J., Huneycutt, A. J., Saykally, R. J., et al. 2003, Nature, 422, 500 NASA ADS  
• McKee, C. F. 1989, ApJ, 345, 782 NASA ADS  
• Nakano, T., & Umebayashi, T. 1980, Publ. Astron. Soc. Japan, 32, 613 NASA ADS  
• Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943 NASA ADS  
• Pagani, L., Wannier, P. G., Frerking, M. A., et al. 1992, A&A, 258, 472 NASA ADS  
• Ramanlal, J., Polyansky, O. L., & Tennyson, J. 2003, A&A, 406, 383 NASA ADS  
• Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41 NASA ADS  
• Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 NASA ADS  
• Tafalla, M., Myers, P. C., Caselli, P., Walmsley, C. M., & Comito, C. 2002, ApJ, 569, 815 NASA ADS  
• Uy, D., Cordonnier, M., & Oka, T. 1997, Phys. Rev. Lett., 78, 3844 NASA ADS  
• Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, Proc. Fourth Cologne-Bonn-Zermatt Symp., The Dense Interstellar Medium in Galaxies (Springer Verlag)  
• Watson, W. D. 1978, ARA&A, 16, 585 NASA ADS

  
6 Online Material

Appendix A: Chemical reactions and rate coefficients

In this Appendix are specified the reactions included and the rate coefficients adopted in our chemical model.

Protons have nuclear spin I = 1/2 and are fermions. The ortho and para forms of H₂, H₂⁺, H₃⁺ and H₂D⁺ have been treated as distinct species. Ortho- and para-H₂ were assumed to form hot, on grains, in the ratio 3:1 of their nuclear spin statistical weights. Subsequent proton-exchange reactions interconvert the ortho and para forms. Such reactions were assumed to involve only the rotational ground states of the ortho and para modifications, and the corresponding rate coefficients take account of the relative statistical weights (2I+1)(2J+1) of the products and reactants.

Deuteration of H₃⁺ in reactions with HD may be viewed as proton transfer (from H₃⁺ to HD) or deuteron-proton exchange; we allowed for both viewpoints in our reaction set. Rate coefficients for dissociative recombination, with electrons in the gas phase or on negatively charged grains, were assumed to be the same for the ortho and para forms of the recombining ion.

We note that deuterons have nuclear spin I = 1 and are bosons. Multiply deuterated species have ortho and para forms. In the present study, we have not attempted to allow for the consequences of Bose-Einstein statistics when evaluating the rates of reactions involving multiply-deuterated species.

 

 

Table A.1: Reactions and rate coefficients adopted in our chemical model. The parameters $\alpha$, $\beta$, and $\gamma$ define the rate coefficients k cm³ s^-1 at temperature T through the relation k = J $(a_{\rm g}, T)\gamma (T/300)^{\alpha } {\rm exp}(-\beta /T)$; J allows for Coulomb focusing in reactions of positive ions and negatively charged grains (Draine & Sutin 1987, Eq. (3.4)); $a_{\rm g}$ = 0.1 $\mu$m is adopted in the table. The rates s^-1 of reactions induced directly by cosmic rays (crp) or indirectly by the secondary electrons, which generate H₂fluorescence photons (secpho), are given by $\gamma \zeta$, where $\zeta$is the rate of cosmic ray ionization of H₂. In the first three reactions (formation of H₂ or HD on grain surfaces), $\gamma$ denotes the fraction of such reactions which form the specified product; the rate is calculated internally in the program from the grain parameters. Key reactions are in bold face. Numbers in parentheses are powers of 10.

Reaction	$\gamma$	$\alpha$	$\beta$
H + H $\rightarrow$ H2(p)	0.25
H + H $\rightarrow$ H2(o)	0.75
H + D $\rightarrow$ HD	1.00
H + crp $\rightarrow$ H+ + e-	0.46
He + crp $\rightarrow$ He+ + e-	0.50
H2(p) + crp $\rightarrow$ H+ + H + e-	0.04
H2(o) + crp $\rightarrow$ H+ + H + e-	0.04
H2(p) + crp $\rightarrow$ H + H	1.50
H2(o) + crp $\rightarrow$ H + H	1.50
H2(p) + crp $\rightarrow$ H2+(p) + e-	0.96
H2(o) + crp $\rightarrow$ H2+(o) + e-	0.96
H+ + H2(o) $\rightarrow$ H+ + H2(p)	2.20(-10)	0.00	0.0
H3+(p) + H2(o) $\rightarrow$ H3+(o) + H2(p)	4.40(-10)	0.00	0.0
H3+(o) + H2(o) $\rightarrow$ H3+(p) + H2(p)	1.10(-10)	0.00	0.0
H+ + H2(p) $\rightarrow$ H+ + H2(o)	1.98(-09)	0.00	170.5
H3+(p) + H2(p) $\rightarrow$ H3+(o) + H2(o)	1.98(-09)	0.00	203.4
H3+(o) + H2(p) $\rightarrow$ H3+(p) + H2(o)	1.98(-09)	0.00	137.6
H2+(o) + H2(o) $\rightarrow$ H3+(p) + H	1.05(-09)	0.00	0.0
H2+(o) + H2(o) $\rightarrow$ H3+(o) + H	1.05(-09)	0.00	0.0
H2+(p) + H2(o) $\rightarrow$ H3+(p) + H	1.05(-09)	0.00	0.0
H2+(p) + H2(o) $\rightarrow$ H3+(o) + H	1.05(-09)	0.00	0.0
H2+(o) + H2(p) $\rightarrow$ H3+(p) + H	1.05(-09)	0.00	0.0
H2+(o) + H2(p) $\rightarrow$ H3+(o) + H	1.05(-09)	0.00	0.0
H2+(p) + H2(p) $\rightarrow$ H3+(p) + H	2.10(-09)	0.00	0.0
He+ + H2(p) $\rightarrow$ H+ + H + He	1.10(-13)	-0.24	0.0
He+ + H2(o) $\rightarrow$ H+ + H + He	1.10(-13)	-0.24	0.0
H+ + e- $\rightarrow$ H + photon	3.61(-12)	-0.75	0.0
H2+(p) + e- $\rightarrow$ H + H	1.60(-08)	-0.43	0.0
H2+(o) + e- $\rightarrow$ H + H	1.60(-08)	-0.43	0.0
He+ + e- $\rightarrow$ He + photon	4.50(-12)	-0.67	0.0
H3+(p) + e- $\rightarrow$ H + H + H	5.10(-08)	-0.52	0.0
H3+(o) + e- $\rightarrow$ H + H + H	5.10(-08)	-0.52	0.0
H3+(p) + e- $\rightarrow$ H2(p) + H	1.13(-08)	-0.52	0.0
H3+(p) + e- $\rightarrow$ H2(o) + H	0.57(-08)	-0.52	0.0
H3+(o) + e- $\rightarrow$ H2(o) + H	1.70(-08)	-0.52	0.0
D + crp $\rightarrow$ D+ + e-	0.46
HD + crp $\rightarrow$ H+ + D + e-	0.02
HD + crp $\rightarrow$ D+ + H + e-	0.02
HD + crp $\rightarrow$ H + D	1.50
HD + crp $\rightarrow$ HD+ + e-	0.96
D+ + H $\rightarrow$ H+ + D	1.00(-09)	0.00	0.0
H+ + D $\rightarrow$ D+ + H	1.00(-09)	0.00	41.0
D+ + H2(p) $\rightarrow$ H+ + HD	2.10(-09)	0.00	0.0
D+ + H2(o) $\rightarrow$ H+ + HD	2.10(-09)	0.00	0.0
H+ + HD $\rightarrow$ D+ + H2(p)	1.00(-09)	0.00	464.0
H+ + HD $\rightarrow$ D+ + H2(o)	1.00(-09)	0.00	634.5
HD+ + H2(p) $\rightarrow$ H2D+(p) + H	5.25(-10)	0.00	0.0
HD+ + H2(p) $\rightarrow$ H2D+(o) + H	5.25(-10)	0.00	0.0
HD+ + H2(o) $\rightarrow$ H2D+(p) + H	5.25(-10)	0.00	0.0
HD+ + H2(o) $\rightarrow$ H2D+(o) + H	5.25(-10)	0.00	0.0
HD+ + H2(p) $\rightarrow$ H3+(p) + D	1.05(-09)	0.00	0.0
HD+ + H2(o) $\rightarrow$ H3+(o) + D	5.25(-09)	0.00	0.0
HD+ + H2(o) $\rightarrow$ H3+(p) + D	5.25(-09)	0.00	0.0
H2+(p) + HD $\rightarrow$ H2D+(p) + H	5.25(-10)	0.00	0.0
H2+(p) + HD $\rightarrow$ H2D+(o) + H	5.25(-10)	0.00	0.0
H2+(o) + HD $\rightarrow$ H2D+(p) + H	5.25(-10)	0.00	0.0
H2+(o) + HD $\rightarrow$ H2D+(o) + H	5.25(-10)	0.00	0.0
H2+(p) + HD $\rightarrow$ H3+(p) + D	1.05(-09)	0.00	0.0
H2+(o) + HD $\rightarrow$ H3+(p) + D	5.25(-10)	0.00	0.0
H2+(o) + HD $\rightarrow$ H3+(o) + D	5.25(-10)	0.00	0.0
H3+(p) + D $\rightarrow$ H2D+(p) + H	0.66(-09)	0.00	0.0
H3+(p) + D $\rightarrow$ H2D+(o) + H	0.33(-09)	0.00	0.0
H3+(o) + D $\rightarrow$ H2D+(o) + H	1.00(-09)	0.00	0.0
H2D+(p) + H $\rightarrow$ H3+(p) + D	1.00(-09)	0.00	632.0
H2D+(o) + H $\rightarrow$ H3+(o) + D	0.50(-09)	0.00	578.4
H2D+(o) + H $\rightarrow$ H3+(p) + D	0.50(-09)	0.00	545.5
H3+(p) + HD $\rightarrow$ H2D+(p) + H2(p)	1.17(-10)	0.00	0.0
H3+(p) + HD $\rightarrow$ H2D+(p) + H2(o)	1.17(-10)	0.00	0.0
H3+(p) + HD $\rightarrow$ H2D+(o) + H2(p)	5.83(-11)	0.00	0.0
H3+(p) + HD $\rightarrow$ H2D+(o) + H2(o)	5.83(-11)	0.00	25.0
H3+(o) + HD $\rightarrow$ H2D+(p) + H2(o)	1.75(-10)	0.00	0.0
H3+(o) + HD $\rightarrow$ H2D+(o) + H2(o)	1.75(-10)	0.00	0.0
H2D+(p) + H2(o) $\rightarrow$ H2D+(o) + H2(p)	1.98(-09)	0.00	0.0
H2D+(o) + H2(p) $\rightarrow$ H2D+(p) + H2(o)	1.98(-09)	0.00	84.0
H2D+(o) + H2(o) $\rightarrow$ H2D+(p) + H2(p)	2.44(-11)	0.00	0.0
H2D+(p) + H2(p) $\rightarrow$ H2D+(o) + H2(o)	1.98(-09)	0.00	257.0
H2D+(p) + H2(p) $\rightarrow$ H3+(p) + HD	1.40(-10)	0.00	232.0
H2D+(p) + H2(o) $\rightarrow$ H3+(p) + HD	0.70(-10)	0.00	61.5
H2D+(p) + H2(o) $\rightarrow$ H3+(o) + HD	0.70(-10)	0.00	94.4
H2D+(o) + H2(p) $\rightarrow$ H3+(p) + HD	1.40(-10)	0.00	145.5
H2D+(o) + H2(o) $\rightarrow$ H3+(p) + HD	0.70(-10)	0.00	0.0
H2D+(o) + H2(o) $\rightarrow$ H3+(o) + HD	0.70(-10)	0.00	8.0
He+ + HD $\rightarrow$ H+ + D + He	5.50(-14)	-0.24	0.0
He+ + HD $\rightarrow$ D+ + H + He	5.50(-14)	-0.24	0.0
D+ + e- $\rightarrow$ D + photon	3.61(-12)	-0.75	0.0
HD+ + e- $\rightarrow$ H + D	3.40(-09)	-0.40	0.0
H2D+(p) + e- $\rightarrow$ H + H + D	4.96(-08)	-0.52	0.0
H2D+(o) + e- $\rightarrow$ H + H + D	4.96(-08)	-0.52	0.0
H2D+(p) + e- $\rightarrow$ HD + H	1.36(-08)	-0.52	0.0
H2D+(o) + e- $\rightarrow$ HD + H	1.36(-08)	-0.52	0.0
H2D+(p) + e- $\rightarrow$ H2(p) + D	4.76(-09)	-0.52	0.0
H2D+(o) + e- $\rightarrow$ H2(o) + D	4.76(-09)	-0.52	0.0
D+ + H $\rightarrow$ HD+ + photon	3.90(-19)	1.80	0.0
H+ + D $\rightarrow$ HD+ + photon	3.90(-19)	1.80	0.0
H2+(p) + D $\rightarrow$ H2D+(p) + photon	7.00(-18)	1.80	0.0
HD+ + H $\rightarrow$ H2D+(p) + photon	1.20(-17)	1.80	0.0
HD+ + H $\rightarrow$ H+ + HD	6.40(-10)	0.00	0.0
H2+(p) + D $\rightarrow$ D+ + H2(p)	6.40(-10)	0.00	0.0
HD+ + H $\rightarrow$ H2+(p) + D	1.00(-09)	0.00	154.0
H2+(p) + D $\rightarrow$ HD+ + H	1.00(-09)	0.00	0.0
D2 + crp $\rightarrow$ D+ + D + e-	0.04
D2 + crp $\rightarrow$ D + D	1.50
D2 + crp $\rightarrow$ D2+ + e-	0.96
H+ + D2 $\rightarrow$ D+ + HD	2.10(-09)	0.00	491.0
D+ + HD $\rightarrow$ H+ + D2	1.00(-09)	0.00	0.0
D2+ + H2(p) $\rightarrow$ H2D+(p) + D	1.05(-09)	0.00	0.0
D2+ + H2(o) $\rightarrow$ H2D+(o) + D	1.05(-09)	0.00	0.0
D2+ + H2(p) $\rightarrow$ HD2+ + H	1.05(-09)	0.00	0.0
D2+ + H2(o) $\rightarrow$ HD2+ + H	1.05(-09)	0.00	0.0
H2+(p) + D2 $\rightarrow$ H2D+(p) + D	1.05(-09)	0.00	0.0
H2+(o) + D2 $\rightarrow$ H2D+(o) + D	1.05(-09)	0.00	0.0
H2+(p) + D2 $\rightarrow$ HD2+ + H	1.05(-09)	0.00	0.0
H2+(o) + D2 $\rightarrow$ HD2+ + H	1.05(-09)	0.00	0.0
HD+ + HD $\rightarrow$ H2D+(p) + D	5.25(-10)	0.00	0.0
HD+ + HD $\rightarrow$ H2D+(o) + D	5.25(-10)	0.00	0.0
HD+ + HD $\rightarrow$ HD2+ + H	1.05(-09)	0.00	0.0
H2D+(p) + D $\rightarrow$ HD2+ + H	1.00(-09)	0.00	0.0
H2D+(o) + D $\rightarrow$ HD2+ + H	1.00(-09)	0.00	0.0
HD2+ + H $\rightarrow$ H2D+(p) + D	1.00(-09)	0.00	600.0
HD2+ + H $\rightarrow$ H2D+(o) + D	1.00(-09)	0.00	513.5
H3+(o) + D2 $\rightarrow$ H2D+(o) + HD	3.50(-10)	0.00	0.0
H3+(o) + D2 $\rightarrow$ HD2+ + H2(o)	3.50(-10)	0.00	0.0
H3+(p) + D2 $\rightarrow$ HD2+ + H2(p)	2.33(-10)	0.00	0.0
H3+(p) + D2 $\rightarrow$ HD2+ + H2(o)	1.17(-10)	0.00	0.0
H3+(p) + D2 $\rightarrow$ H2D+(p) + HD	2.33(-10)	0.00	0.0
H3+(p) + D2 $\rightarrow$ H2D+(o) + HD	1.17(-10)	0.00	0.0
H2D+(p) + HD $\rightarrow$ HD2+ + H2(p)	2.60(-10)	0.00	0.0
H2D+(o) + HD $\rightarrow$ HD2+ + H2(o)	2.60(-10)	0.00	0.0
H2D+(o) + HD $\rightarrow$ H3+(o) + D2	1.30(-10)	0.00	108.0
H2D+(o) + HD $\rightarrow$ H3+(p) + D2	1.30(-10)	0.00	75.0
H2D+(p) + HD $\rightarrow$ H3+(p) + D2	2.60(-10)	0.00	161.5
HD2+ + H2(p) $\rightarrow$ H2D+(p) + HD	1.00(-10)	0.00	190.0
HD2+ + H2(p) $\rightarrow$ H2D+(o) + HD	1.00(-10)	0.00	277.0
HD2+ + H2(o) $\rightarrow$ H2D+(o) + HD	1.00(-10)	0.00	106.0
HD2+ + H2(o) $\rightarrow$ H2D+(p) + HD	1.00(-10)	0.00	20.1
HD2+ + H2(p) $\rightarrow$ H3+(p) + D2	2.00(-10)	0.00	320.0
HD2+ + H2(o) $\rightarrow$ H3+(p) + D2	1.00(-10)	0.00	149.5
HD2+ + H2(o) $\rightarrow$ H3+(o) + D2	1.00(-10)	0.00	182.4
He+ + D2 $\rightarrow$ D+ + D + He	1.10(-13)	-0.24	0.0
He+ + D2 $\rightarrow$ D2+ + He	2.50(-14)	0.00	0.0
D2+ + e- $\rightarrow$ D + D	3.40(-09)	-0.40	0.0
HD2+ + e- $\rightarrow$ D + D + H	4.96(-08)	-0.52	0.0
HD2+ + e- $\rightarrow$ HD + D	1.36(-08)	-0.52	0.0
HD2+ + e- $\rightarrow$ D2 + H	4.76(-09)	-0.52	0.0
D2+ + H $\rightarrow$ H+ + D2	6.40(-10)	0.00	0.0
HD+ + D $\rightarrow$ D+ + HD	6.40(-10)	0.00	0.0
HD+ + D $\rightarrow$ D2+ + H	1.00(-09)	0.00	0.0
D2+ + H $\rightarrow$ HD+ + D	1.00(-09)	0.00	472.0
D2+ + D $\rightarrow$ D+ + D2	6.40(-10)	0.00	0.0
D2+ + HD $\rightarrow$ HD2+ + D	1.05(-09)	0.00	0.0
D2+ + HD $\rightarrow$ D3+ + H	1.05(-09)	0.00	0.0
HD+ + D2 $\rightarrow$ HD2+ + D	1.05(-09)	0.00	0.0
HD+ + D2 $\rightarrow$ D3+ + H	1.05(-09)	0.00	0.0
HD2+ + D $\rightarrow$ D3+ + H	1.00(-09)	0.00	0.0
D3+ + H $\rightarrow$ HD2+ + D	1.00(-09)	0.00	655.0
HD2+ + HD $\rightarrow$ D3+ + H2(p)	1.00(-10)	0.00	0.0
HD2+ + HD $\rightarrow$ D3+ + H2(o)	1.00(-10)	0.00	0.0
HD2+ + HD $\rightarrow$ H2D+(p) + D2	1.00(-10)	0.00	110.0
HD2+ + HD $\rightarrow$ H2D+(o) + D2	1.00(-10)	0.00	196.5
H2D+(p) + D2 $\rightarrow$ HD2+ + HD	8.50(-10)	0.00	0.0
H2D+(o) + D2 $\rightarrow$ HD2+ + HD	8.50(-10)	0.00	0.0
H2D+(p) + D2 $\rightarrow$ D3+ + H2(p)	8.50(-10)	0.00	0.0
H2D+(o) + D2 $\rightarrow$ D3+ + H2(o)	8.50(-10)	0.00	0.0
D3+ + H2(p) $\rightarrow$ H2D+(p) + D2	1.50(-09)	0.00	352.0
D3+ + H2(o) $\rightarrow$ H2D+(o) + D2	1.50(-09)	0.00	268.0
D3+ + H2(p) $\rightarrow$ HD2+ + HD	1.50(-09)	0.00	241.0
D3+ + H2(o) $\rightarrow$ HD2+ + HD	1.50(-09)	0.00	70.5
D3+ + e- $\rightarrow$ D + D + D	2.03(-08)	-0.50	0.0
D3+ + e- $\rightarrow$ D2 + D	6.75(-09)	-0.50	0.0
D2+ + D2 $\rightarrow$ D3+ + D	2.10(-09)	0.00	0.0
D3+ + HD $\rightarrow$ HD2+ + D2	1.50(-09)	0.00	161.5
g0 + secpho $\rightarrow$ g+ + e-	0.63(08)
g- + secpho $\rightarrow$ g0 + e-	0.41(09)
g0 + e- $\rightarrow$ g- + photon	6.90(-05)	0.50	0.0
g- + H+ $\rightarrow$ g0 + H	1.60(-06)	0.50	0.0
g- + H3+(p) $\rightarrow$ g0 + H2(p) + H	3.07(-07)	0.50	0.0
g- + H3+(p) $\rightarrow$ g0 + H2(o) + H	1.54(-07)	0.50	0.0
g- + H3+(o) $\rightarrow$ g0 + H2(o) + H	4.61(-07)	0.50	0.0
g- + H3+(p) $\rightarrow$ g0 + 3H	4.61(-07)	0.50	0.0
g- + H3+(o) $\rightarrow$ g0 + 3H	4.61(-07)	0.50	0.0
g- + He+ $\rightarrow$ g0 + He	8.00(-07)	0.50	0.0
g0 + H+ $\rightarrow$ g+ + H	1.60(-06)	0.50	0.0
g0 + H3+(p) $\rightarrow$ g+ + H2(p) + H	3.07(-07)	0.50	0.0
g0 + H3+(p) $\rightarrow$ g+ + H2(o) H	1.54(-07)	0.50	0.0
g0 + H3+(o) $\rightarrow$ g+ + H2(o) H	4.61(-07)	0.50	0.0
g0 + H3+(p) $\rightarrow$ g+ + 3H	4.61(-07)	0.50	0.0
g0 + H3+(o) $\rightarrow$ g+ + 3H	4.61(-07)	0.50	0.0
g0 + He+ $\rightarrow$ g+ + He	8.00(-07)	0.50	0.0
g+ + e- $\rightarrow$ g0 + photon	6.90(-05)	0.50	0.0
g- + D+ $\rightarrow$ g0 + D	1.13(-06)	0.50	0.0
g- + H2D+(p) $\rightarrow$ g0 + H2(p) + D	1.33(-07)	0.50	0.0
g- + H2D+(p) $\rightarrow$ g0 + HD + H	2.66(-07)	0.50	0.0
g- + H2D+(o) $\rightarrow$ g0 + H2(o) + D	1.33(-07)	0.50	0.0
g- + H2D+(o) $\rightarrow$ g0 + HD + H	2.66(-07)	0.50	0.0
g- + H2D+(p) $\rightarrow$ g0 + D + 2H	3.99(-07)	0.50	0.0
g- + H2D+(o) $\rightarrow$ g0 + D + 2H	3.99(-07)	0.50	0.0
g0 + D+ $\rightarrow$ g+ + D	1.13(-06)	0.50	0.0
g0 + H2D+(p) $\rightarrow$ g+ + H2(p) + D	1.33(-07)	0.50	0.0
g0 + H2D+(p) $\rightarrow$ g+ + HD + H	2.66(-07)	0.50	0.0
g0 + H2D+(o) $\rightarrow$ g+ + H2(o) + D	1.33(-07)	0.50	0.0
g0 + H2D+(o) $\rightarrow$ g+ + HD + H	2.66(-07)	0.50	0.0
g0 + H2D+(p) $\rightarrow$ g+ + D + 2H	3.99(-07)	0.50	0.0
g0 + H2D+(o) $\rightarrow$ g+ + D + 2H	3.99(-07)	0.50	0.0
g- + HD2+ $\rightarrow$ g0 + HD + D	2.38(-07)	0.50	0.0
g- + HD2+ $\rightarrow$ g0 + D2 + H	1.19(-07)	0.50	0.0
g- + HD2+ $\rightarrow$ g0 + 2D + H	3.57(-07)	0.50	0.0
g0 + HD2+ $\rightarrow$ g+ + HD + D	2.38(-07)	0.50	0.0
g0 + HD2+ $\rightarrow$ g+ + D2 + H	1.19(-07)	0.50	0.0
g0 + HD2+ $\rightarrow$ g+ + 2D + H	3.57(-07)	0.50	0.0
g- + D3+ $\rightarrow$ g0 + D2 + D	3.26(-07)	0.50	0.0
g- + D3+ $\rightarrow$ g0 + 3D	3.26(-07)	0.50	0.0
g0 + D3+ $\rightarrow$ g+ + D2 + D	3.26(-07)	0.50	0.0
g0 + D3+ $\rightarrow$ g+ + 3D	3.26(-07)	0.50	0.0

Appendix B: Ice Mantle Thickness

In this Appendix, we derive estimates of the ice-mantle thickness in completely depleted cores. We follow the approach of FPdesF.

Following FPdesF, we assume the refractory cores to be composed of silicates (in practice, olivine, MgFeSiO₄) or graphite. The mass of the cores, per unit volume of gas, is  $\sum n({\rm X}) M$, where the summation extends over the elements X, and M denotes the elemental mass. We define $\alpha_{\rm c}$ through the relation $\sum n({\rm X}) M = \alpha_{\rm c}~ n_{\rm H} m_{\rm H}$, and $\alpha_{\rm c} = 7.84\times 10^{-3}$, which corresponds to a fraction of 0.0056 by mass relative to the gas-phase H and He.

Assuming that all elements heavier than He are in the grain cores or frozen in the mantles, the total dust-to-gas mass density ratio is 0.0125, and the mass per unit volume of the gas of the ice in the mantle, is $\alpha_{\rm i}~ n_{\rm H}m_{\rm H}$, where $\alpha_{\rm i} = 9.63 \times 10^{-3}$.

The volume $V_{\rm i}$of a grain ice mantle is determined by the relation

$$n_{\rm c}~ V_{\rm i}~ \rho_{\rm i} = \alpha_{\rm i}~ n_{\rm H} m_{\rm H}$$ (B.1)

where $n_{\rm c}$ is the number of grain cores (i.e. the number of grains) per unit volume of the gas and $\rho_{\rm i}$ is the density of the mantle material; $n_{\rm c}$ is defined by

$$n_{\rm c}~ V_{\rm c}~ \rho_{\rm c} = \alpha_{\rm c}~ n_{\rm H} m_{\rm H}$$ (B.2)

where $\rho_{\rm c}$ is the density and $V_{\rm c}$ is the volume of refractory core material. We adopt $\rho_{\rm c}$ = 3 gm cm^-3 and $\rho_{\rm i}$ = 1 gm cm^-3. Hence, the ratio of mantle to core volume is

$$V_{\rm i}/V_{\rm c} = \frac{\alpha_{\rm i}\rho_{\rm c}} {\alpha_{\rm c}\rho_{\rm i}}\cdot$$ (B.3)

Substituting the above numbers, we find that the ratio $V_{\rm i}/V_{\rm c}$ = 3.72, or a fraction 0.79 of the grain volume is occupied by the mantle. The ratio of mantle thickness to core radius is 0.67, and the ratio of mantle thickness to grain (core + mantle) radius is 0.40, in the limit of complete heavy element depletion. Thus, for a grain of radius $a_{\rm g}=0.1$ $\mu$m, the mantle thickness is 400 Å.

This analysis can be generalized to the case of a MRN power-law size distribution, using the formalism of Le Bourlot et al. (1995). If the lower limit of the size distribution is taken to be a_m = 0.01 $\mu$m (and the upper limit is a_M = 0.3 $\mu$m), then the thickness $\delta$of the mantle is obtained by solving the equation

$$f_4(V_{\rm i}/V_{\rm c}) = \left(f_1\delta ^3 + 3f_2\delta ^2 + 3f_3\delta \right)$$ (B.4)

where f_n is defined by

$$f_n = \left(a_{M}^{n-r} - a_{m}^{n-r}\right)/(n-r)$$ (B.5)

and r = 3.5 for the MRN grain size distribution. Solving for $\delta$, we find that the mantle thickness is 260 Å, independent of a for the specified values of a_m (and a_M). Thus, very small grains are not present under conditions of complete heavy element depletion, owing to the formation of an ice mantle around the core.