© ESO, 2017

1. Introduction

The ion is a universal proton donor and has therefore a central role in ion-molecule chemistry (Herbst & Klemperer 1973). In cold interstellar clouds, where the abundance of CO is reduced owing to its accretion onto dust, becomes the most abundant cation and the principal distributor of deuterium from HD to other species (Dalgarno & Lepp 1984; Roberts et al. 2003) in reactions such as (1)which is exothermic by 232 K when the reactants and products lie in their ground (para) states. The ion also contributes strongly to the ortho-para conversion of H₂ in these regions (Flower et al. 2006; Pagani et al. 2009). (In what follows, the ortho and para states of each species are simply referred to as o and p, respectively.) Spin states play an important role in the development of deuterium chemistry, because the high energy of the rotational ground state of oH₂ compared to that of pH₂ (~ 170 K) can cause reaction (1) to proceed in the backward direction even at low temperatures, transferring deuterium back to HD.

The isotopologs are among the most important tracers of the high-density and low-temperature regions of pre-stellar cores, where species heavier than He may be highly depleted because of freeze-out (Caselli et al. 2003; Walmsley et al. 2004; Friesen et al. 2014). Also, a good understanding of spin state chemistry is fundamental because the spin ratios (in particular of H₂D⁺) and the linked deuteration fraction measured in, for example, N₂H⁺, have been shown to be sensitive to the chemical age of dense clouds (Brünken et al. 2014; Kong et al. 2015; Pagani et al. 2011). An accurate derivation of the chemical age of a cloud requires accurate chemical codes.

Investigation of chemical reaction dynamics is being pursued through quantum mechanical scattering calculations and state-to-state experiments (e.g., Teslja & Valentini 2006; Zhang & Guo 2016). Besides providing microscopic characterization of collisions between molecules, some of the results of these endeavors are directly applicable to astrophysical environments. One such work is the study of the isotopic system by Hugo et al. (2009). The reaction between the trihydrogen cation, , and the hydrogen molecule, H₂, is of fundamental importance for the physics and chemistry of interstellar clouds.

Hugo et al. (2009) derived state-to-state thermal rate coefficients for inelastic and reactive collisions between all the isotopic variants of and H₂ and their different spin modifications. They compiled a table of ground-state-to-species rate coefficients for reactive collisions, where it is assumed that all reacting isotopologs of are in their ground rotational states. These data have been used in many studies of the chemistry of the reacting system (e.g., Pagani et al. 2009; Sipilä et al. 2010; Albertsson et al. 2014; Le Gal et al. 2014; Furuya et al. 2015; Majumdar et al. 2017). The state-to-state rate coefficients can also be used to estimate species-to-species rate coefficients, where possible rotational excitation of the reactant species is explicitly taken into account. The assumption that the reactant species lie in their ground states is likely to be valid in the dense starless cores of molecular clouds with average densities of up to n(H₂) ~ 10⁵ cm^-3, but in pre-stellar and star-forming cores with central densities exceeding 10⁶ cm^-3, some of the lowest rotationally excited levels of , H₂D⁺, D₂H⁺, and should be excited and contribute to the total species-to-species reaction rates. Besides the construction of rate coefficients for chemical reactions, the state-to-state coefficients given by Hugo et al. (2009) are needed for calculating the populations of the rotationally excited levels of H₂D⁺ and other isotopologs in connection with radiative transfer calculations.

A discrepancy between the predictions of chemical models utilizing either the ground-state-to-species or species-to-species rate coefficients was recently suggested to arise in the protostellar core IRAS 16293-2422, where the model using the ground-state-to-species rate coefficients seems to overpredict the oH₂D⁺ abundance (Harju et al. 2017). In the present paper, we discuss a chemical model which makes an effort to correct this effect by calculating the species-to-species rate coefficients from the state-to-state coefficients of Hugo et al. (2009), considering the density of the gas where the model is to be applied. This is done by estimating the populations of excited rotational states of , H₂D⁺, D₂H⁺, and at different densities. For H₂D⁺ and D₂H⁺, the populations are obtained through radiative transfer modeling. For the symmetric ions and , which have no electric dipole moments, the populations are assumed to obey the Boltzmann distribution.

The paper is organized as follows. We describe our models in Sect. 2. In Sect. 3 we present our results, which are further discussed in Sect. 4. We give our conclusions in Sect. 5. Appendices A and B contain additional discussion on critical densities, and rate coefficient tables.

2. Model description

2.1. Chemical model

In this paper we aim to compare the abundances of the spin states of the isotopologs as calculated with a model using either newly calculated species-to-species rate coefficients or the previously used ground-state-to-species rate coefficients. We employ the gas-grain chemical model presented in detail in Sipilä et al. (2015a,b). Unless otherwise noted, we use the same physical parameters and initial abundances as given in Tables 1 and 3 in Sipilä et al. (2015a). Here we use the KIDA gas-phase network (Wakelam et al. 2015) as the basis upon which the deuterium and spin-state chemistry is added according to the procedures discussed in detail in Sipilä et al. (2015a,b).

2.2. New fits to the ground state-to-species rate coefficients

Hugo et al. (2009) calculated the state-to-state rate coefficients for the reacting system in the temperature range 5–50 K. The data shown in their Table VIII assumes that the reactants are in their rotational ground states. These so-called ground-state-to-species rate coefficients were produced by fitting a two-parameter curve of the form αexp(−γ/T) in the temperature range 5–20 K¹ to data obtained from the state-to-state calculations. Later, Pagani et al. (2013) extended the fit to the range 5–50 K for all included reactions using the modified Arrhenius rate law k = α (T/ 300)^β exp(−γ/T).

Because the species-to-species rate coefficients (see Sect. 2.3) pertain to the temperature range 5–50 K, for consistency we made a new fit to the ground-state-to-species rates calculated by Hugo et al. (2009) using the modified Arrhenius rate law as Pagani et al. (2013) did. The resulting rate coefficients are given in Table B.1. Unlike Pagani et al. (2013), we did not correct any of the fits by hand, which leads to different values for the rate coefficients in some cases. We point out that there is a mistake in Table VIII in Hugo et al. (2009). The reaction is erroneously marked as forbidden, while it is in fact allowed by spin selection rules (as also confirmed by the Hugo et al. state-so-state data). However, this reaction is insignificant as its rate coefficient is of the order of 10^-18 cm³ s^-1 at 10 K.

Figure 1 shows a comparison of the abundances of the spin states of the isotopologs calculated with the ground-state-to-species rate coefficient fit from Table VIII in Hugo et al. (2009), and our new fit. Evidently the total abundances of the isotopologs are unaffected by the choice of rate coefficients, while the spin-state abundance ratios present slight deviations for the doubly and triply deuterated species (the variations for and H₂D⁺ are not clearly visible in the plot). Calculations at different temperatures yield the same conclusions.

Fig. 1 Abundances (upper row) and abundance ratios (lower row) of the various isotopologs as functions of time. The left-hand panels correspond to n(H2) = 105 cm-3, while the right-hand panels correspond to n(H2) = 106 cm-3. The temperature is set to Tgas = Tdust = 10 K. Solid lines represent calculations using the Hugo et al. (2009) ground-state-to-species rate coefficients, while dashed lines represent calculations using our new fit to the same coefficients (see text).

2.3. Species-to-species rate coefficients

Given that H₂D⁺ emission is observed to be ubiquitous (low-mass starless, pre-stellar, and protostellar cores, Caselli et al. 2008; high-mass star forming regions, Pillai et al. 2012), it does not appear plausible to assume that the reactants lie only in their respective ground states in chemical reactions. However, up to now, the effect of the higher-lying states on chemical reactions (for the system) has not been studied in the context of chemical modeling. It is this point that we study in the present paper.

Rate coefficients assuming that higher rotational levels can be populated – the so-called species-to-species rate coefficients – can be constructed from the state-to-state rate coefficients calculated by Hugo et al. (2009). Below, we show how this is accomplished in practice.

When the quantum states of the reactants and products are resolved, a chemical reaction can be written as where the states are labeled with i, j, m, and n. The states considered in Hugo et al. (2009) are the rotational levels of the isotopologs and H₂ in their ground vibrational states. The ground state-to-species coefficients pertain to the reaction where the species A and B are in their ground states and C and D can enter into any state upon formation. The coefficient k₀₀ is obtained through summation over the possible product states: The species-to-species rate coefficient, , is defined in terms of the total formation rates of C and D in reactions between A and B: where [A_i] is the number density (in units of cm^-3) of species A in state i, and so on, and k_ij is the sum over all the product states m and n. The populations of the energy levels of A and B depend on the gas density and temperature, the cross-sections for collisional transitions, and the Einstein coefficients of the radiative transitions between the energy levels. For example, if we have a reason, based on these parameters, to believe that only the two lowest energy levels of A and only the ground state of B are populated, the species-to-species rate coefficient can be calculated from At very high densities, or when the Einstein coefficients are very small, collisional excitation overpowers radiative transitions, and the level populations follow the Boltzmann distribution. In this case, the species-to-species rate coefficient can be written as where T is the kinetic temperature, and are the state energies (in K), and are the statistical weights of the levels (a product of the spin and rotational statistical weights), and Q^A(T) and Q^B(T) are the partition functions, The energy levels of the various rotational states, and the nuclear spin and rotational weights, can be read off Table II in Hugo et al. (2009).

Fig. 2 Total abundances (sums over spin states) of the various isotopologs as functions of time. The medium density is n(H2) = 106 cm-3 (upper row) or n(H2) = 107 cm-3 (lower row). From left to right, the panels show calculations assuming Tgas = Tdust = 10 , 15, or 20 K. Species-to-species rate coefficients are adopted in two of the models (method 1, dashed lines; method 2, solid lines). The dotted lines show the results of calculations using the ground-state-to-species rate coefficients.

The use of the species-to-species rate coefficients is only sensible if the medium density is high enough so that rotational states above the ground state can be assumed to be (significantly) populated, that is, if the medium density is comparable to or above the critical density of a given rotational transition. Setting the appropriate values of the critical densities is not straightforward because they depend on the temperature. Furthermore, the relative strengths of collisional and radiative excitation change gradually as the density increases. In this work, we define the critical density as a value of medium density for which a rotationally excited level population is 0.8 times the value expected from the Boltzmann distribution. This choice is arbitrary, and its effect is discussed in Sect. 4. The level populations were determined using the radiative transfer program of Juvela (1997). In this scheme the critical densities also depend on the abundances of the various species.

The calculation of the critical densities is discussed in greater detail in Appendix A, where we present the critical densities of H₂D⁺ and D₂H⁺ as functions of temperature and abundance. Table 1 shows the critical densities of the various rotational energy levels of H₂D⁺ and D₂H⁺ adopted in this paper. These values have been collected from the tables given in Appendix A and are shown here for convenience. Evidently, for a medium density of n(H₂) = 10⁶ cm^-3, for example, we can expect the lowest two rotational levels of oH₂D⁺ to be populated at low temperature², while pH₂D⁺, pD₂H⁺, and oD₂H⁺ should all lie mainly in their rotational ground states.

The final value of the species-to-species rate coefficient depends on which energy levels are taken into account. In this paper we consider two different approaches to choosing the included levels. We discuss these approaches next.

Fig. 3 Spin-state abundance ratios of the various isotopologs as functions of time. The medium density is n(H2) = 106 cm-3 (upper row) or n(H2) = 107 cm-3 (lower row). From left to right, the panels show calculations assuming Tgas = Tdust = 10 , 15, or 20 K. Species-to-species rate coefficients are adopted in two of the models (method 1, dashed lines; method 2, solid lines). The dotted lines show the results of calculations using the ground-state-to-species rate coefficients.

2.3.1. Local thermal equilibrium

We consider first a scheme where all of the excited rotational states are accessible as long as the medium density is higher than the critical density of the first excited rotational state of (o, p)H₂D⁺ or (o, p)D₂H⁺. This situation corresponds to local thermodynamic equilibrium (LTE). One great advantage of this approach is that it allows the construction of a reaction set that can be easily read into a chemical model. We calculated the species-to-species rate coefficients for all reactions included in the reacting system and fitted the results with a modified Arrhenius rate law in the temperature range 5–50 K. The resulting reaction set is given in Table B.2. We stress that the rate coefficients given in this table are only applicable for (o, p)H₂D⁺ or (o, p)D₂H⁺ depending on the medium density as explained above. Because and are homonuclear molecules and do not have a permanent dipole moment, we assume that the species-to-species rate coefficients can be used at all medium densities for these species. In what follows, we refer to this LTE-based approach as “method 1”.

2.3.2. Restricted states

A more careful treatment of the (o, p)H₂D⁺ or (o, p)D₂H⁺ species-to-species rate coefficients involves selecting only those states that have a critical density below the medium density. Further restrictions apply: for example if the medium density is n(H₂) = 5 × 10⁷ cm^-3, we include only the 1₀₁, 1₁₀, and 2₁₂ rotational levels of pD₂H⁺ and not the 3₀₃ level, even though it is allowed by the medium density, because the 2₂₁ level is not accessible owing to our assumptions (see Table 1). However, for and , we again assume that all levels can be populated regardless of the medium density.

Because the values of the rate coefficients are now strictly tied to the density, the rate coefficients need to be calculated on a case-by-case basis and the construction of a ready-made reaction set is not practical. Instead, the calculation of the rate coefficients is performed internally in our chemical code. We call this restricted-state approach “method 2”.

3. Results

3.1. Single-point models

Figure 2 shows the abundances (sums over spin states) of the isotopologs as calculated with single-point chemical models assuming different values of medium density and temperature (T_gas = T_dust). Figure 3 shows the spin-state abundance ratios in the same models. One feature of the models is immediately evident: the difference between methods 1 and 2 is small, that is, one can employ the species-to-species rate coefficients given in Table B.2 with good confidence when modeling cold and dense environments. We checked that at T = 50 K the difference between methods 1 and 2 remains smaller than a factor of two. However, at such a high temperature the abundances of the isotopologs are so low that the chosen method is of no practical significance.

The total abundances of the isotopologs are unaffected by the changes in the rate coefficients at T = 10 K. The difference between the ground-state-to-species model and the species-to-species model increases slightly with temperature, but more strongly with density, which is expected since the effect of the higher-lying rotational states becomes larger as consecutively higher states are populated. The general tendency in the species-to-species model is that the deuteration degree decreases with respect to the ground-state-to-species model toward higher temperatures and densities. This is due to the activation of several backward reaction channels that are suppressed when using the ground state-to-species coefficients. The simplistic approach of method 1 slightly underestimates the deuteration degree at n(H₂) = 10⁷ cm^-3 and T = 20 K. We discuss the rate coefficients of key reactions in Sect. 4.

The spin-state abundance ratios (Fig. 3) are naturally also affected by the changes in the rate coefficients. At n(H₂) = 10⁶ cm^-3 and T = 10 K, only the ortho/para ratio of H₂D⁺ is modified, which is expected based on the critical densities (Table 1) as oH₂D⁺ is the only species with accessible excited rotational states at n(H₂) = 10⁶ cm^-3. Notably, the ortho/para ratios of and are almost unchanged at these conditions even though the species-to-species rate coefficients are assumed to be applicable at all densities for these two species. This is because of the large energy differences between the first excited rotational states and the ground states for and (Hugo et al. 2009), which hinders rotational excitation at low temperature. The overall difference between the ground-state-to-species and species-to-species models increases with temperature and density, as was the case with the total abundances (Fig. 2). However, the difference between the two models for any given spin-state abundance ratio is typically less than a factor of two.

3.2. Source models: IRAS 16293 and L1544

Fig. 4 Left: radial distributions of the fractional abundances of selected species at t = 5 × 105 yr in a protostellar core resembling IRAS 16293 according to the model of Crimier et al. (2010). Solid lines correspond to the species-to-species model (method 1), while dashed lines correspond to the ground-state-to-species model. Middle: radial distributions of the o/p ratios of selected species, and the meta/ortho ratio of . The thin solid lines show the thermal spin-state ratios of the plotted species. Right: density and temperature distributions of the IRAS 16293 core model. The blue and red horizontal lines mark the critical densities of the first excited rotational transitions of oH2D+ and pD2H+, respectively.

It is not obvious based on the single-point models how and if the abundances and spin-state ratios of the isotopologs change in the context of source models with radially-varying density and temperature profiles, when one switches from ground-state-to-species rate coefficients to species-to-species rate coefficients. Potential problems arising from the use of ground-state-to-species rate coefficients for the system in the interpretation of observational data was discussed by Harju et al. (2017). Following the discussion in that paper, we carried out chemical calculations using a physical model for the protostellar system IRAS 16293-2422 A/B (Crimier et al. 2010; see also Brünken et al. 2014; Harju et al. 2017). The chemical modeling setup is essentially the same as described in Brünken et al. (2014), where the physical model is also described in detail. The source model, which assumes spherical symmetry, is separated into concentric shells and the chemical evolution in each shell is tracked using our pseudo-time-dependent chemical code (Sect. 2.1). The model outputs the abundances of the various species as functions of radius and time, allowing us to reconstruct the time-evolution of the isotopologs (among others) in the model core. We carried out chemical calculations using the ground-state-to-species and species-to-species rate coefficients for the system and compared the difference in predicted radial abundance profiles for the isotopologs and specifically the o/p ratios of H₂D⁺ and D₂H⁺.

Fig. 5 Left: radial distributions of the fractional abundances of the spin states of H2D+ and D2H+ at t = 1 × 106 yr in the innermost 10 000 AU of the L1544 model. Solid lines correspond to the species-to-species model (method 1), while dashed lines correspond to the ground-state-to-species model. Middle: radial distributions of the o/p ratios of H2D+ and D2H+. Linestyles are the same as in the left panel. Right: density and temperature distributions of the L1544 core model. The blue and red horizontal lines mark the critical densities of the first excited rotational transitions of oH2D+ and pD2H+, respectively.

Fig. 6 Rate coefficients of the pH2D+ + oH2− → oH2D+ + pH2 (left), oH2D+ + pH2− → pH2D+ + oH2 (middle), and (right) reactions as functions of temperature. Red lines represent ground-state-to-species rate coefficients, while blue lines represent species-to-species rate coefficients. Solid lines represent the raw data from Hugo et al. (2009), and dashed lines represent our fits using the modified Arrhenius rate law (see text). We highlight the different y-axis scaling in the panels.

Figure 4 shows the total abundances (summed over spin states) of the isotopologs and the o/p ratios of H₂D⁺, D₂H⁺, and H₂, and the meta/ortho ratio of , at t = 5.0 × 10⁵ yr in the IRAS 16293 model. The figure concentrates on the radius range where most of the H₂D⁺ and D₂H⁺ emission/absorption originates (Brünken et al. 2014; Harju et al. 2017). The total abundances are not significantly affected by the choice of the rate coefficients, but the spin-state abundance ratios show differences depending on the rate coefficients used. In particular, we recover exactly the type of behavior expected on the basis of the discussion in Harju et al. (2017): The H₂D⁺ o/p ratio decreases by a factor of ~1.5 and the D₂H⁺ o/p ratio increases by a similar factor when one switches from the ground-state-to-species to species-to-species rate coefficients. The sharp features in the o/p ratios evident at R ~ 2500 AU for H₂D⁺ and at R ~ 1000 AU for D₂H⁺ are a result of the activation of the species-to-species rate coefficients at the densities corresponding to these radii (Table 1). We note that for H₂D⁺ for example, the solid and dashed lines do not overlap below the critical density of the first excited rotational state, even though the ground-state-to-species rate coefficients apply in this regime. This is because of the effect of and , for which the species-to-species rate coefficients apply in all conditions.

To investigate the effect of the species-to-species rates in prestellar core conditions, we performed another test using the physical model for the prestellar core L1544 discussed in Keto & Caselli (2010) and Keto et al. (2014). The chemical calculation procedure was identical to that used in the IRAS 16293 modeling. Figure 5 shows the results of the calculations at t = 10⁶ yr, zoomed in to the innermost 10 000 AU of the core model. The density is high enough only in the central 2000 AU to activate the species-to-species rate coefficients for oH₂D⁺, while pD₂H⁺ is only affected in the central ~500 AU. The results from the two types of model are however more or less identical, which is caused by the low temperature in the central regions. This test confirms the tendencies shown in Fig. 2, that is, that the species-to-species rate coefficients become important at high density only if the temperature is higher than 10 K.

4. Discussion

The results presented in Sect. 3 clearly show that the switch from ground-state-to-species rate coefficients to species-to-species rate coefficients can have a marked impact on the chemistry of the isotopologs. As discussed in Brünken et al. (2014) and Harju et al. (2017), the pH₂D⁺ + oH₂ ↔ oH₂D⁺ + pH₂ reaction, endothermic in the backward direction, plays a major role in the evolution of the H₂D⁺ o/p ratio. Figure 6 shows the rate coefficient of this reaction in the forward and backward directions. The ground-state-to-species model slightly overestimates the destruction of pH₂D⁺ at T> 10 K, although the difference between the two cases is only ~20%. The backward rate coefficient is however clearly higher in the species-to-species model (calculated with method 1) throughout the temperature range considered; the difference is about a factor of 1.8 at T = 10 K, and ~1.5–2 throughout the temperature range considered. This difference in the backward rate coefficient – although only of a factor of two – can translate to a decrease in the H₂D⁺ o/p ratio as evident in the figures presented in Sect. 3, and can also be observable (Harju et al. 2017).

For many reactions pertaining to the system, the ground-state-to-species and species-to-species rate coefficients are very similar, showing maximum differences on the order of 10%. The general tendency is that large differences are usually seen in the rates of endothermic reactions, which are underestimated by the ground-state-to-species model, as displayed in Fig. 6. The reaction (shown in the right panel of Fig. 6) was included to demonstrate that sometimes the difference between the rate coefficients may be larger than one order of magnitude (for this reaction about a factor of 12 at 20 K). It is clear from the above that one should employ the species-to-species rate coefficients in a chemical model to obtain a better estimate of the abundances of the spin states of the isotopologs.

The modeling estimates for the abundances of the isotopologs are tied to the critical densities, as this affects when the species-to-species rate coefficients are assumed to apply. The critical densities depend relatively strongly on the abundances of the species at high density where optical depth effects become important. However, it is indeed the abundances that we wish to obtain from the chemical modeling. To derive truly consistent values for the various critical densities, one should iterate the chemical calculations in order to obtain successively better estimates for the critical densities and, in turn, the abundances. In this paper we settled on a simple scheme where the critical densities (Table 1) were chosen on the basis of steady-state abundances obtained from our rate coefficient fitting test (Fig. 1).

The choice of the cutoff for when the species-to-species rates are assumed apply (see Appendix A for more details) is arbitrary, and it is understood that the use of species-to-species rate coefficients adds to the parameter space in the modeling, and so is an additional source of uncertainty. Altering the critical densities would affect the radii where the discontinuities in Figs. 4 and 5 appear. Here we considered that a rotationally excited level is taken into account in the calculation of the rate coefficient if its population is ≳0.8 times that of the thermal population. The goodness of this choice should be tested through a comparison against a state-to-state chemistry model, which is, however, beyond the scope of the present paper. The species most affected by the choice of the threshold value is oH₂D⁺. Using a threshold ratio of 0.9 would shift the discontinuity from 2500 AU to 1500 AU, whereas with a threshold of 0.7, the discontinuity would appear at a distance of 3500 AU from the center. The threshold 0.9 gives practically the same abundance profile as the ground-state-to-species model does. The mentioned changes alter the line-of-sight average abundance of oH₂D⁺ by less than 10%. Combined with the steep density and temperature gradients of this model, the small changes in the oH₂D⁺ abundance profile have, however, a noticeable effect on the observed line intensity. We note that the uncertainty associated with the threshold is not likely to be greater than, for example, the uncertainty in the temperature in the modeled objects (typically ~1–2 K for the cold gas). To summarize the experiences from the present study, an accurate model of the chemistry of the isotopologs should include the effect of excited rotational states, but the outcome also depends on how well the density and temperature distributions of the object are known.

5. Conclusions

We studied the chemistry of the spin states of the isotopologs using a model where rotational excitation of the reactants is taken into account (the species-to-species model), as opposed to previous studies of the subject which have thus far assumed that only the rotational ground state is populated (the ground-state-to-species model). We considered two different methods for constructing the species-to-species rate coefficients: the first method assumes that all rotational states can be populated once the medium density is higher than the critical density of the first excited rotational state, while in the second method only those rotational states whose critical densities lie below the medium density are included.

We found that the two methods for constructing the species-to-species rate coefficients produce similar results in the physical conditions considered here, which means that the species-to-species rate coefficients can be conveniently read from a precalculated table with good confidence, and are thus easily implementable in a chemical model. On the other hand, the abundances and spin-state ratios predicted by the ground-state-to-species and species-to-species models are different from each other at high density, depending also on the temperature. Notable differences between the models are seen for T> 10 K where the ground-state-to-species model overestimates the abundances of the isotopologs; by a factor of ~2 for H₂D⁺, and even up to an order of magnitude for . The species-to-species model is the most realistic of the two. We note however that the species-to-species rate coefficients introduce an additional uncertainty into the modeling through the critical densities of the rotational transitions of H₂D⁺ and D₂H⁺, which determine when the species-to-species model is applicable.

The new model is a step toward full state-to-state modeling of the isotopolog chemistry. A state-to-state model would not suffer from the problem of setting the appropriate values of critical densities in the way that the present model does, but the implementation of state-to-state rates into a complete gas-grain chemical model including elements heavier than hydrogen is a monumental task. Such an effort, although certainly called for, is left for future work. In the meantime, the species-to-species rate coefficients represent the next best thing, and should be implemented in models of the chemistry of the reacting system.

Acknowledgments

P.C. acknowledges financial support of the European Research Council (ERC; project PALs 320620).

References

Appendix A: Determining the critical densities

In the present study, the densities where the rotational levels of H₂D⁺ and D₂H⁺ become significantly populated were estimated through radiative transfer calculations. The core model had a steep density gradient, increasing from n(H₂) = 10⁵ cm^-3 at the outer boundary to 3 × 10⁹ cm^-3 in the center. The core was assumed to be isothermal, and we varied the temperature in the range 5–25 K. Different fractional abundances were also tested to examine optical thickness effects. The Monte Carlo program of Juvela (1997) was used in the simulations. The results for the temperature T = 15 K are shown in Fig. A.1. The curves show the ratios of the populations to the values expected from the Boltzmann distribution. In the figure, the fractional abundances of all the four species, oH₂D⁺, pH₂D⁺, oD₂H⁺, and pD₂H⁺, are assumed to be X = 10^-10. A level is thermalized when the ratio is unity. In “method 2” of the chemistry model described in Sect. 2.3.2, a level is taken into account, and its population relative to the ground level is calculated by the Boltzmann factor, at densities where the ratio to the thermal value exceeds 0.8. This limit is shown with dashed horizontal lines in Fig. A.1. The critical densities of oH₂D⁺, pH₂D⁺, oD₂H⁺, and pD₂H⁺ calculated with the radiative transfer model for different temperatures and for different abundances of the species are collected in Tables A.1 to A.3, from which the values used in our model calculations, given in Table 1, were compiled.

For completeness we show in Fig. A.2 similar results as in Fig. A.1, but assuming X = 10^-8 for H₂D⁺ and D₂H⁺, which can be considered as a strong upper limit as the abundances are in reality unlikely to rise this high. The increase in the fractional abundance increases the optical thickness of the rotational transitions. The effect is that photons are absorbed and re-emitted several times before escaping the cloud, and the role of collisional transitions in determining the level populations becomes more important. This so-called “radiative trapping” influences the level populations in the same way as an increase in the gas density, and results in a greater thermalisation of the populations (Walmsley 1987). Comparison between Figs. A.1 and A.2 shows that assuming fractional abundances of X = 10^-8 causes the thermalization of the lowest rotational levels to occur at clearly lower densities than in the case of X = 10^-10. The levels 2₁₂ and 2₁₁ of oH₂D⁺ become mildly supra-thermally excited (the maximum excitation temperatures are T_ex ~ 15.4–15.5 K) when X = 10^-8 is assumed. This is caused by rapid radiative decay from the collisionally excited higher rotational levels. We note that for canonical abundances (X = 10^-11−10^-9) supra-thermal populations are never found for rotationally excited levels of H₂D⁺ and D₂H⁺.

Fig. A.1 Excitation of the lowest rotational levels of oH2D+, pH2D+, oD2H+, and pD2H+ as functions of the gas density at T = 15 K, assuming a fractional abundance X = 10-10 for each species. The curves show the ratios of the populations to the values expected from the Boltzmann distribution. The dashed horizontal lines indicates the limit where a level is considered to be significantly populated and taken into account in the calculation of the species-to-species rate coefficients.

Fig. A.2 As Fig. A.1, but assuming X = 10-8 for each species.

Appendix B: Ground state-to-species and species-to-species rate coefficients

All Tables

All Figures

Fig. 1 Abundances (upper row) and abundance ratios (lower row) of the various isotopologs as functions of time. The left-hand panels correspond to n(H2) = 105 cm-3, while the right-hand panels correspond to n(H2) = 106 cm-3. The temperature is set to Tgas = Tdust = 10 K. Solid lines represent calculations using the Hugo et al. (2009) ground-state-to-species rate coefficients, while dashed lines represent calculations using our new fit to the same coefficients (see text).

In the text

Fig. 2 Total abundances (sums over spin states) of the various isotopologs as functions of time. The medium density is n(H2) = 106 cm-3 (upper row) or n(H2) = 107 cm-3 (lower row). From left to right, the panels show calculations assuming Tgas = Tdust = 10 , 15, or 20 K. Species-to-species rate coefficients are adopted in two of the models (method 1, dashed lines; method 2, solid lines). The dotted lines show the results of calculations using the ground-state-to-species rate coefficients.

In the text

Fig. 3 Spin-state abundance ratios of the various isotopologs as functions of time. The medium density is n(H2) = 106 cm-3 (upper row) or n(H2) = 107 cm-3 (lower row). From left to right, the panels show calculations assuming Tgas = Tdust = 10 , 15, or 20 K. Species-to-species rate coefficients are adopted in two of the models (method 1, dashed lines; method 2, solid lines). The dotted lines show the results of calculations using the ground-state-to-species rate coefficients.

In the text

Fig. 4 Left: radial distributions of the fractional abundances of selected species at t = 5 × 105 yr in a protostellar core resembling IRAS 16293 according to the model of Crimier et al. (2010). Solid lines correspond to the species-to-species model (method 1), while dashed lines correspond to the ground-state-to-species model. Middle: radial distributions of the o/p ratios of selected species, and the meta/ortho ratio of . The thin solid lines show the thermal spin-state ratios of the plotted species. Right: density and temperature distributions of the IRAS 16293 core model. The blue and red horizontal lines mark the critical densities of the first excited rotational transitions of oH2D+ and pD2H+, respectively.

In the text

Fig. 5 Left: radial distributions of the fractional abundances of the spin states of H2D+ and D2H+ at t = 1 × 106 yr in the innermost 10 000 AU of the L1544 model. Solid lines correspond to the species-to-species model (method 1), while dashed lines correspond to the ground-state-to-species model. Middle: radial distributions of the o/p ratios of H2D+ and D2H+. Linestyles are the same as in the left panel. Right: density and temperature distributions of the L1544 core model. The blue and red horizontal lines mark the critical densities of the first excited rotational transitions of oH2D+ and pD2H+, respectively.

In the text

Fig. 6 Rate coefficients of the pH2D+ + oH2− → oH2D+ + pH2 (left), oH2D+ + pH2− → pH2D+ + oH2 (middle), and (right) reactions as functions of temperature. Red lines represent ground-state-to-species rate coefficients, while blue lines represent species-to-species rate coefficients. Solid lines represent the raw data from Hugo et al. (2009), and dashed lines represent our fits using the modified Arrhenius rate law (see text). We highlight the different y-axis scaling in the panels.

In the text

Fig. A.1 Excitation of the lowest rotational levels of oH2D+, pH2D+, oD2H+, and pD2H+ as functions of the gas density at T = 15 K, assuming a fractional abundance X = 10-10 for each species. The curves show the ratios of the populations to the values expected from the Boltzmann distribution. The dashed horizontal lines indicates the limit where a level is considered to be significantly populated and taken into account in the calculation of the species-to-species rate coefficients.

In the text

	Fig. A.2 As Fig. A.1, but assuming X = 10-8 for each species.
In the text