---

Chemical modeling of L183 (L134N): an estimate of the ortho/para H₂ ratio

L. Pagani¹ - C. Vastel² - E. Hugo³ - V. Kokoouline⁴ - C. H. Greene⁵ - A. Bacmann⁶ - E. Bayet⁷ - C. Ceccarelli⁶ - R. Peng⁸ - S. Schlemmer³

1 - LERMA & UMR8112 du CNRS, Observatoire de Paris, 61 Av. de l'Observatoire, 75014 Paris, France
2 - CESR, 9 avenue du colonel Roche, BP 44348, Toulouse Cedex 4, France
3 - I. Physikalisches Institut, Universität zu Köln, Zülpicher Strasse 77, 50937 Köln, Germany
4 - Department of Physics, University of Central Florida, Orlando, 32816, USA
5 - Department of Physics and JILA, University of Colorado, Boulder, Colorado 80309-0440, USA
6 - Laboratoire d'Astrophysique de Grenoble, Université Joseph Fourier, UMR 5571 du CNRS, 414 rue de la Piscine, 38041 Grenoble Cedex 09, France
7 - Department of physics and astronomy, University College London, Gower street, London WC1E 6BT, UK
8 - Caltech Submillimeter Observatory, 111 Nowelo Street, Hilo, HI 96720, USA

Received 14 July 2008 / Accepted 6 October 2008

Abstract
Context. The high degree of deuteration observed in some prestellar cores depends on the ortho-to-para H₂ ratio through the H₃⁺ fractionation.
Aims. We want to constrain the ortho/para H₂ ratio across the L183 prestellar core. This is required to correctly describe the deuteration amplification phenomenon in depleted cores such as L183 and to relate the total (ortho+para) H₂D⁺ abundance to the sole ortho-H₂D⁺ column density measurement.
Methods. To constrain this ortho/para H₂ ratio and derive its profile, we make use of the N₂D⁺/N₂H⁺ ratio and of the ortho-H₂D⁺ observations performed across the prestellar core. We use two simple chemical models limited to an almost totally depleted core description. New dissociative recombination and trihydrogen cation-dihydrogen reaction rates (including all isotopologues) are presented in this paper and included in our models.
Results. We estimate the H₂D⁺ ortho/para ratio in the L183 cloud, and constrain the H₂ ortho/para ratio: we show that it varies across the prestellar core by at least an order of magnitude, being still very high (0.1) in most of the cloud. Our time-dependent model indicates that the prestellar core is presumably older than 1.5-2  10⁵ years but that it may not be much older. We also show that it has reached its present density only recently and that its contraction from a uniform density cloud can be constrained.
Conclusions. A proper understanding of deuteration chemistry cannot be attained without taking into account the whole ortho/para family of molecular hydrogen and trihydrogen cation isotopologues as their relations are of utmost importance in the global scheme. Tracing the ortho/para H₂ ratio should also place useful constraints on the dynamical evolution of prestellar cores.

Key words: ISM: abundances - ISM: clouds - ISM: structure - astrochemistry - ISM: individual objects: L183 - molecular processes

1 Introduction

Studies of the earliest stages of star formation have increased in the last ten years with the advent of new receivers acquiring better spatial and spectral resolution. Prestellar cores are dense and cold cores where gravitational collapse has not yet occured. In the densest regions of the core (where is greater than a few 10⁴ cm^-3) most heavy molecules containing carbon, nitrogen and oxygen seem to deplete onto the dust grains and only light ions remain in the gas phase. There has been extensive observational evidence of CO and CS depletion in the center of prestellar cores (e.g. Tafalla et al. 2004; Pagani et al. 2005; Bacmann et al. 2002; Caselli et al. 1999; Bergin et al. 2002) which seems to be typical of the majority of dense cores. Nitrogen-bearing species like CN, NH₃ and N₂H⁺ appear to subsist longer before freezing-out onto the dust grains (Hily-Blant et al. 2008; Pagani et al. 2007; Tafalla et al. 2006).

In conditions under which heavy species are depleted, H⁺ and H₃⁺ (and its deuterated counterparts) are the most abundant ions subsisting in the gas phase. H₂D⁺ has been widely detected and mapped in protostars and prestellar cores (Caselli et al. 2003; Vastel et al. 2006; Caselli et al. 2008) through its ortho fundamental line, with abundances high enough to be explained by the CO depletion theory/observations. Although difficult to observe from Earth, the D₂H⁺ molecule has been detected with its para line towards the 16293E prestellar core in the L1689N molecular cloud (Vastel et al. 2004) with an abundance similar to the ortho-H₂D⁺ molecule as suggested by Phillips & Vastel (2003).

Consequently many theoretical studies were performed that included all the deuterated forms of the H₃⁺ ion (e.g. Roberts et al. 2003; Walmsley et al. 2004). However all nuclear spin states (ortho, meta, para, corresponding to the spin state of the protons or deuterons) must be taken into account in order to compare with the observational sets. Moreover the thermicity of the forward/backward reactions strongly depends on the symmetric state of the species. Though the influence of the ortho/para H₂ ratio on the chemistry of H₂D⁺ (Pagani et al. 1992) had been described several years before the first detection of the ion (Stark et al. 1999), it is only recently that this specific spin-state dependent chemistry has been studied in detail (Flower et al. 2004,2006a).

The motivation for our study is the many deuterated observations performed in the L183 prestellar core (PSC) and the main aim of this paper is to study the ortho/para H₂ ratio from some of these deuterated species observations. To this effect, we constrain two chemical models including all the symmetric states of H₂, D₂, H₃⁺ and its deuterated counterparts with observations of ortho-H₂D⁺, combined with previous N₂H⁺ and N₂D⁺ observations. These models have been set up using recent dissociative recombination rates computed for H₃⁺ and its isotopologues as well as all non negligible reaction rates between H₂ and H₃⁺ and their isotopologues (both presented for the first time in this paper).

2 Observations

2.1 Deuterated H₃⁺

We first observed the ortho-H₂D⁺ (and para-D₂H⁺) with the Caltech Submillimeter Observatory (CSO) monopixel receiver and subsequently took advantage of the newly installed Heterodyne ARray Program, 16 channel 350 GHz band (``B-band'') (HARP-B) camera on the James Clerk Maxwell Telescope (JCMT) to fully map its emission.

2.1.1 CSO observations

The H₂D⁺ and D₂H⁺ observations were carried out at the CSO, between August 2004 and April 2005, under good weather conditions (225 GHz zenith opacity always less than 0.06). Scans were taken, using the chopping secondary with a throw of 3, using the reference position:  = ,  = .

Table 1:   Line parameters from the JCMT and CSO observations. The positions are offsets to the dust peak emission at  = ,  = . For non-detected lines we give the 3 upper limit. For JCMT and p-D₂H⁺ at CSO, we give both the Monte Carlo (MC) and the LTE column density estimates.

The 345 GHz (respectively 650 GHz) sidecab receiver with a 50 MHz acousto-optical spectrometer backend was used for all observations with an average velocity resolution of 0.04 km s^-1 (respectively 0.02 km s^-1). At the observed frequencies of 372.421385(10) GHz for the H₂D⁺ (1₁₀-1₁₁) and 691.660483(20) for the D₂H⁺ (1₁₀-1₀₁) lines (Amano et al. 2005), the CSO 10.4 m antenna has a HPBW of about 20 and 11 respectively. We performed a cut in declination across the main dust peak in L183 between (0, -40 ) and (0, 130 ) for H₂D⁺ and concentrated on the reference position for D₂H⁺. The system temperature was typically 1000 to 2000 K for H₂D⁺ and 2500 to 3500 K for D₂H⁺.

The beam efficiency at 372 GHz (respectively 692 GHz) was measured on Saturn and Jupiter and found to be 60% (respectively 40%) for point sources. Pointing was monitored every 1.5 h and found to be better than 3 . If the emission is extended compared to the beam size of CSO, as appears to be the case for H₂D⁺ then the efficiency is about 70% at 372 GHz and 60% at 692 GHz. The data reduction was performed using the CLASS program of the GILDAS software developed at IRAM and the Observatoire de Grenoble and the LTE data analysis with CASSIS developed at CESR (http://www.cassis.cesr.fr).

These CSO observations confirmed that the H₂D⁺ emission was extended in this source as we could hypothesize from several tracers and we carried on the project at the JCMT, recently equipped with a new 16-pixel camera to be able to fully map the source in a reasonable amount of time.

2.1.2 JCMT observations

The JCMT observations were obtained during semester 07A in service mode, using the HARP-B 16 pixel camera (one pixel, located in a corner, was unavailable). A third of the observations was obtained in jiggle-chop mode and two thirds in position switch (PSw) mode. The jiggle-chop mode appeared to be no faster, the displacement of the telescope in PSw mode seeming minor compared to other overheads, and because the jiggle-chop mode works in the Nyquist regime, each pixel receives much less integration time than in PSw mode. As adjacent pixels had the same off spectrum subtracted, the spatial average did not bring much improvement and we subsequently changed to position switch mode because deep integration on weak signal appears more important than Nyquist sampling for this work. In PSw mode, we made 2  2 pointings to fill the gaps in the camera, thus achieving a full beam sampling. Two such sets were performed to cover the main dust peak and its northern extension (Fig. 2) with one pixel row overlap between the two.

Most of the observations were run in band 1 weather ( ) while a few were done in band 2 weather leading to rapid degradation of the system temperature. The source was observed only above 40 elevation and the band 1 weather system temperature was in the range 500-1000 K depending on the pixels and on the elevation. To observe both H₂D⁺ and N₂H⁺ (J:4-3) (at 372.672509 GHz), we tuned the receiver half way between the two lines and used a frequency resolution of 61 kHz so that the backend could cover both lines at once.

Data pre-reduction was done with the Starlink software (KAPPA, SMURF and STLCONVERT packages) and subsequently translated into CLASS format for final reduction.

2.2 CO depletion and dust content

All other observations used in this paper have been obtained and published previously. The dust content of L183 both in emission and in absorption has been described in Pagani et al. (2004,2003). The source size is half a degree and contains a long filament extending 6 from north to south. Two peaks are clearly visible, one just south of the middle of the filament (which we call the main peak) with an opacity of 150  and a second one, 3 north of the first one (the north peak) with an estimated opacity of 70 . These peaks have the characteristics of prestellar cores. Most of the filament have an opacity above 40 .

Two large scale C¹⁸O and C¹⁷O maps obtained with the Kitt Peak 12-m telescope fail to trace the dense filament (Pagani et al. 2005). It is now well established that this is due to depletion of CO onto grains. Surprisingly, the opacity at which the depletion begins (20 ) is two times higher than what is usually observed in other clouds (e.g. Alves et al. 1999; Bergin et al. 2002; Kramer et al. 1999) though it still appears at a density (3  10⁴ cm^-3) that is a typical depletion density threshold (Pagani et al. 2005). Possibly, the low density envelope where depletion has not yet occurred is very extended in this cloud (which is confirmed by its large influence on the C¹⁸O (J:2-1) line intensity (Pagani et al. 2002). The depletion factor in volume for CO has been estimated to be 43 on average (Pagani et al. 2005) on the line of sight of the main dust peak and is probably much higher in the inner part of this core where density is above 1  10⁵ cm^-3.

2.3 N₂H⁺ and N₂D⁺

N₂H⁺ and N₂D⁺ have been mapped at both low (Kitt Peak 12-m, Pagani et al. 2005) and high (IRAM 30-m) resolutions. From the high resolution data, a strip crossing the main dust peak has been extracted and published (Pagani et al. 2007). In that paper, we performed a detailed analysis of the N₂H⁺ and N₂D⁺ emission with a Monte Carlo model treating exactly the hyperfine structure and line overlap of these species. We derived several physical properties, namely a maximum density of 2  10⁶ cm^-3, with a radial dependence proportional to r^-1 up to 4000 AU and proportional to r^-2 beyond, a kinetic temperature of 7 (1) K, a slight depletion of N₂H⁺ in the inner 3000 AU of the core and a deuterium fractionation which is non-measurable at 10⁴ AU (<0.03) and reaches 0.7 (0.15) in the center. As far as we know this is the highest fractionation reported yet for a singly deuterated species. However this may not be exceptional when compared to the detection of triply deuterated species, like ND₃ (van der Tak et al. 2002; Lis et al. 2002) and CD₃OH (Parise et al. 2004) or to the fact that most reported fractionations are line-of-sight averages and are not derived from detailed profiles.

3 Analysis

We present 3 models in this paper: a Monte Carlo model to compute H₂D⁺ and D₂H⁺ line intensities, a chemical steady-state model and a chemical time-dependent model.

Compared to previous works, we benefit here from two new sets of coefficients and a large set of observations in a PSC. The H₃⁺ + H₂ (+isotopologues) set of rate coefficients are extracted from the PhD work of Hugo in advance of publication and the H₃⁺ (and deuterated counterparts) recombination rates have been computed for this work by Kokoouline and Greene and are presented in Appendix B. Rate coefficients, as computed by Hugo, describe all possible interactions between H₃⁺ and H₂ isotopologues, including reactive and non-reactive, elastic and inelastic collisional rates, while recombination coefficients describe the dissociative recombination (DR) rates of trihydrogen cation isotopologues. Both sets take into account all the ortho, para and (D₃⁺) meta states.

3.1 Line emission

We have analysed the line emission of both ortho-H₂D⁺ and para-D₂H⁺ using a two-level Monte Carlo code (adapted from Bernes 1979) with our new collisional coefficients. Because the temperature of the PSC is around 7 K for both the gas (Pagani et al. 2007) and the dust (Pagani et al. 2004,2003), the possibility of populating the next rotational level, = 2₁₂ at 114 K and 75 K above ground level for ortho-H₂D⁺ and para-D₂H⁺ respectively, is so low that it can be safely ignored. With the new coefficients, the critical densities are 1.1  10⁵ and 4.9  10⁵ cm^-3 respectively and thus the lines are close to thermal equilibrium in the inner core. In the case of ortho-H₂D⁺, the difference between LTE (2.0  0.25  10¹³ cm^-2) and Monte Carlo (2.3  0.25  10¹³ cm^-2) column density estimates is typically 10-15% in the direction of the dust peak.

The para-D₂H⁺ line has not been detected (see Caselli et al. 2008) and the 3 upper limit corresponds to a total column density of 2.4  10¹³ cm^-2 using the Monte Carlo code. The LTE estimate yields 1.5  10¹³ cm^-2.

3.2 Deuteration

H₃⁺ ions are formed at a rate 0.96  by cosmic ray ionization of H₂ (Walmsley et al. 2004), rapidly followed by reaction with another H₂ to form H₃⁺, and are destroyed in reactions with neutral species and in dissociative recombination with free electrons, negatively charged grains and possibly negatively charged polycyclic aromatic hydrocarbons (PAHs^-). In prestellar cores, the primary reservoirs of hydrogen and deuterium are H₂ and HD, with HD/H₂ = 2(D/H)   3.2  10^-5 (Linsky 2007). The proton exchanging reaction of H₃⁺ with H₂ partly regulates the H₂ ortho-to-para ratio but has no effect on the H₃⁺ abundance. Concurrently, the reaction with HD forms H₂D⁺ and this primary fractionation is then followed by the subsequent fractionations and produces D₂H⁺ and D₃⁺ (Phillips & Vastel 2003; Roberts et al. 2003):

The backward reactions are endothermic with an energy barrier of about 200 K (when considering only the ground level for each species) and were thought to be negligible at the low temperatures found in prestellar cores (20 K), in which case the abundance ratios would be greatly enhanced. However, such enhancement can be limited by various processes (see Fig. 4):

• dissociative recombination of H₃⁺ (and its deuterated counterparts) with free electrons or negatively charged grains (and PAHs^-?);

• reactions with ``proton-friendly'' molecules such as CO and N₂ which destroy the trihydrogen cations to produce HCO⁺ and N₂H⁺;

• ortho-H₂ which can react with the deuterated trihydrogen cation and remove the deuterium (see below).

In this modeling, we introduce the backward reactions to Eqs. (1)-(3) as we distinguish between ortho, meta and para states of the different species. When these reactions are completely neglected, the deuteration fractionation is considerably enhanced and observations towards pre-stellar cores cannot be reproduced (Roberts et al. 2003). Indeed, the importance of considering ortho and para states of various H/D carriers in the chemistry of the trihydrogen cation and isotopologues was first discussed by Pagani et al. (1992) and subsequently expanded in a series of papers by Flower and coworkers (Flower et al. 2006b,2004,2006a; Walmsley et al. 2004, hereafter collectively referred to as FPdFW). Not only is this important in comparing the chemical model predictions on the abundance of H₂D⁺ to the observations of the ortho-H₂D⁺ species alone but also because some important reactions are much faster with ortho-H₂ than with para-H₂ hence no longer negligible.. Indeed, the 170 K internal energy of the lowest ortho-H₂ level (J=1) is large with respect to the temperatures of interest and can significantly enhance the Boltzmann factor of endothermic reactions. In some cases, reactions which are endothermic with para-H₂ can turn out to be exothermic with ortho-H₂ i.e. fast and temperature independent. In fact, the ortho-to-para ratio of H₂ is found to be a crucial parameter for the entire chemistry of deuterium.

The key reactions involving ortho-H₂ are essentially with meta-D₃⁺, para-D₂H⁺ and ortho-H₂D⁺ as well as para-H₂D⁺ (see rates in Appendix A) because the internal energy of the ortho-H₂ alone is not enough to overcome the endothermicity of reactions 1 to 3. Thus only those species which also have an internal energy high enough (so that the sum of the two internal energies is higher than the endothermicity of reactions 1 to 3) can react with ortho-H₂ at the Langevin rate in cold gas. Thus, the three reactions involving ortho-H₂ with meta-D₃⁺, para-D₂H⁺ and ortho-H₂D⁺ present exothermic or thermoneutral channels to rehydrogenate the cations forming ortho-D₂H⁺, para-H₂D⁺ and ortho/para-H₃⁺ while the reaction between ortho-H₂ and para-H₂D⁺ can efficiently convert the latter to ortho-H₂D⁺. The ortho-H₂ molecule thus opens a path to climb the 4 step energy ladder back from para-D₂H⁺ to H₃⁺ via para-H₂D⁺ and ortho-H₂D⁺ which can be very efficient in the presence of large ortho-H₂ fractions. However, this efficient ladder scheme does not include D₃⁺ because the conversion of ortho-D₂H⁺ to para-D₂H⁺ is strictly forbidden in collisions with H₂ and very inefficient in collisons with HD. Nevertheless, these reactions can be a strong limit to the isotopic fractionation of H₃⁺ hence of other species. Any chemical model which includes deuterium chemistry must distinguish between ortho and para states of dihydrogen and trihydrogen cation isotopologues and include reactions between the different spin states following Flower et al. (2006b); Pagani et al. (1992); Flower et al. (2006a) and the present work.

3.3 CO and N₂ chemistry

The CO and N₂ abundances are critical parameters in the deuteration of the H₃⁺ ion. CO is expected to freeze-out onto the grain mantles at high densities (a few 10⁴ cm^-3) and low temperatures (20 K) (e.g. Tafalla et al. 2002; Pagani et al. 2005; Bacmann et al. 2002; Bergin et al. 2001; Caselli et al. 1999). With an N₂ binding energy similar to the CO binding energy (Bisschop et al. 2006; Öberg et al. 2005), these two molecules are expected to behave similarly. However, observations towards prestellar cores prove that N₂H⁺ (produced from N₂) remains in the gas phase at higher densities than CO. This can be explained by the fact that N₂H⁺ is mainly destroyed by CO (Aikawa et al. 2005; Pagani et al. 2005; Caselli 2002), so that the CO freeze-out implies a drop in the destruction rate of N₂H⁺. This would partially balance the lower formation rate due to the N₂ freeze-out. Consequently, N₂H⁺ is observed to survive in the gas phase at higher densities (10⁶ cm^-3). In the case of L183, we have shown that N₂H⁺ partially survives but suffers from depletion at densities starting at 5  10⁵ cm^-3 to reach a factor 6 ⁺¹³_-3 at the core centre (2  10⁶ cm^-3). Because of growing deuterium fractionation, N₂D⁺ abundance still increases towards the PSC center until the N₂ depletion becomes predominant over the deuterium enhancement, and in turn, the N₂D⁺ abundance slightly decreases in the inner most part of the core (Pagani et al. 2007).

The N₂D⁺ and N₂H⁺ ions can be produced via the following routes:

We assumed that all the H₃⁺ isotopologues react at the Langevin rate  with N₂ (which is inversely proportional to the square root of the reduced mass of the colliding system, hence to the mass of the H₃⁺ isotopologue) and that deuterium and hydrogen nuclei are equiprobably transferred. Consequently, H₃⁺, H₂D⁺, D₂H⁺ and D₃⁺ respectively produce an N₂D⁺:N₂H⁺ ratio of 0:3, 1:2, 2:1 and 3:0. The measured ratio of 0.7  0.15 in the center of L183 thus implies significant fractions of D₂H⁺ and D₃⁺. It has been shown (Walmsley et al. 2004) that in the case of complete depletion of heavy species (C, N, O ...), D₃⁺ would be the dominant trihydrogen cation isotopologue which would imply that N₂D⁺ is more abundant than N₂H⁺. This is not the case here; nevertheless the N₂H⁺ deuterium fractionation is a good constraint to the abundance of the four trihydrogen cation isotopologues in our chemical model.

At steady state (d[N₂H⁺]/dt = 0 and d[N₂D⁺]/dt = 0), reaction 4 and its isotopic variants (5 to 9) being the main path to produce N₂D⁺ and N₂H⁺, we obtain:

This ratio has been measured in the cut through the L183 main PSC. We describe in the following how our method can provide an estimate of the ortho/para H₂D⁺ ratio using this variable and subsequently, of the ortho/para H₂ ratio itself as well as some indication of the cosmic ionization rate and mean grain size.

3.4 Grain distribution

Recombination of ions with electrons on negatively charged grain surfaces is an important process since it can be much faster than in the gas phase, especially in the case of H⁺ (Draine & Sutin 1987). The negatively charged grain surface area is therefore a crucial parameter (we can safely ignore positively and multiply negatively charged grains, considered to be very rare in cold environments, Draine & Sutin 1987). The grain size distribution in prestellar cores is unknown since it mostly depends on grain condensation and also on ice condensation (e.g. Vastel et al. 2006). We thus treat the grain radius  as a parameter of the model, assuming all the grains to have the same size and the dust to gas mass ratio to be 0.01. Different values could be advocated (for example, in places where depletion is important, the ices increase the dust mass by up to 25%, Walmsley et al. 2004) but the net result is to change only slightly the average grain size which is not well constrained in the PSCs in any case.

The focusing effect of the Coulomb attraction between charged particles and oppositely charged grains has been included using the Draine & Sutin (1987) formalism:

where is the reduced temperature ( , e being the electron charge, k the Boltzman constant). Therefore the recombination rate of the H⁺ ion on a negatively charged grain can be expressed as:

where is the grain radius, is the hydrogen mass, T the kinetic temperature and S is the sticking coefficient (S  1). The latter represents the probability that a colliding species will stick onto the grain surface. For ions, Draine & Sutin (1987) concluded that the sticking coefficient should be unity. The same computation can be made to estimate the recombination rate of other ions, H₃⁺, H₂D⁺, D₂H⁺, D₃⁺, by a simple correction of the atomic mass of the ions (respectively /, /, /, /, ). In the case of collisions between charged particles and neutral grains, the attraction due to the polarization of the grain by the charged particle can be expressed through:

Therefore the sticking rate of electrons on neutral grains can be expressed as:

where is the electron mass and S is the sticking coefficient. S is about unity (Umebayashi & Nakano 1980) for a planar surface but a curvature of the grain surface will tend to reduce this parameter. However in the following we will use a factor of about unity as this parameter did not seem to have a large influence on the results in our runs.

The grain abundance [gr] can be expressed using:

where is the mean grain density (assumed to be 3 g cm^-3, is the dust-to-gas mass ratio, and is the mass of molecular hydrogen. Another important parameter in our model is the abundance of the negatively charged grains ( ). At steady-state, assuming partial depletion of CO and N₂ and total depletion of all the other heavy species:

(here we have neglected HD⁺, D₂⁺, He⁺, ...).

3.5 Steady-state chemical model

The code we describe in the following is used to calculate the steady-state abundances of the chemical species found in the different layers of the L183 prestellar core as listed in Table 2.

In the steady-state approximation the abundance species are interlinked via their production rates and their destruction rates (production = destruction).

Since H₃⁺ is produced at a rate 0.96, the H⁺ abundance can be expressed as (including only the main reactions):

The main production path of H₃⁺ is via cosmic ray ionization of H₂ and proceeds in two steps:

and H₂⁺ rapidly reacts with another H₂ to form H₃⁺ but the branching ratios between different combinations of spin states are non-trivial (Oka 2004):

These are different from those advocated by FPdFW who took branching ratios of 1/2 for both species. The ortho-H₃⁺ formation rate from cosmic ray ionization is therefore the sum of several terms:

The production rate for ortho-H₃⁺ can be expressed as (including only the main reactions. The rates are listed in Table A.1):

which represent respectively the formation from cosmic ray ionization, backward destruction of ortho-H₂D⁺ with ortho-H₂, spin conversion of para-H₃⁺ with ortho-H₂ and finally, spin conversion of para-H₃⁺ with HD.

Table 2:   Source parameters: distance from the core center, H₂ density, N₂H⁺ abundance and N₂D⁺/N₂H⁺ ratio. The position is measured away from the PSC center along the RA axis (from Pagani et al. 2007).

The destruction rate for ortho-H₃⁺ can be expressed as (including only the main reactions):

which respectively represents its destruction by dissociative recombination with electrons, spin conversion with ortho-H₂, spin conversion and deuteration with HD, proton transfer reactions with CO and N₂ and dissociative recombination on grains.

Similarly, the para-H₃⁺ formation from cosmic ray ionization can be expressed as:

The production rate for para-H₃⁺ is:

and the destruction rate is:

The N₂ abundance has been solved numerically to obtain the observed N₂H⁺ abundance. Electronic abundance is adjusted to reach equilibrium.

In our steady-state model, the H₂ ortho/para ratio, the average grain radius and the cosmic ionization rate  are the varying input parameters. Within each layer of the PSC model (Table 2), these parameters are adjusted to match the following:

• the H₂ density;

• the N₂D⁺/N₂H⁺ ratio at 7 K;

• the observed ortho-H₂D⁺ column density;

• the upper limit on the p-D₂H⁺ column density.

Although the full range of grain sizes and ortho-to-para H₂ ratios have been explored for each H₂ density, we have not allowed solutions in which, for example, the grain size would oscillate from one layer to the next. We have searched for solutions throughout the layers in two different scenarios: (i) the grain size and the ortho-to-para ratio were both kept constant; or (ii) the grain size was allowed to increase and the ortho-to-para H₂ ratio to decrease with the H₂ density. In both scenarios,  was kept constant throughout the layers. We neglected detailed reactions with D₂ as Flower et al. (2004) have shown that its role is negligible in general and we have kept the HD abundance constant which is generally a good approximation.

3.6 Time-dependent chemistry

In a second step, we have constructed a pseudo time-dependent model based on NAHOON, a chemical model, a version of which has been made publicly available by Wakelam. We have modified this model in two ways: 1) we have replaced electron (resp. ion) reactions with neutral (resp. charged) grains as provided in the Ohio State University (OSU) reaction file (delivered with NAHOON) by the set of equations described above (Sects. 3.4 and 3.5), which we have directly included in the program, to take into account Coulomb focusing; 2) we have included the formation of HD and D₂ on grain surfaces and we have introduced the spin state of H₂ and D₂ with the usual assumption that they are formed with the statistical ortho/para spin state ratio of 3 and 2 respectively. We have used the formation rate provided in the OSU reaction file (5  10^-17 cm³ s^-1) for the formation of molecular hydrogen. Because grains are covered by ice in the environments concerned here, we consider that the only interaction between the atoms and the surface is physisorption. In this case, the formation rates of HD and D₂ (in cm³ s^-1) is half for HD and 10⁵ times lower for D₂ with respect to H₂ formation (as dicussed in ). In environments where grains are not covered by icy mantles, on the other hand, one would have to consider chemisorption, which strongly changes the efficiencies of the formation of HD and D₂ (Cazaux et al. 2008). We have reduced the set of species and reactions to our needs, limiting ourselves to the most important reactions (see below) but differentiating all ortho and para (plus meta-D₃⁺) species as independent species and including all the detailed rates between the trihydrogen cation and dihydrogen isotopologues (see Sect. 3.7.1) including spin state conversions. However we have included more reactions than in the steady-state model, taking into account reactions with D₂, H₂⁺, He⁺, etc. and allowing the ortho/para H₂ ratio and the HD abundance to vary.

The main path to convert ortho-H₂ into para-H₂ is via the reaction

which proceeds seven orders of magnitude faster at 7 K than the reverse reaction.

3.7 Rate coefficients

Many groups have made available gas-phase rate coefficients. The University of Manchester Institute of Science and Technology (UMIST) Database for Astrochemistry contains information on 4500 reactions of which 35% have been measured experimentally, some at temperatures down to 20 K (Woodall et al. 2007). The OSU group provides approximately the same database but focuses more on low temperature chemistry. We accordingly use in our modeling some of the reactions in the latter (with the most recent version OSU2007), considering the low temperatures found in L183. Apart from the rates presented below, all reaction rates involving deuterium have been taken from FPdFW (except recombination on grains for which we have used different sticking probabilities). We disregard the odd branching ratio of N₂H⁺ dissociative recombination reported by Geppert et al. (2004) to consider a single possibility, namely the liberation of dinitrogen (Molek et al. 2007).

3.7.1 H₃⁺ + H₂ isotopologue reaction rates

Phase space theory (PST) was used to derive thermal state-to-state rate coefficients for the whole H₃⁺ + H H₃⁺ + H₂ system and isotopic variants in the temperature range 5-50 K. This statistical method accounts for such quantities as mass, energy, rotational angular momentum, nuclear spin symmetry and their respective conservation laws. The ergodic hypothesis which is a requisite for PST as well as the full-scrambling hypothesis are assumed according to the topology of the PES (Xie et al. 2005; Yamaguchi et al. 1987). Reactant (product) trajectories are treated with the classical Langevin model. The resulting set of state-to-state rate coefficients deviates from the detailed balance principle by a few percent at worst and is consistent with thermodynamical equilibrium constants. Details will be given in a forthcoming publication (Hugo et al., in prep.)

In the present astrochemical model, nuclear spin states of the different molecules are treated as distinct species but their rotational states are not considered individually. We thus made the assumption that only the rotational ground states of each nuclear spin species were populated and used the ground state-to-species thermal rate coefficients obtained by summing ground state-to-state thermal rate coefficients over the product channels.

3.7.2 H₃⁺ isotopologue dissociative recombination rates

The dissociative recombination (DR) rate of H₃⁺ is difficult to measure experimentally. Since the early measurements e.g. by Adams et al. (1984), numerous attempts have been made and are summarized in two papers (Fonseca dos Santos et al. 2007; Bates et al. 1993). In parallel, theoretical work has also been developed with the latest one published by Fonseca dos Santos et al. (2007). Extending upon that work, we present here, in Table B.1, the H₃⁺ updated DR rate (Fonseca dos Santos et al. 2007) along with newly calculated H₂D⁺, D₂H⁺, and D₃⁺ DR rates (see Appendix B for more details). These calculations do not predict the branching ratio of the DR products. We have thus adopted the branching ratios published elsewhere (Datz et al. 1995a; Zhaunerchyk et al. 2008; Strasser et al. 2004; Datz et al. 1995b) which we have applied to the calculated rates. The resulting DR rates at 7 K are listed in Table A.1. Several remarkable effects at low temperature are visible (see Figs. B.1-B.4):

• the strong departure of the ortho-H₃⁺ DR rate from that of para-H₃⁺. The difference is a factor of 10 at 10 K;

• D₂H⁺ shows a large DR rate drop, by a factor of 10 at 10 K for both ortho and para species compared to the extrapolated value used by FPdFW;

• on the contrary, a large increase of the D₃⁺ DR rate is predicted to occur but mostly at temperatures where deuteration is low and therefore the consequences for D₃⁺ abundance are limited.

4 Results and discussion

4.1 H₂D⁺ linewidth

In order to fit the observed ortho-H₂D⁺ line profiles, we run the Monte Carlo model using the ``best model'' velocity, density and temperature profiles derived from the N₂H⁺ and N₂D⁺ data analyzed in Pagani et al. (2007). However, the linewidth for the two stronger spectra (offsets (0, 0) and (15 , 0), Fig. 3) is too wide to be reproduced with the same micro-turbulent width which we have used for N₂H⁺ ( (FWHM)  0.14 km s^-1). Indeed, the central spectrum linewidth measured by fitting a Gaussian yields 0.45 (0.03) km s^-1. The thermal contribution is

at 7 K, which implies that (FWHM) contribution is 0.35 km s^-1, 2.5 times larger than for N₂H⁺. Here k is the Boltzmann constant and m is the mass of H₂D⁺. If we impose a turbulent velocity similar to the one modeled for N₂H⁺, then the temperature needed to obtain such a wide line would be 16 K which is completely ruled out by N₂H⁺ observations (Pagani et al. 2007). Infall motion limited to the inner core could have explained the H₂D⁺ width if no H₂D⁺ spectrum other than the central one had been wide combined with a large depletion of N₂H⁺. This is not possible because the H₂D⁺ spectrum at (15 , 0) has the largest width (0.49  0.03 km s^-1) in a region where N₂H⁺ is hardly depleted. This remains therefore a pending problem.

4.2 N₂H⁺ deuteration

4.2.1 Requested conditions

We next discuss the main parameters that control the N₂H⁺ deuteration using the steady-state chemical model.

The models have been run for a temperature of 7 K which prevails in all the layers where N₂D⁺ has been detected in the PSC cut presented in Pagani et al. (2007). We have also run the models for the corresponding density, N₂H⁺ abundance and N₂D⁺/N₂H⁺ ratio of each layer (the parameters are listed in Table 2).

As discussed above, the abundance of ortho-H₂ is the main controlling factor of the trihydrogen cation isotopologue abundances and therefore of the N₂D⁺/N₂H⁺ ratio (and similarly of the DCO⁺/HCO⁺ ratio, see e.g. Maret & Bergin 2007). We have therefore explored the range of possible solutions for the ortho/para H₂ ratio in the two extreme layers of our core profile (n(H₂) = 7  10⁴ cm^-3 and 2  10⁶ cm^-3) for which we have a N₂D⁺/N₂H⁺ ratio of <0.03 and 0.7  0.15 respectively (Fig. 5). We have done this for three cosmic ray ionization rates (10^-17, 10^-16, and 10^-15 s^-1) covering the values generally discussed in the literature (e.g. Maret & Bergin 2007) and for all possible electronic abundances (or average grain sizes as they are linked via the abundance of H⁺ which is mostly controlled by the grain surface area). In this first run, we have simulated total depletion by adjusting the CO and N₂ abundance to 10^-8. We have also indicated the range of N₂D⁺/N₂H⁺ ratio measured in both layers. The average grain radius has been varied from 0.01 m to 5 m and electronic abundance from 10^-11 to 10^-6 which cover the usually accepted values. We can see that N₂D⁺/N₂H⁺ ratios above 100 are possible in dense gas though they require very low electronic abundances and therefore very small grains which are probably absent from these dense and cold regions due to grain coagulation (see e.g. Stepnik et al. 2003).

In the lower density outer layer where no N₂D⁺ has been detected, the ortho/para H₂ ratio must be high enough, i.e. above 0.01, to prevent any deuteration from occuring whatever the cosmic ray ionization rate. On the contrary, the dense, strongly deuterated layer has solutions only below a maximum ortho/para H₂ ratio of 0.01 (or lower for high  rates). Thus the ortho/para H₂ ratio across the PSC clearly must vary from above 0.01 to below 0.01. In the case of low cosmic ray ionization rates (10^-17 s^-1), though the ortho/para H₂ ratio of 0.01 seems to be a common solution for both layers, it requires a high electronic abundance (and large grains) in the outer layer and a low electronic abundance (and small grains) in the inner dense part. This is clearly improbable. The temperature being low enough in all the layers, warm layers (above 20 K) cannot be invoked instead of a high ortho/para H₂ ratio to limit the deuteration in the outer parts of the PSC. CO total depletion is however questionable and in the model we also used a CO depletion factor of 10 (abundance of 10^-5) in the outer layer and a CO depletion factor of 100 (abundance of 10^-6) in the inner layer (Fig. 6). This only limits the maximum N₂D⁺/N₂H⁺ ratio which decreases by one order of magnitude. Indeed, the destruction of H₃⁺ by CO dominates over recombination with electrons when their abundance is very low and vice versa. However, the conditions to reach the observed N₂D⁺/N₂H⁺ ratio remain unchanged and therefore only a variable ortho/para H₂ ratio can be invoked. Such a variable ortho/para H₂ ratio cannot be investigated with a steady-state model because in all layers, the ortho-H₂ abundance would eventually decrease to values about of 10^-3-10^-4 as discussed by FPdFW.

4.2.2 The ortho/para H₂ variation

We discuss here the possibilities of making the ortho/para H₂ ratio vary across a single PSC.

It is commonly accepted that H₂ is formed on grain surfaces with an ortho/para ratio of 3 because of spin statistics and the exothermicity of the reaction H + H  H₂. Subsequently, the ortho-H₂ is converted into para-H₂ following Eq. (29) and to a lesser extent with reactions involving H₃⁺ and its isotopologues. As already discussed by FPdFW, this conversion is slow and has probably not reached steady state in clouds with ages between 10⁵ and 10⁶ years.

We have therefore run the modified Nahoon model to search for a simultaneous solution for all layers which would get the necessary ortho/para H₂ gradient across the PSC to achieve both the observed N₂D⁺/N₂H⁺ ratio and ortho-H₂D⁺ abundance profiles. In this model, we have no direct measurement of the CO abundance but for the credibility of the model, we have set the CO abundance to 10^-5 in the outer layer and increased the depletion to reach a factor of 100 (i.e. a CO abundance of 10^-6) in the core. The results are presented in Figs. 7 and 8. We have searched for a solution where each modeled layer meets the two observational constraints (N₂D⁺/N₂H⁺ ratio and ortho-H₂D⁺ abundance) at the same time, but these solutions are not simultaneous between the different layers. We varied the average grain radius and the cosmic ray ionization rate, . We could find solutions for grains of average radius 0.025 to 0.3 m. No solution has been found for grains above 0.3 m. Figures 7 and 8 show the case for which the grain average radius is 0.1 m and  = 2  10^-17 s^-1. In this case, the time range inside which all layers meet the requested conditions is 0.6 to 1.7  10⁵ years. Figure 8 shows how the ortho/para H₂ ratio evolves for 3 selected layers. We have marked the appropriate time which is the solution for each of these layers as established from Fig. 7. In that figure, we can see that the o/p H₂ ratio is below 0.01 for the dense layers and still above 0.01 for the outer layer for which no H₂D⁺ has been detected, as expected from the steady-state model. We can also see that the full o/p H₂ relaxation has not yet occurred even for the densest part of the cloud. Smaller grains have a larger interacting surface and therefore lower the abundance of H⁺ ions which preferentially recombine on negatively charged grains (or PAHs^-), consequently slowing down the dominant ortho-H₂ relaxation reaction (Eq. (29)). Though smaller grains also imply a lower electronic density, therefore favouring a higher deuteration of N₂H⁺ as shown in Fig. 5, the slower disappearance of ortho-H₂ is the dominant process here and finally, smaller grains slow down the deuteration process. For grains of average radius 0.025 m, the range of ages matching the range of N₂H⁺/N₂D⁺ observed ratios is 2.7-3.8  10⁵ years, while for grains of average radius 0.3 m, the time range is only 3.8-7.2  10⁴  years. Figure 8 suggests that D₂ should become a sizeable fraction of available deuterium a short while after the present state (typically 2-3  10⁵  years) and that HD correspondingly should drop slightly.

4.3 Age of the core and collapse

Though it is normal that dense layers evolve faster than less dense ones, at least to account for a differential ortho/para H₂ ratio, the layers are evolving too fast in our model. The densest layers would have reached their present status 2 to 3 times faster than the outer, less dense layers. This could be possible only if the denser layers had reached their steady-state equilibrium. This is not the case here where the densest layer would reach its steady-state equilibrium only after 2  10⁵  years and this would imply a N₂D⁺/N₂H⁺ ratio of 6, almost an order of magnitude higher than observed. The most probable reason for this time discrepancy is that the core has undergone a contraction and therefore all layers were not so dense in the past. While the outer layers have only little evolved in density (the most external one probably started at 0.5-1  10⁴ cm^-3 to reach 7  10⁴ cm^-3 today), the inner ones have undergone a much greater density increase. As constant density through the core would give no chemical differentiation while a time-frozen density profile as measured here gives too much differentiation, the solution is in between the two. Starting from a uniform gas, the chemical differential evolution of the core should therefore help us to constrain the duration of the contraction and the type of contraction. Of course the model should also include the evolution of depletion which also plays a role in the acceleration of the process.

As the core must have started to contract from a lower density region, typically 10⁴ cm^-3, it is clear that all layers have accelerated their chemical evolution while their density was increasing. Therefore, the layer with the longest time to reach the observed N₂D⁺/N₂H⁺ ratio gives a lower limit to the age of the cloud. Depending on the exact average size of the grains, this lower limit is 1.5-2  10⁵ years here. It is even larger because before the cloud underwent contraction, depletion had not yet occured and therefore species like atomic sulfur, S, must have been present in quantities large enough to transfer notable quantities of electric charges from H⁺ to S⁺ (H⁺ + S  H + S⁺) and PAHs must have also been abundant enough to help destroy H⁺ ions (Wakelam & Herbst 2008). All these phenomena contribute to the diminution of the H⁺ abundance, therefore slowing down the ortho-H₂ relaxation process. Indeed, Flower et al. (2006b) show that the relaxation process in some cases takes 3  10⁷  years, typically 50 times slower than in the case presented here (and 15 times slower for similar conditions of grain size and cosmic ray ionization rate but without depletion in their case).

4.4 Para-D₂H⁺

At 7 K, the strongest possible line intensity (LTE case) for the ground transition of para-D₂H⁺ is below 0.3 K because of the Rayleigh-Jeans correction at 691 GHz which becomes large. Moreover the thermalization of the line is difficult to obtain in this source beyond the radius of 3000 AU because the PSC density drops below the para-D₂H⁺ critical density (  = 4.9  10⁵ cm^-3) and a slight drop of the excitation temperature turns into an exponential decrease of the brightness temperature. For  = 6 K,   0.13 K. Therefore searches for para-D₂H⁺ must reach very low noise levels to have a chance of detection. From the models we present here, we predict an integrated line intensity of 11 mK km s^-1 (34 mK peak) with the Monte Carlo model, and an upper limit of 16 mK km s^-1 (48 mK peak) in the case of LTE. This is a factor of 3 to 4 below the upper limit we obtained from the observations.

4.5 The chemical profile

Finally, we obtain a detailed profile of the PSC which we present in Fig. 9. It represents the solution for the model we presented in Figs. 7 and 8 taking for each layer the values at their respective best-fit time. The large variation of the ortho-H₂ species across the core (a factor of 15) makes D₃⁺ change by an even larger amount (2 orders of magnitude) but it does not become the most abundant trihydrogen cation isotopologue in the core center at this stage because the ortho-H₂ abundance is not yet low enough. In the present case, for a density of 2  10⁶ cm^-3, the inversion between H₃⁺ and D₃⁺ occurs when the abundance of ortho-H₂ drops below 3  10^-3.

The ortho/para H₂ ratio in the outer layer is 0.04 (as expected from Fig. 6 which indicates a lower limit of 0.01). As discussed above, the ortho/para H₂ ratio evolution speed is linked to density and grain size, both of which are lower outside the PSC, in its embedding parental cloud. We can thus expect this ratio to be at least 0.05 and probably above 0.1 in the envelope of the cloud.

N₂ is an input parameter in our model because we have not included all the nitrogen chemistry. As reactions with H₃⁺ isotopologues are the main path to destroy this molecule (Flower et al. 2006b), we do not introduce a large error in determining its abundance directly from the N₂H⁺ abundance itself and obtain a N₂ profile which is probably closer to reality than if we had let the whole N-chemistry freely establish its abundance because of too many unknowns. The abundance profile thus starts at 1.5  10^-7 with respect to H₂ in the low density layer to diminish to 3  10^-8 in the densest layer. The undepleted N₂ abundance after attaining the steady state is 3  10^-5 (Flower et al. 2006b) but this is reached only after 5  10⁶ years. As depletion of N₂ in the outer layer is possibly still small, we can conclude that N₂ has not yet reached its steady state abundance which puts an upper limit on the age of the cloud of about 1  10⁵ years following the estimate of Flower et al. (2006b). However, this depends very much upon several factors, for example the C:O elemental abundance ratio. Consequently this information is only indicative.

Though our model does not deal with the nitrogen chemistry, it seems to indicate that low abundances of N₂ are sufficient to explain the observed N₂H⁺ abundance and therefore, the N₂ depletion seems not to be a critical factor as long as CO is also depleted. The much debated contradiction between the presence of N₂H⁺ in depleted cores while N₂ should deplete like CO would thus be a false problem. Low CO and electronic abundances, limiting the destruction rate of N₂H⁺, seem to be sufficient to compensate for the N₂ depletion itself to a large extent. Finally, in the inner core, N₂ and N₂H⁺ (+N₂D⁺) follow a similar decreasing trend, suggesting that N₂ depletion eventually forces N₂H⁺ to decrease.

While ortho-H₂D⁺ is 83 K above the para ground state, it is more abundant all over the core profile by almost an order of magnitude whereas the thermal equilibrium ratio would be ortho/para H₂D⁺ 2  10^-5 at 7 K. This demonstrates the efficiency of the ortho-H₂ in converting para-H₂D⁺ into ortho-H₂D⁺ and therefore limits the total abundance of H₂D⁺, as the backward channel to H₃⁺ remains open even at 7 K. This ortho/para population inversion does not occur for D₂H⁺ as the species needed to perform this inversion is no longer ortho-H₂ but the much rarer ortho-D₂. Therefore, the para-D₂H⁺ remains the least abundant of the two spin state species which, combined with the fact that its ground transition is higher in frequency than the one of ortho-H₂D⁺, makes its detection extremely difficult.

5 Conclusions

We have presented a pair of simple chemical models restricted to H-carriers, He plus CO and N₂ to account for the observed HCO⁺, DCO⁺ (not discussed in this paper), N₂H⁺, and N₂D⁺ ions. We have benefited from new computed reactions rates for both the H₃⁺ + H₂ isotopologue combinations and for the H₃⁺ isotopologue dissociative recombination rates which explicitly take into account the nuclear states individually.

With the steady-state model we have shown that the ortho/para ratio of H₂ must vary from above 0.01 in the outer parts of the L183 PSC to less than 0.01 in the inner parts to explain the variation of deuteration across the core. Checking the reality of the ortho/para H₂ variation with a time-dependent model, we have also found that if the present PSC density profile is static, then the inner layer would have reached its present status 2 to 3 times faster than the outer layers. Because the present status is not in steady-state, the layers should evolve at a similar rate and therefore the density must have been lower in the past. The most probable explanation is that the core has probably evolved from a uniform density cloud to the present centrally condensed PSC. The time-dependent model also suggests that the ortho/para H₂ ratio changes by one order of magnitude from 5% at a density of 7  10⁴ cm^-3 down to a few  10^-3 in the inner dense core. This has two important consequences:

• It is most probable that most of the cloud, outside the densest regions (i.e. the two PSC and the ridge in between) have an ortho/para H₂ ratio also above 5%, and possibly 10%, contrarily to what is usually assumed in models.

• In principle, it should be possible to fit the PSC profile with this chemical model combined with a dynamical model including depletion, to set an age to this PSC and possibly discriminate between several types of collapse but it is beyond the scope of this paper.

We already have some indication that the age of the PSC is somewhat above 1.5-2  10⁵ years though the N₂ abundance suggests a relatively short time (10⁵ years, except if depletion is compensating for its formation) when we compare our adjusted abundance to the formation rate of N₂ given by Flower et al. (2006b). The low abundance of N₂ needed to explain the observed N₂H⁺ abundance indicates that its depletion is not a real problem, though, obviously, N₂H⁺ would be much more abundant if N₂ was not being depleted but CO still was.

Finally, we stress the importance of considering ortho/meta/para chemistry when dealing with the deuteration of the interstellar medium. The importance of the ortho-H₂ on the amount of deuteration and the observations limited to only ortho-H₂D⁺ make this inclusion compulsory. Moreover, a complete state-to-state chemical model should be developed to take into account rotational pumping, leading to a higher destruction rate of the deuterated trihydrogen cation and possibly explaining the observed linewidth of ortho-H₂D⁺.

Detecting para-D₂H⁺ would be highly desirable to help constrain the models, but the high frequency and limited transparency of the atmosphere make it a difficult tool to use. Though observable from the ground, because of its weakness in cold dark clouds, which are the only places where it should be found, direct para-D₂H⁺ observations should be made on a large number of Galactic lines of sight using the HIFI receiver on board the Heschel Space Observatory.

Acknowledgements

We would like to thank D. Flower, G. Pineau des Forêts, M. Walmsley, S. Cazaux and V. Wakelam for fruitful discussions and an anonymous referee for her/his careful reading. Part of this work was supported by the National Science Foundation (NSF) grant AST 05-40882 to the CSO. The authors are grateful to the CSO and JCMT staffs for their support at the telescopes. We would like to thank Samantha Santos for the help provided in numerical calculations of thermal DR rate coefficients. This work has been supported by the NSF under Grants No. PHY-0427460 and PHY-0427376, by an allocation of NERSC and NCSA (project # PHY-040022) supercomputing resources. This work has benefited from research funding from the European Community's Sixth Framework Programme under RadioNet contract R113CT 2003 5158187.

References

• Adams, N. G., Smith, D., & Alge, E. 1984, J. Chem. Phys., 81, 1778 [CrossRef] (In the text)
• Aikawa, Y., Herbst, E., Roberts, H., & Caselli, P. 2005, ApJ, 620, 330 [NASA ADS] [CrossRef]
• Alves, J., Lada, C. J., & Lada, E. A. 1999, ApJ, 515, 265 [CrossRef]
• Amano, T., & Hirao, T. 2005, J. Mol. Spectrosc., 233, 7 [NASA ADS] [CrossRef] (In the text)
• Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2002, A&A, 389, L6 [NASA ADS] [CrossRef] [EDP Sciences]
• Bates, D. R., Guest, M. F., & Kendall, R. A. 1993, Planet. Space Sci., 41, 9 [NASA ADS] [CrossRef]
• Bergin, E. A., Ciardi, D. R., Lada, C. J., Alves, J., & Lada, E. A. 2001, ApJ, 557, 209 [NASA ADS] [CrossRef]
• Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101 [NASA ADS] [CrossRef]
• Bernes, C. 1979, A&A, 73, 67 [NASA ADS] (In the text)
• Bisschop, S. E., Fraser, H. J., Öberg, K. I., van Dishoeck, E. F., & Schlemmer, S. 2006, A&A, 449, 1297 [NASA ADS] [CrossRef] [EDP Sciences]
• Bunker, P., & Jensen, P. 1998, Molecular Symmetry and Spectroscopy (NRC Research Press) (In the text)
• Caselli, P. 2002, Planet. Space Sci., 50, 1133 [NASA ADS] [CrossRef]
• Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C. 1999, ApJ, 523, L165 [CrossRef]
• Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, 37 [NASA ADS] [CrossRef]
• Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 707 [NASA ADS]
• Cazaux, S., Caselli, P., Cobut, V., & Le Bourlot, J. 2008, A&A, 483, 495 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Chang, E. S., & Fano, U. 1972, Phys. Rev. A, 6, 173 [NASA ADS] [CrossRef]
• Chang, Q., Cuppen, H. M., & Herbst, E. 2006, A&A, 458, 497 [NASA ADS] [CrossRef] [EDP Sciences]
• Datz, S., Larsson, M., Stromholm, C., et al. 1995a, Phys. Rev. A, 52, 2901 [NASA ADS] [CrossRef]
• Datz, S., Sundström, G., Biedermann, Ch., et al. 1995b, , 74, 896
• Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 [NASA ADS] [CrossRef] (In the text)
• Flower, D. R., Pineau des Forêts, G., & Walmsley, C. M. 2004, A&A, 427, 887 [NASA ADS] [CrossRef] [EDP Sciences], FPdFW
• Flower, D. R., Pineau des Forêts, G., & Walmsley, C. M. 2006a, A&A, 449, 621 [NASA ADS] [CrossRef] [EDP Sciences], FPdFW
• Flower, D. R., Pineau des Forêts, G., & Walmsley, C. M. 2006b, A&A, 456, 215 [NASA ADS] [CrossRef] [EDP Sciences], FPdFW
• Fonseca dos Santos, S., Kokoouline, V., & Greene, C. H. 2007, J. Chem. Phys., 127, 124309 [NASA ADS] [CrossRef]
• Geppert, W. D., Thomas, R., Semaniak, J., et al. 2004, ApJ, 609, 459 [NASA ADS] [CrossRef] (In the text)
• Gerlich, D. 1990, J. Chem. Phys., 92, 2377 [NASA ADS] [CrossRef] (In the text)
• Harju, J., Juvela, M., Schlemmer, S., et al. 2008, A&A, 482, 535 [NASA ADS] [CrossRef] [EDP Sciences]
• Hily-Blant, P., Walmsley, M., Pineau des Forêts, G., & Flower, D. 2008, A&A, 480, L5 [NASA ADS] [CrossRef] [EDP Sciences]
• Kokoouline, V., & Greene, C. H. 2003a, , 90, 133201
• Kokoouline, V., & Greene, C. H. 2003b, Phys. Rev. A, 68, 012703 [NASA ADS] [CrossRef]
• Kokoouline, V., & Greene, C. H. 2004, Faraday Discuss., 127, 413 [CrossRef]
• Kokoouline, V., & Greene, C. H. 2005, Phys. Rev. A, 72, 022712 [NASA ADS] [CrossRef]
• Kokoouline, V., Greene, C. H., & Esry, B. D. 2001, Nature, 412, 891 [NASA ADS] [CrossRef] (In the text)
• Kramer, C., Alves, J., Lada, C. J., et al. 1999, A&A, 342, 257
• Kreckel, H., Motsch, M., Mikosch, J., et al. 2005, , 95, 263201 (In the text)
• Landau, L., & Lifshitz, E. 2003, Quantum Mechanics: Non-relativistic Theory (Burlington MA: Butterworth Heinemann) (In the text)
• Larsson, M., Danared, H., Larson, Å., et al. 1997, , 79, 395 (In the text)
• Linsky, J. L. 2007, Space Sci. Rev., 130, 367 [NASA ADS] [CrossRef] (In the text)
• Lipshtat, A., Biham, O., & Herbst, E. 2004, MNRAS, 348, 1055 [NASA ADS] [CrossRef]
• Lis, D. C., Roueff, E., Gerin, M., et al. 2002, ApJ, 571, L55 [NASA ADS] [CrossRef]
• Maret, S., & Bergin, E. A. 2007, ApJ, 664, 956 [NASA ADS] [CrossRef] (In the text)
• Mistrik, I., Reichle, R., Helm, H., & Muller, U. 2001, Phys. Rev. A, 63 (In the text)
• Molek, C. D., McLain, J. L., Poterya, V., & Adams, N. G. 2007, J. Phys. Chem. A, 111, 6760 [CrossRef] (In the text)
• Öberg, K. I., van Broekhuizen, F., Fraser, H. J., et al. 2005, ApJ, 621, L33 [NASA ADS] [CrossRef]
• Oka, T. 2004, J. Mol. Spectrosc., 228, 635 [NASA ADS] [CrossRef] (In the text)
• Pagani, L., Salez, M., & Wannier, P. G. 1992, A&A, 258, 479 [NASA ADS] (In the text)
• Pagani, L., Pardo, J.-R., Fich, M., Motte, F., & Stepnik, B. 2002, in Infrared and Submillimeter Space Astronomy, ed. M. Giard et al., held in Toulouse, June 2001 (Les Ulis: EDP Sciences), EAS Ser., 4, 145 (In the text)
• Pagani, L., Lagache, G., Bacmann, A., et al. 2003, A&A, 406, L59 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Pagani, L., Bacmann, A., Motte, F., et al. 2004, A&A, 417, 605 [NASA ADS] [CrossRef] [EDP Sciences]
• Pagani, L., Pardo, J.-R., Apponi, A. J., Bacmann, A., & Cabrit, S. 2005, A&A, 429, 181 [NASA ADS] [CrossRef] [EDP Sciences]
• Pagani, L., Bacmann, A., Cabrit, S., & Vastel, C. 2007, A&A, 467, 179 [NASA ADS] [CrossRef] [EDP Sciences]
• Parise, B., Castets, A., Herbst, E., et al. 2004, A&A, 416, 159 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Phillips, T. G., & Vastel, C. 2003, in Chemistry as a Diagnostic of Star Formation, ed. C. L. Curry, & M. Fich (Ottawa: NRC Press), 3 (In the text)
• Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41 [NASA ADS] [CrossRef]
• Stark, R., van der Tak, F. F. S., & van Dishoeck, E. F. 1999, ApJ, 521, L67 [CrossRef] (In the text)
• Stepnik, B., Abergel, A., Bernard, J.-P., et al. 2003, A&A, 398, 551 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Strasser, D., Lammich, L., Kreckel, H., et al. 2004, Phys. Rev. A, 69, 064702 [NASA ADS] [CrossRef]
• Tafalla, M., Myers, P. C., Caselli, P., Walmsley, C. M., & Comito, C. 2002, ApJ, 569, 815 [NASA ADS] [CrossRef]
• Tafalla, M., Myers, P. C., Caselli, P., & Walmsley, C. M. 2004, ApJ, 416, 191 [NASA ADS]
• Tafalla, M., Santiago-García, J., Myers, P. C., et al. 2006, A&A, 455, 577 [NASA ADS] [CrossRef] [EDP Sciences]
• Umebayashi, T., & Nakano, T. 1980, PASJ, 32, 405 [NASA ADS] (In the text)
• van der Tak, F. F. S., Schilke, P., Müller, H. S. P., et al. 2002, A&A, 388, L53 [NASA ADS] [CrossRef] [EDP Sciences]
• Vastel, C., Phillips, T. G., & Yoshida, H. 2004, ApJ, 606, L127 [NASA ADS] [CrossRef] (In the text)
• Vastel, C., Caselli, P., Ceccarelli, C., et al. 2006, ApJ, 645, 1198 [NASA ADS] [CrossRef]
• Wakelam, V., & Herbst, E. 2008, ApJ, 680, 371 [NASA ADS] [CrossRef] (In the text)
• Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035 [NASA ADS] [CrossRef] [EDP Sciences], FPdFW
• Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Xie, Z., Braams, B. J., & Bowman, J. M. 2005, J. Chem. Phys., 122, 224307 [NASA ADS] [CrossRef]
• Yamaguchi, Y., Gaw, J. F., Remington, R. B., & Schaefer III, H. F. 1987, J. Chem. Phys., 86, 5072 [NASA ADS] [CrossRef]
• Zhaunerchyk, V., Thomas, R. D., Geppert, W. D., et al. 2008, Phys. Rev. A, 77, 034701 [NASA ADS] [CrossRef]

Online Material

Appendix A: The reaction rate table used in the Nahoon modified chemical model

Table A.1:   Reactions used in the model. The rate coefficients are given for 7 K. Reaction rates less than 10^-15 cm^-3 s^-1 are not taken into account in our models. Reference (1) corresponds to Gerlich (1990), references (2)-(4) correspond to Hugo, OSU 07, and this paper respectively. For OSU 07, branching ratios involving spin states have been infered from quantum mechanical rules. For reactions involving grains, a grain radius of 0.1 m and a sticking coefficient of 1 have been considered. (5) Datz et al. (1995b) (6) Datz et al. (1995a), (7) Zhaunerchyk et al. (2008), (8) Larsson et al. (1997), (9) Molek et al. (2007).

Appendix B: On the rate coefficients for dissociation recombination of H₃⁺ and its isotopologues

The dissociative recombination (DR) rate coefficients for ortho- and para-H₃⁺ have been published by Fonseca dos Santos et al. (2007). Here, we present the results obtained for all four H₃⁺ isotopologues. The DR rate coefficients for different species of the nuclear spin are calculated using the approach described in a series of papers devoted to DR theory for triatomic molecular ions. See Fonseca dos Santos et al. (2007); Kokoouline & Greene (2003a,b) for H₃⁺ and D₃⁺ calculations and Kokoouline & Greene (2004,2005) for H₂D⁺ and D₂H⁺. The scope of this paper does not allow us to review the theoretical approach in detail. We only list its main ingredients.

The theoretical approach is fully quantum mechanical and incorporates no adjustable parameters. It relies on ab initio calculations of potential surfaces for the ground electronic state of the H₃⁺ ion and several excited states of the neutral molecule H₃, performed by Mistrik et al. (2001).

The total wave function of the system is constructed by an appropriate symmetrization of products of vibrational, rotational, electronic, and nuclear spin factors. Therefore, rovibronic and nuclear spin degrees of freedom are explicitly taken into account.

The electronic Born-Oppenheimer potentials for the four H₃⁺ (and H₃) isotopologues have the C_3v symmetry group. The C_3v symmetry group has a two-dimensional irreducible representation E. The ion has a closed electronic shell. The lowest electronic state of the outer electron in H₃ has p-wave character. The p-wave state of the electron also belongs to the E representation. Due to the Jahn-Teller theorem (Landau & Lifshitz 2003), this leads to a strong non-adiabatic coupling between the E-degenerate vibrational modes of the ion and the p-wave states of the incident electron. The coupling is responsible for the fast DR rate (Kokoouline et al. 2001) in H₃⁺. In the present model, only the p-wave electronic states are included because other partial waves have a much smaller effect on the DR probability: the s-wave states have no E-type character and, therefore, are only weakly coupled to the dissociative electronic states of H₃; d-wave electronic states are coupled to the E-vibrational modes, but the coupling is rather small because the d-wave of the incident electron does not penetrate sufficiently close to the ionic core owing to the d-wave centrifugal potential barrier.

All three internal vibrational coordinates are taken into account. Vibrational dynamics of the ionic core are described using the hyper-spherical coordinates, which represent the three vibrational degrees of freedom by a hyperradius and two hyperangles. The hyperradius is treated as a dissociation coordinate that represents uniformly the two possible DR channels, the three-body (such as H+H+H) and two-body (such as H₂+H). Although the initial vibrational state of the ion is the ground state, after recombination with the electron, other vibrational states of the ionic target molecule can be populated. Therefore, in general, many vibrational states have to be included in the treatment. In particular, the states of the vibrational continuum have to be included, because only such states can lead to the dissociation of the neutral molecule. The vibrational states of the continuum are obtained using a complex absorbing potential placed at a large hyperradius to absorb the flux of the outgoing dissociative wave.

Since the rovibrational symmetry is D_3h for H₃⁺ and D₃⁺ and C_2v for H₂D⁺ and D₂H⁺, the rovibrational functions are classified according to the irreducible representations of the corresponding symmetry groups, i.e. A₁', A₁'', A₂', A₂'', E', and E'' for D_3h and A₁, A₂, B₁, and B₂ for C_2v. We use the rigid rotor approximation, i.e. the vibrational and rotational parts of the total wave function are calculated independently by diagonalizing the corresponding Hamiltonians. In our approach, the rotational wave functions must be obtained separately for the ions and for the neutral molecules. They are constructed in a different way for the D_3h and C_2v cases. The rotational eigenstates and eigenenergies of the D_3h molecules are symmetric top wave functions (see, for example Bunker & Jensen 1998). They can be obtained analytically if the rotational constants are known. The rotational constants are obtained numerically from vibrational wave functions, i.e. they are calculated separately for each vibrational level of the target molecule. The rotational functions for the C_2v ions are obtained numerically by diagonalizing the asymmetric top Hamiltonian (Kokoouline & Greene 2005; Bunker & Jensen 1998).

Once the rovibrational wave functions are calculated, we construct the electron-ion scattering matrix (S-matrix). The S-matrix is calculated in the framework of quantum defect theory (QDT) (see, for example, Chang & Fano 1972; Kokoouline & Greene 2005,2003b) using the quantum defect parameters obtained from the ab initio calculation (Mistrik et al. 2001). The constructed scattering matrix accounts for the Jahn-Teller effect and is diagonal with respect to the different irreducible representations  and the total angular momentum N of the neutral molecule. Thus, the actual calculations are made separately for each and N. Elements of the matrix describe the scattering amplitudes for the change of the rovibrational state of the ion after a collision with the electron. However, the S-matrix is not unitary due to the presence of the dissociative vibrational channels (i.e. continuum vibrational states of the ion, discussed above), which are not explicitly listed in the computed S-matrix. The ``defect'' from unitarity of each column of this S-matrix is associated with the dissociation probability of the neutral molecule formed during the scattering process. The dissociation probability per collision is then used to calculate the DR cross-sections and rate coefficients.

The nuclear spin states are characterized by one of the A₁, A₂, or E irreducible representations of the symmetry group S₃ for D_3h molecules and by the A or B irreducible representations of the symmetry group S₂ for C_2v molecules. The irreducible representation  of a particular nuclear spin state determines its statistical weight and is related to the total nuclear spin  of the state. Here,  is the vector sum of spins  of identical nuclei.

For H₃⁺, the  states (A₂' and A₂'' rovibrational states) correspond to I=3/2 (ortho); the  states (E' and E'' rovibrational states) correspond to I=1/2 (para). The statistical ortho:para weights are 2:1.

For H₂D⁺, the states (B₁ and B₂ rovibrational states) correspond to I=1 (ortho); the  states (A₁ and A₂ rovibrational states) correspond to I=0 (para). The statistical ortho:para weights are 3:1.

For D₂H⁺, the states (A₁ and A₂ rovibrational states) correspond to I=0,2 (ortho); the  states (B₁ and B₂ rovibrational states) correspond to I=1 (para). The statistical ortho:para weights are 2:1.

Finally, for D₃⁺, the states (A₁' and A₁'' rovibrational states) correspond to I=1,3 (ortho); the  states (A₂' and A₂'' rovibrational states) correspond to I=0 (para); the  states (E' and E'' rovibrational states) correspond to I=1,2 (meta). The statistical ortho:para:meta weights are 10:1:8.

Figures (B.1)-(B.4) summarize the obtained DR thermal rate coefficients calculated separately for each nuclear spin species of the four H₃⁺ isotopologues and the numerical values are listed in Table B.1. For comparison, the figures also show the analytical dependences used in previous models of prestellar core chemistry (FPdFW). As one can see, the rates for different nuclear spin species are similar to each other (for a given isotopologue) at high temperatures. However, for lower temperatures, the rates for different ortho/para/meta-nuclear spin species significantly differ from each other. The difference in behavior at low temperatures is explained by different energies of Rydberg resonances present in DR cross-sections at low electron energies. The actual energies of such resonances are important for the thermal average at temperatures below or similar to the energy difference between ground rotational levels of different nuclear spin species. At higher temperatures, the exact energy of the resonances is not important. The averaged rate is determined by the density and the widths of the resonances, which are similar for all nuclear spin species over a large range of collision energies.

Table B.1:   Dissociative recombination rates of H₃⁺, H₂D⁺, D₂H⁺, and D₃⁺ for each individual nuclear spin state species.

Footnotes

... ratio

Appendices A and B are only available in electronic form at http://www.aanda.org

... jiggle-chop

http://www.jach.hawaii.edu/software/jcmtot/het_obsmodes.html

... Wakelam

http://www.obs.u-bordeaux1.fr/radio/VWakelam/Valentine%20 Wakelam/Downloads.html

... abundance

Here, we look for general solutions of the N₂D⁺/N₂H⁺ ratio as a function of several parameters, we thus do not try to fit the N₂ abundance to obtain the observed N₂H⁺ abundance.

All Tables

Table 1:   Line parameters from the JCMT and CSO observations. The positions are offsets to the dust peak emission at  = ,  = . For non-detected lines we give the 3 upper limit. For JCMT and p-D₂H⁺ at CSO, we give both the Monte Carlo (MC) and the LTE column density estimates.

Table 2:   Source parameters: distance from the core center, H₂ density, N₂H⁺ abundance and N₂D⁺/N₂H⁺ ratio. The position is measured away from the PSC center along the RA axis (from Pagani et al. 2007).

Table A.1:   Reactions used in the model. The rate coefficients are given for 7 K. Reaction rates less than 10^-15 cm^-3 s^-1 are not taken into account in our models. Reference (1) corresponds to Gerlich (1990), references (2)-(4) correspond to Hugo, OSU 07, and this paper respectively. For OSU 07, branching ratios involving spin states have been infered from quantum mechanical rules. For reactions involving grains, a grain radius of 0.1 m and a sticking coefficient of 1 have been considered. (5) Datz et al. (1995b) (6) Datz et al. (1995a), (7) Zhaunerchyk et al. (2008), (8) Larsson et al. (1997), (9) Molek et al. (2007).

Table B.1:   Dissociative recombination rates of H₃⁺, H₂D⁺, D₂H⁺, and D₃⁺ for each individual nuclear spin state species.

All Figures

---

Copyright ESO 2009