Issue 
A&A
Volume 607, November 2017



Article Number  A26  
Number of page(s)  21  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201731039  
Published online  31 October 2017 
Speciestospecies rate coefficients for the H_{3}^{+} + H_{2} reacting system
^{1} MaxPlanckInstitute for Extraterrestrial Physics (MPE), Giessenbachstr. 1, 85748 Garching, Germany
email: osipila@mpe.mpg.de
^{2} Department of Physics, PO Box 64, 00014 University of Helsinki, Finland
Received: 25 April 2017
Accepted: 10 July 2017
Aims. We study whether or not rotational excitation can make a large difference to chemical models of the abundances of the H_{3}^{+} isotopologs, including spin states, in physical conditions corresponding to starless cores and protostellar envelopes.
Methods. We developed a new rate coefficient set for the chemistry of the H_{3}^{+} isotopologs, allowing for rotational excitation, using previously published statetostate rate coefficients. These new socalled speciestospecies rate coefficients are compared with previouslyused groundstatetospecies rate coefficients by calculating chemical evolution in variable physical conditions using a pseudotimedependent chemical code.
Results. We find that the new speciestospecies model produces different results to the ground statetospecies model at high density and toward increasing temperatures (T> 10 K). The most prominent difference is that the speciestospecies model predicts a lower H_{3}^{+} deuteration degree at high density owing to an increase of the rate coefficients of endothermic reactions that tend to decrease deuteration. For example at 20 K, the groundstatetospecies model overestimates the abundance of H_{2}D^{+} by a factor of about two, while the abundance of D_{3}^{+} can differ by up to an order of magnitude between the models. The spinstate abundance ratios of the various H_{3}^{+} isotopologs are also affected, and the new model better reproduces recent observations of the abundances of ortho and para H_{2}D^{+} and D_{2}H^{+}. The main caveat is that the applicability regime of the new rate coefficients depends on the critical densities of the various rotational transitions which vary with the abundances of the species and the temperature in dense clouds.
Conclusions. The difference in the abundances of the H_{3}^{+} isotopologs predicted by the speciestospecies and ground statetospecies models is negligible at 10 K corresponding to physical conditions in starless cores, but inclusion of the excited states is very important in studies of deuteration at higher temperatures, for example in protostellar envelopes. The speciestospecies rate coefficients provide a more realistic approach to the chemistry of the H_{3}^{+} isotopologs than the groundstatetospecies rate coefficients do, and so the former should be adopted in chemical models describing the chemistry of the H_{3}^{+} + H_{2} reacting system.
Key words: astrochemistry / ISM: clouds / ISM: molecules / ISM: abundances
© ESO, 2017
1. Introduction
The ion is a universal proton donor and has therefore a central role in ionmolecule 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 orthopara conversion of H_{2} 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_{2} compared to that of pH_{2} (~ 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 highdensity and lowtemperature regions of prestellar cores, where species heavier than He may be highly depleted because of freezeout (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_{2}D^{+}) and the linked deuteration fraction measured in, for example, N_{2}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 statetostate 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_{2}, is of fundamental importance for the physics and chemistry of interstellar clouds.
Hugo et al. (2009) derived statetostate thermal rate coefficients for inelastic and reactive collisions between all the isotopic variants of and H_{2} and their different spin modifications. They compiled a table of groundstatetospecies 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 statetostate rate coefficients can also be used to estimate speciestospecies 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_{2}) ~ 10^{5} cm^{3}, but in prestellar and starforming cores with central densities exceeding 10^{6} cm^{3}, some of the lowest rotationally excited levels of , H_{2}D^{+}, D_{2}H^{+}, and should be excited and contribute to the total speciestospecies reaction rates. Besides the construction of rate coefficients for chemical reactions, the statetostate coefficients given by Hugo et al. (2009) are needed for calculating the populations of the rotationally excited levels of H_{2}D^{+} and other isotopologs in connection with radiative transfer calculations.
A discrepancy between the predictions of chemical models utilizing either the groundstatetospecies or speciestospecies rate coefficients was recently suggested to arise in the protostellar core IRAS 162932422, where the model using the groundstatetospecies rate coefficients seems to overpredict the oH_{2}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 speciestospecies rate coefficients from the statetostate 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_{2}D^{+}, D_{2}H^{+}, and at different densities. For H_{2}D^{+} and D_{2}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 speciestospecies rate coefficients or the previously used groundstatetospecies rate coefficients. We employ the gasgrain 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 gasphase network (Wakelam et al. 2015) as the basis upon which the deuterium and spinstate chemistry is added according to the procedures discussed in detail in Sipilä et al. (2015a,b).
2.2. New fits to the ground statetospecies rate coefficients
Hugo et al. (2009) calculated the statetostate 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 socalled groundstatetospecies rate coefficients were produced by fitting a twoparameter curve of the form αexp(−γ/T) in the temperature range 5–20 K^{1} to data obtained from the statetostate 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 speciestospecies rate coefficients (see Sect. 2.3) pertain to the temperature range 5–50 K, for consistency we made a new fit to the groundstatetospecies 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. statesostate data). However, this reaction is insignificant as its rate coefficient is of the order of 10^{18} cm^{3} s^{1} at 10 K.
Figure 1 shows a comparison of the abundances of the spin states of the isotopologs calculated with the groundstatetospecies 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 spinstate abundance ratios present slight deviations for the doubly and triply deuterated species (the variations for and H_{2}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 lefthand panels correspond to n(H_{2}) = 10^{5} cm^{3}, while the righthand panels correspond to n(H_{2}) = 10^{6} cm^{3}. The temperature is set to T_{gas} = T_{dust} = 10 K. Solid lines represent calculations using the Hugo et al. (2009) groundstatetospecies rate coefficients, while dashed lines represent calculations using our new fit to the same coefficients (see text). 
Critical densities (n_{c}) of the excited rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009), used in the models presented in this paper.
2.3. Speciestospecies rate coefficients
Given that H_{2}D^{+} emission is observed to be ubiquitous (lowmass starless, prestellar, and protostellar cores, Caselli et al. 2008; highmass 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 higherlying 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 socalled speciestospecies rate coefficients – can be constructed from the statetostate 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_{2} in their ground vibrational states. The ground statetospecies 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_{00} is obtained through summation over the possible product states: The speciestospecies 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 crosssections 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 speciestospecies 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 speciestospecies 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(H_{2}) = 10^{6} cm^{3} (upper row) or n(H_{2}) = 10^{7} cm^{3} (lower row). From left to right, the panels show calculations assuming T_{gas} = T_{dust} = 10 , 15, or 20 K. Speciestospecies 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 groundstatetospecies rate coefficients. 
The use of the speciestospecies 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_{2}D^{+} and D_{2}H^{+} as functions of temperature and abundance. Table 1 shows the critical densities of the various rotational energy levels of H_{2}D^{+} and D_{2}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_{2}) = 10^{6} cm^{3}, for example, we can expect the lowest two rotational levels of oH_{2}D^{+} to be populated at low temperature^{2}, while pH_{2}D^{+}, pD_{2}H^{+}, and oD_{2}H^{+} should all lie mainly in their rotational ground states.
The final value of the speciestospecies 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 Spinstate abundance ratios of the various isotopologs as functions of time. The medium density is n(H_{2}) = 10^{6} cm^{3} (upper row) or n(H_{2}) = 10^{7} cm^{3} (lower row). From left to right, the panels show calculations assuming T_{gas} = T_{dust} = 10 , 15, or 20 K. Speciestospecies 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 groundstatetospecies 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_{2}D^{+} or (o, p)D_{2}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 speciestospecies 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_{2}D^{+} or (o, p)D_{2}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 speciestospecies rate coefficients can be used at all medium densities for these species. In what follows, we refer to this LTEbased approach as “method 1”.
2.3.2. Restricted states
A more careful treatment of the (o, p)H_{2}D^{+} or (o, p)D_{2}H^{+} speciestospecies 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_{2}) = 5 × 10^{7} cm^{3}, we include only the 1_{01}, 1_{10}, and 2_{12} rotational levels of pD_{2}H^{+} and not the 3_{03} level, even though it is allowed by the medium density, because the 2_{21} 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 casebycase basis and the construction of a readymade reaction set is not practical. Instead, the calculation of the rate coefficients is performed internally in our chemical code. We call this restrictedstate approach “method 2”.
3. Results
3.1. Singlepoint models
Figure 2 shows the abundances (sums over spin states) of the isotopologs as calculated with singlepoint chemical models assuming different values of medium density and temperature (T_{gas} = T_{dust}). Figure 3 shows the spinstate 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 speciestospecies 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 groundstatetospecies model and the speciestospecies model increases slightly with temperature, but more strongly with density, which is expected since the effect of the higherlying rotational states becomes larger as consecutively higher states are populated. The general tendency in the speciestospecies model is that the deuteration degree decreases with respect to the groundstatetospecies model toward higher temperatures and densities. This is due to the activation of several backward reaction channels that are suppressed when using the ground statetospecies coefficients. The simplistic approach of method 1 slightly underestimates the deuteration degree at n(H_{2}) = 10^{7} cm^{3} and T = 20 K. We discuss the rate coefficients of key reactions in Sect. 4.
The spinstate abundance ratios (Fig. 3) are naturally also affected by the changes in the rate coefficients. At n(H_{2}) = 10^{6} cm^{3} and T = 10 K, only the ortho/para ratio of H_{2}D^{+} is modified, which is expected based on the critical densities (Table 1) as oH_{2}D^{+} is the only species with accessible excited rotational states at n(H_{2}) = 10^{6} cm^{3}. Notably, the ortho/para ratios of and are almost unchanged at these conditions even though the speciestospecies 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 groundstatetospecies and speciestospecies 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 spinstate 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 × 10^{5} yr in a protostellar core resembling IRAS 16293 according to the model of Crimier et al. (2010). Solid lines correspond to the speciestospecies model (method 1), while dashed lines correspond to the groundstatetospecies 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 spinstate 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 oH_{2}D^{+} and pD_{2}H^{+}, respectively. 
It is not obvious based on the singlepoint models how and if the abundances and spinstate ratios of the isotopologs change in the context of source models with radiallyvarying density and temperature profiles, when one switches from groundstatetospecies rate coefficients to speciestospecies rate coefficients. Potential problems arising from the use of groundstatetospecies 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 162932422 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 pseudotimedependent 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 timeevolution of the isotopologs (among others) in the model core. We carried out chemical calculations using the groundstatetospecies and speciestospecies 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_{2}D^{+} and D_{2}H^{+}.
Fig. 5 Left: radial distributions of the fractional abundances of the spin states of H_{2}D^{+} and D_{2}H^{+} at t = 1 × 10^{6} yr in the innermost 10 000 AU of the L1544 model. Solid lines correspond to the speciestospecies model (method 1), while dashed lines correspond to the groundstatetospecies model. Middle: radial distributions of the o/p ratios of H_{2}D^{+} and D_{2}H^{+}. 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 oH_{2}D^{+} and pD_{2}H^{+}, respectively. 
Fig. 6 Rate coefficients of the pH_{2}D^{+} + oH_{2}− → oH_{2}D^{+} + pH_{2} (left), oH_{2}D^{+} + pH_{2}− → pH_{2}D^{+} + oH_{2} (middle), and (right) reactions as functions of temperature. Red lines represent groundstatetospecies rate coefficients, while blue lines represent speciestospecies 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 yaxis scaling in the panels. 
Figure 4 shows the total abundances (summed over spin states) of the isotopologs and the o/p ratios of H_{2}D^{+}, D_{2}H^{+}, and H_{2}, and the meta/ortho ratio of , at t = 5.0 × 10^{5} yr in the IRAS 16293 model. The figure concentrates on the radius range where most of the H_{2}D^{+} and D_{2}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 spinstate 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_{2}D^{+} o/p ratio decreases by a factor of ~1.5 and the D_{2}H^{+} o/p ratio increases by a similar factor when one switches from the groundstatetospecies to speciestospecies rate coefficients. The sharp features in the o/p ratios evident at R ~ 2500 AU for H_{2}D^{+} and at R ~ 1000 AU for D_{2}H^{+} are a result of the activation of the speciestospecies rate coefficients at the densities corresponding to these radii (Table 1). We note that for H_{2}D^{+} for example, the solid and dashed lines do not overlap below the critical density of the first excited rotational state, even though the groundstatetospecies rate coefficients apply in this regime. This is because of the effect of and , for which the speciestospecies rate coefficients apply in all conditions.
To investigate the effect of the speciestospecies 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^{6} 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 speciestospecies rate coefficients for oH_{2}D^{+}, while pD_{2}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 speciestospecies 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 groundstatetospecies rate coefficients to speciestospecies 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_{2}D^{+} + oH_{2} ↔ oH_{2}D^{+} + pH_{2} reaction, endothermic in the backward direction, plays a major role in the evolution of the H_{2}D^{+} o/p ratio. Figure 6 shows the rate coefficient of this reaction in the forward and backward directions. The groundstatetospecies model slightly overestimates the destruction of pH_{2}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 speciestospecies 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_{2}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 groundstatetospecies and speciestospecies 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 groundstatetospecies 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 speciestospecies 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 speciestospecies 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 steadystate abundances obtained from our rate coefficient fitting test (Fig. 1).
The choice of the cutoff for when the speciestospecies rates are assumed apply (see Appendix A for more details) is arbitrary, and it is understood that the use of speciestospecies 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 statetostate 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_{2}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 groundstatetospecies model does. The mentioned changes alter the lineofsight average abundance of oH_{2}D^{+} by less than 10%. Combined with the steep density and temperature gradients of this model, the small changes in the oH_{2}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 speciestospecies model), as opposed to previous studies of the subject which have thus far assumed that only the rotational ground state is populated (the groundstatetospecies model). We considered two different methods for constructing the speciestospecies 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 speciestospecies rate coefficients produce similar results in the physical conditions considered here, which means that the speciestospecies 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 spinstate ratios predicted by the groundstatetospecies and speciestospecies 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 groundstatetospecies model overestimates the abundances of the isotopologs; by a factor of ~2 for H_{2}D^{+}, and even up to an order of magnitude for . The speciestospecies model is the most realistic of the two. We note however that the speciestospecies rate coefficients introduce an additional uncertainty into the modeling through the critical densities of the rotational transitions of H_{2}D^{+} and D_{2}H^{+}, which determine when the speciestospecies model is applicable.
The new model is a step toward full statetostate modeling of the isotopolog chemistry. A statetostate 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 statetostate rates into a complete gasgrain 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 speciestospecies rate coefficients represent the next best thing, and should be implemented in models of the chemistry of the reacting system.
This is reinforced by the fact that the separation between the two lowest levels is only ~20 K (Hugo et al. 2009).
Acknowledgments
P.C. acknowledges financial support of the European Research Council (ERC; project PALs 320620).
References
 Albertsson, T., Indriolo, N., Kreckel, H., et al. 2014, ApJ, 787, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Brünken, S., Sipilä, O., Chambers, E. T., et al. 2014, Nature, 516, 219 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crimier, N., Ceccarelli, C., Maret, S., et al. 2010, A&A, 519, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dalgarno, A., & Lepp, S. 1984, ApJ, 287, L47 [NASA ADS] [CrossRef] [Google Scholar]
 Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2006, A&A, 449, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Friesen, R. K., Di Francesco, J., Bourke, T. L., et al. 2014, ApJ, 797, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Furuya, K., Aikawa, Y., Hincelin, U., et al. 2015, A&A, 584, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harju, J., Sipilä, O., Brünken, S., et al. 2017, ApJ, 840, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Hugo, E., Asvany, O., & Schlemmer, S. 2009, J. Chem. Phys., 130, 164302 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Juvela, M. 1997, A&A, 322, 943 [NASA ADS] [Google Scholar]
 Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625 [NASA ADS] [CrossRef] [Google Scholar]
 Keto, E., Rawlings, J., & Caselli, P. 2014, MNRAS, 440, 2616 [NASA ADS] [CrossRef] [Google Scholar]
 Kong, S., Caselli, P., Tan, J. C., Wakelam, V., & Sipilä, O. 2015, ApJ, 804, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Le Gal, R., HilyBlant, P., Faure, A., et al. 2014, A&A, 562, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Majumdar, L., Gratier, P., Ruaud, M., et al. 2017, MNRAS, 466, 4470 [Google Scholar]
 Pagani, L., Vastel, C., Hugo, E., et al. 2009, A&A, 494, 623 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pagani, L., Roueff, E., & Lesaffre, P. 2011, ApJ, 739, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Pagani, L., Lesaffre, P., Jorfi, M., et al. 2013, A&A, 551, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pillai, T., Caselli, P., Kauffmann, J., et al. 2012, ApJ, 751, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Sipilä, O., Hugo, E., Harju, J., et al. 2010, A&A, 509, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sipilä, O., Caselli, P., & Harju, J. 2015a, A&A, 578, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sipilä, O., Harju, J., Caselli, P., & Schlemmer, S. 2015b, A&A, 581, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teslja, A., & Valentini, J. J. 2006, J. Chem. Phys., 125, 132304 [NASA ADS] [CrossRef] [Google Scholar]
 Wakelam, V., Loison, J.C., Herbst, E., et al. 2015, ApJS, 217, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Walmsley, C. M. 1987, in Physical Processes in Interstellar Clouds, eds. G. E. Morfill, & M. Scholer, NATO ASI Ser., 210, 161 [Google Scholar]
 Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhang, D. H., & Guo, H. 2016, Ann. Rev. Phys. Chem., 67, 135 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Determining the critical densities
In the present study, the densities where the rotational levels of H_{2}D^{+} and D_{2}H^{+} become significantly populated were estimated through radiative transfer calculations. The core model had a steep density gradient, increasing from n(H_{2}) = 10^{5} cm^{3} at the outer boundary to 3 × 10^{9} 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_{2}D^{+}, pH_{2}D^{+}, oD_{2}H^{+}, and pD_{2}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_{2}D^{+}, pH_{2}D^{+}, oD_{2}H^{+}, and pD_{2}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_{2}D^{+} and D_{2}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 reemitted several times before escaping the cloud, and the role of collisional transitions in determining the level populations becomes more important. This socalled “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_{12} and 2_{11} of oH_{2}D^{+} become mildly suprathermally 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}) suprathermal populations are never found for rotationally excited levels of H_{2}D^{+} and D_{2}H^{+}.
Fig. A.1 Excitation of the lowest rotational levels of oH_{2}D^{+}, pH_{2}D^{+}, oD_{2}H^{+}, and pD_{2}H^{+} 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 speciestospecies rate coefficients. 
Critical densities (n_{c}) of the rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009) assuming T = 10 K.
Critical densities (n_{c}) of the rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009) assuming T = 15 K.
Critical densities (n_{c}) of the rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009) assuming T = 20 K.
Appendix B: Ground statetospecies and speciestospecies rate coefficients
Ground statetospecies rate coefficients for the system calculated from the statetostate data by Hugo et al. (2009).
Speciestospecies rate coefficients for the system calculated from the statetostate data by Hugo et al. (2009).
All Tables
Critical densities (n_{c}) of the excited rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009), used in the models presented in this paper.
Critical densities (n_{c}) of the rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009) assuming T = 10 K.
Critical densities (n_{c}) of the rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009) assuming T = 15 K.
Critical densities (n_{c}) of the rotational levels of H_{2}D^{+} and D_{2}H^{+}, in order of increasing energy (Hugo et al. 2009) assuming T = 20 K.
Ground statetospecies rate coefficients for the system calculated from the statetostate data by Hugo et al. (2009).
Speciestospecies rate coefficients for the system calculated from the statetostate data by Hugo et al. (2009).
All Figures
Fig. 1 Abundances (upper row) and abundance ratios (lower row) of the various isotopologs as functions of time. The lefthand panels correspond to n(H_{2}) = 10^{5} cm^{3}, while the righthand panels correspond to n(H_{2}) = 10^{6} cm^{3}. The temperature is set to T_{gas} = T_{dust} = 10 K. Solid lines represent calculations using the Hugo et al. (2009) groundstatetospecies 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(H_{2}) = 10^{6} cm^{3} (upper row) or n(H_{2}) = 10^{7} cm^{3} (lower row). From left to right, the panels show calculations assuming T_{gas} = T_{dust} = 10 , 15, or 20 K. Speciestospecies 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 groundstatetospecies rate coefficients. 

In the text 
Fig. 3 Spinstate abundance ratios of the various isotopologs as functions of time. The medium density is n(H_{2}) = 10^{6} cm^{3} (upper row) or n(H_{2}) = 10^{7} cm^{3} (lower row). From left to right, the panels show calculations assuming T_{gas} = T_{dust} = 10 , 15, or 20 K. Speciestospecies 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 groundstatetospecies rate coefficients. 

In the text 
Fig. 4 Left: radial distributions of the fractional abundances of selected species at t = 5 × 10^{5} yr in a protostellar core resembling IRAS 16293 according to the model of Crimier et al. (2010). Solid lines correspond to the speciestospecies model (method 1), while dashed lines correspond to the groundstatetospecies 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 spinstate 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 oH_{2}D^{+} and pD_{2}H^{+}, respectively. 

In the text 
Fig. 5 Left: radial distributions of the fractional abundances of the spin states of H_{2}D^{+} and D_{2}H^{+} at t = 1 × 10^{6} yr in the innermost 10 000 AU of the L1544 model. Solid lines correspond to the speciestospecies model (method 1), while dashed lines correspond to the groundstatetospecies model. Middle: radial distributions of the o/p ratios of H_{2}D^{+} and D_{2}H^{+}. 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 oH_{2}D^{+} and pD_{2}H^{+}, respectively. 

In the text 
Fig. 6 Rate coefficients of the pH_{2}D^{+} + oH_{2}− → oH_{2}D^{+} + pH_{2} (left), oH_{2}D^{+} + pH_{2}− → pH_{2}D^{+} + oH_{2} (middle), and (right) reactions as functions of temperature. Red lines represent groundstatetospecies rate coefficients, while blue lines represent speciestospecies 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 yaxis scaling in the panels. 

In the text 
Fig. A.1 Excitation of the lowest rotational levels of oH_{2}D^{+}, pH_{2}D^{+}, oD_{2}H^{+}, and pD_{2}H^{+} 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 speciestospecies rate coefficients. 

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

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.