Abstract
The triatomic hydrogen ion is one of the most important species for the gas phase chemistry of the interstellar medium. Observations of are used to constrain important physical and chemical parameters of interstellar environments. However, the temperatures inferred from the two lowest rotational states of in diffuse lines of sight – typically the only ones observable – appear consistently lower than the temperatures derived from H observations in the same sightlines. All previous attempts at modelling the temperatures of in the diffuse interstellar medium failed to reproduce the observational results. Here we present new studies, comparing an independent master equation for H level populations to results from the Meudon PDR code for photon dominated regions. We show that the populations of the lowest rotational states of are strongly affected by the formation reaction and that H ions experience incomplete thermalisation before their destruction by free electrons. Furthermore, we find that for quantitative analysis more than two levels of H have to be considered and that it is crucial to include radiative transitions as well as collisions with H2. Our models of typical diffuse interstellar sightlines show very good agreement with observational data, and thus they may finally resolve the perceived temperature difference attributed to these two fundamental species.
1. Introduction
The triatomic hydrogen ion is one of the main drivers of interstellar chemistry in the gas phase [Citation1]. It is formed very efficiently in the interstellar medium by collisions between hydrogen molecules and hydrogen molecular ions (1) (1) Owing to the comparatively low proton affinity of , triatomic hydrogen readily reacts with many of the neutral atomic and molecular species present in interstellar environments. By donating a proton in exothermic ion-neutral collisions of the type (2) (2) triatomic hydrogen often initiates gateway processes, which subsequently enable the formation of more complex molecules in interstellar space (here the X stands for any neutral atomic or molecular collision partner). In particular, reactions between ions and neutral , and atoms will lead to the formation of and , respectively, and thus facilitate the introduction of the heavier atomic species into the chemical networks.
The relevance of for interstellar chemistry was recognised already in early quantitative models of interstellar clouds [Citation2,Citation3]. After the breakthrough work of Oka [Citation4], who identified the infrared spectrum of the fundamental vibrational band in the laboratory, was found in both dense [Citation5] and diffuse lines of sight [Citation6], confirming the role of ion-neutral chemistry in space.
In the meantime, triatomic hydrogen has been detected in various interstellar sightlines, including the Galactic centre, as well as extra-galactic sources and planetary atmospheres (see [Citation7] for a recent review of astronomy). Moreover, owing to its seemingly simple formation and destruction mechanisms, observations have been used to constrain important astrophysical parameters like, e.g. the cosmic ray ionisation rate [Citation8–11] and the temperatures and densities of the molecular gas in the vicinity of the Galactic centre [Citation12]. Submillimeter observations of HD, an isotopic variant of triatomic hydrogen, have been used to infer the minimum age of a star-forming molecular cloud [Citation13].
While the hydrogen molecule is by far the most abundant molecule in space, the bulk of molecules in colder environments is difficult to observe with ground-based telescopes. Direct observation of in diffuse and translucent interstellar clouds are obtained either from satellite absorption observations of electronic transitions in the ultraviolet regime (see, e.g. [Citation14,Citation15]) towards bright stellar sources or by infra-red quadrupolar electric emission of its rovibrational spectrum in bright and dense photodissociation regions (PDRs).
However, attempts to understand the observed column densities and of the two lowest rotational states of (with J = 0 and J = 1, respectively) and the column densities and of the two lowest states of (with and , respectively) in the same lines of sight led to a surprise [Citation16]. The temperature derived from the lowest states of , appeared systematically lower than the temperature of the lowest states of , , which is usually found to be in equilibrium with the gas kinetic temperature [Citation15,Citation16]. While the initial studies used very simple model calculations for the abundances and thermalisation processes, later models, employing large chemical networks, focussed principally on the chemical evolution of the ortho/para forms of and abundancesFootnote1 without considering detailed collisional excitation mechanisms nor introducing thermal balance considerations [Citation17].
Here, we present a different approach and introduce first a master equation describing the evolution of rotational levels of , based on updated rate coefficients for collisional and chemical processes at fixed density and temperature, which can be solved for steady state at definite physical conditions. This procedure is able to reproduce the observational trends and temperature differences on a quantitative level, and it allows for an analysis of the contributions of the various processes. We further introduce the same mechanisms in our Meudon PDR code [Citation18], which solves both the chemical and thermal equilibrium of the cloud and add their contribution to the equilibrium state of , where photodissociation, collisional excitation and chemical formation/destruction mechanisms are considered together.
The paper is organised as follows: Section 2 presents a brief overview of nuclear spin and the rotational levels of and . In Section 3, we describe the most important processes driving the ortho-para ratio of . A master equation for H state populations is presented in Section 4 together with the corresponding results on the excitation temperature of . We compare our model results to astronomical observations in Section 5, where we also introduce the modifications included in the PDR code to account for the ortho-para character of . The paper concludes with a brief discussion in Section 6.
2. Nuclear spin and rotational states of and
Both and exist in two different nuclear spin configurations. For the lowest rotational state of , with J = 0, the proton spins are anti-parallel and add up to I = 0. This configuration is denoted as (or ). The next highest level is J = 1, and, owing to the requirement that the total wave function has to change sign under the permutation of both protons, all states with odd rotational quantum numbers have parallel nuclear spin configurations with I = 1, which is denoted as (or ). Likewise, all even rotational quantum numbers in H can be attributed to with I = 0. Figure shows the three lowest levels of (with ) and their respective energies expressed as (in units of kelvin).
The excitation temperature of the two lowest states of is defined as (3) (3) where stands for the energy difference between the states, is the ratio of the multiplicities of the states with J = 1 and J = 0, respectively, and and denote their populations. An analogous equation can be given for the excitation temperature between the states with J = 2 and J = 0, for which K and , and both temperatures are usually found to be consistent with one another [Citation15], giving a good proxy of the kinetic gas temperature. Typical values for in these diffuse sightlines range from 50 to [Citation15,Citation19–21].
The nuclear spin configurations of are also denoted by para and ortho for and , respectively. As with , the symmetry of the nuclear spin wave function imposes restrictions on the rotational quantum numbers. The relevant quantum number is usually denoted by G, which is connected to the projection of the angular momentum (from both vibrational and rotational motion) onto the molecular symmetry axis (see the review by [Citation22] for quantum numbers, symmetries and selection rules). Since we are only concerned with molecules in the vibrational ground state here, the quantum number G can be regarded as equivalent with the projection K of the rotational angular momentum (denoted by quantum number J) onto the normal axis of the molecular plane, implying G = K. The symmetry of the total wave function requires G = 3n (where n is an integer) for the or ortho levels of , while all other levels (effectively fulfilling ) are of the para configuration, with . The right-hand side of Figure shows all rotational levels with and their respective energies. Because of the Pauli exclusion principle, pure rotational levels with even J and G = 0 do not exist. This rules out the nominal ground state and the state, which are both indicated in the graph by dotted lines for completeness. Since we are almost exclusively concerned with rotational states in the vibrational ground state of , we will from now on refer to all states by giving their quantum numbers only, e.g. and for the lowest para and ortho states, respectively.
With the notable exception of the Galactic centre [Citation12,Citation23,Citation24], only the lowest para-state of and the lowest ortho-state have ever been detected in the interstellar medium. Consequently, the excitation temperature, , derived from observations is usually calculated using the column density ratio of these two states (4) (4) where , and and denote the total degeneracies of the ortho- and para-state, respectively. Most observations yield values for between 20 and [Citation16,Citation17], systematically lower than the excitation temperature . As para- and ortho-states of and can not be inter-converted by radiative transitions, we now describe the various collisional and reactive processes allowing ortho/para exchange.
3. Processes controlling the para fraction of
3.1. Formation of
formation via the + reaction (Equation (Equation1(1) (1) )) has been studied for many years by various experimental techniques. A recent compilation of the results can be found in [Citation25]. Particularly noteworthy is a recent study that employs excited Rydberg molecules in a merged beams approach to reach collision energies between 5 and 60 [Citation26]. The results are in very good agreement with the previous measurements [Citation27] conducted at somewhat higher energies. A recommended fit to the overall cross section can be found in [Citation25], yielding values at low temperature that are slightly higher than the corresponding values derived from the classical Langevin collision rate. We have converted this cross section to a thermal rate coefficient (5) (5) with T, the gas kinetic temperature, in kelvin. This rate shows only a weak dependence on temperature. Its value is slightly larger than the constant rate of that is reported in most astrochemistry databases (and based on a previous experimental study [Citation28], which was conducted at room temperature).
To determine the nuclear spin of the ions created by reaction (Equation1(1) (1) ), we refer to the selection rules outlined by Oka [Citation29]. While this scheme is based on pure angular momentum algebra, and thus does not take any energetic barriers or other restrictions of the reaction into account, we consider this to be a very good approximation for the highly exothermic barrier-less ion neutral reaction between and . To quantify the outcome of the reaction in terms of nuclear spin, it is convenient to use the para-fractions of both molecular species derived from the densities and of the para- and ortho-species of , and and of , respectively. The para-fraction for both species is then denoted by (6) (6) and (7) (7) Following the argumentation by Crabtree et al. [Citation16], we assume that the cosmic ray ionisation of (constituting the main source of ionisation that initiates the formation) does not affect the nuclear spin of the molecule, and therefore the para-fraction of has the same value .
With these assumptions the fraction at formation, , can be derived as a function of , using the nuclear spin branching fractions given in [Citation29] (for details see Table 4 in [Citation16]), resulting in a simple linear dependence of the form (8) (8) The values of and are displayed as a function of temperature in Figure for in thermal equilibrium.
We further assume that the rotational levels of and are each populated at formation according to their Boltzmann distribution at a temperature corresponding to 2/3 of the reaction exothermicity [Citation30]. Given the low energy of the relevant levels, this amounts in effect at using rates proportional to the statistical weight of the level.
3.2. Thermalizing collisions of
Collisions with electrons, , and contribute to the energy exchange between the rovibrational levels of , but only collisions with and may change the nuclear spin of the ions. To the best of our knowledge, of the above collision partners, only for collisions with and electrons detailed information is available in the literature.
3.2.1. Collisions between and
Before any detailed study of collision rates with was available, Oka and Epp [Citation31] suggested to use the Langevin expressions to describe the excitation of detected towards the galacatic centre. However, as a result of the five identical Fermion nuclei involved in these collisions, [Citation29,Citation32,Citation33] pointed out that one should consider the nuclear spin dependence of the total wavefunction. Three different possibilities have to be accounted for (9a) (9a) (9b) (9b) (9c) (9c) Specific selection rules on the nuclear spins of the products can be derived, following [Citation29,Citation34], and they have been used to model hydrogen plasma experiments [Citation16]. The calculations of [Citation32] compared favourably to ion trap measurements of the nuclear spin equilibrium of and at low temperature [Citation35], although these studies did not allow for detailed comparisons of the absolute rate coefficients.
For our models, we have adopted the rate coefficients resulting from the most advanced theoretical treatment of the – reaction so far, which was presented by [Citation36]. Their approach is similar to the method of [Citation32], but refined by a dynamical bias that is introduced through a scrambling matrix, accounting for the relative probabilities of the identity/hop/exchange channels. The probabilities are calculated using quasi-classical trajectory calculations, based on a global potential energy surface [Citation37]. We employ a set of rate coefficients that covers collisions of the lowest 24 rotational states of with ortho- and para- (in their lowest rotational states) and temperatures up to , which were kindly provided by O. Roncero. Those rate coefficients are also currently used in our PDR model calculations [Citation18].
3.2.2. Collisions between , and
and are additional collision partners that should be taken into account when describing the excitation equilibrium of , but, to the best of our knowledge, no detailed studies are presently available on these systems. Collisions of with can not modify the nuclear spin configuration of , and we approximate the corresponding collision rates by taking the rates for and scaling them with the reduced mass factor (while vetoing channels that might change the nuclear spin of ), as described previously [Citation38].
Collisions of with may lead to proton exchange or inelastic collisions, similar to the channels discussed for the reaction. However, since we are not aware of more detailed information on this process, we will use the collisional rates with as a proxy for collisions with . We checked that this choice has no major influence on our model calculations for the conditions studied here.
3.2.3. Inelastic collisions between and electrons
Electronic collisions with preserve the ortho-para character of . They have been computed by [Citation39] and are included in the present model. Until very recently no laboratory measurements of the change of rotational states in electron collisions were available for any molecular ion to benchmark the theoretical approach. But a recent measurement of low-energy electron collisions with allowed for a comparison between experiment and theory for the lowest rotational states [Citation40], which revealed very good agreement.
3.3. Dissociative recombination of
In diffuse gas, the main destruction reaction of is the dissociative recombination (DR) with electrons. (10a) (10a) (10b) (10b) This is an important chemical reaction for interstellar chemistry, and as such, has attracted a lot of attention, as well as some controversy. More than 30 independent experimental studies of the DR rate coefficient of have been published, with outcomes that differ by orders of magnitude (there are a number of reviews on this topic, see, e.g, [Citation41–43]). In summary, the present consensus is that the absolute rate of the DR rate coefficient is correctly derived from the storage ring merged beams measurements [Citation43]. Theoretical studies identified the Jahn-Teller effect as a driver for the recombination process at low temperatures [Citation44–46].
Both astrochemistry databases, UMIST Footnote2 and KIDAFootnote3, report rate coefficients with branching ratios of 2/3 for reaction (Equation10a(10a) (10a) ) and 1/3 for reaction (Equation10b(10b) (10b) ), based on storage ring experiments, and the total DR reaction rate coefficient is given as (11) (11) However, there are recent suggestions concerning a possible difference between the DR rate coefficients of the two nuclear spin modifications of . These were first pointed out in the storage ring experiments of [Citation47,Citation48], but the values are dependent on the actual rotational level populations, which could not be determined precisely. Subsequent plasma experiments found an even stronger effect at low temperature [Citation49], and updated theoretical studies [Citation46] predict that at low temperature the difference between the rate coefficients for the two nuclear spin modifications may exceed an order of magnitude, with recombining much faster than . Two different theoretical values are explicitly reported in [Citation50], supporting plasma studies. We have derived an analytic expression from these results, which we can describe by the formulae (12a) (12a) (12b) (12b)
Table summarises the values of the principal reaction rate coefficients introduced in the present study, where we have assumed the same branching ratios for reactions (Equation10a(10a) (10a) ) and (Equation10b(10b) (10b) ) as described above.
4. Master equation for state populations
In this section, we describe the master equation for the level populations of , and we show that chemistry and specific molecular properties result in an excitation temperature that is lower than the kinetic temperature of the gas.
The demonstration starts from a full treatment of the differential state equations, including all possible processes, i.e. collisional and radiative transitions as well as chemical state-to-state formation and destruction reactions. Solving these equations, the next section (4.1) will show that we can indeed recover a value for the excitation temperature that is systematically below the gas kinetic temperature, similar to observational results. The next sub-section explores how sensitive these equations are to the number of levels included in the computation. This allows to understand why at least 5 levels must be included for quantitative results. It also shows that the range is much less sensitive to the number of levels included. By restricting our analysis to this temperature range, we can derive a qualitative analytical approximation with only two levels (Section 4.3). The resulting expression shows that selective formation of p- by p- followed by incomplete thermalisation is responsible for the deviation from thermal equilibrium.
4.1. Differential equations
The variation over time of the density of a level i () of , , can be described by (13) (13) with and denoting the state-dependent chemical formation and destruction rates, the collisional excitation/de-excitation rates with and , and denote collision rates where X stands for , , . are the radiative emission transition probabilities.
We underline two important points:
The state-dependent chemical formation and destruction rates are introduced in the master equation and have to be evaluated. It is important to note that the formation rate of a particular level can differ from its destruction rate (see Section 3.3).
The main formation process of is a highly exothermic reaction, which preferentially populates high energy levels. Thus, chemical formation can be regarded as an excitation mechanism.
Introducing , the relative populations of , so that , the total destruction rate is . The total formation rate is . Since the relative populations are unknown until the differential equations are solved, it is not possible to compute the destruction rate beforehand, as long as the differ for individual states, and thus cannot be computed independently. At steady state, and the system of equations to solve is (14) (14) with . The value of R is initially unknown, as the total destruction rate of requires the knowledge of the relative level populations of the molecular ion. So, we explicitly add the conservation equation (15) (15) We get a system of N + 1 equations with N + 1 unknowns, which is easily solved if the densities of , , and the temperature are known or assumed. All other quantities are rate coefficients or parameters that are derived from experiment or theory.
4.1.1. Derivation of density and electronic fraction
It is possible to reduce the set of equations further by estimating and . The main formation and destruction reactions of are: (16a) (16a) (16b) (16b) (16c) (16c) where represents cosmic ray particles and ζ the corresponding ionisation rate of . At steady state, this leads to (with the dissociative recombination rate of ) (17) (17) This removes one parameter from the system.
As for , in diffuse gas it is often assumed to be equal to the density of . However, for high cosmic ray ionisation rates, as, e.g. in the Central Molecular Zone (CMZ) of our galaxy, Le Petit et al. [Citation11] have shown that protons may contribute significantly to the ionisation fraction. Here we compute the electronic density considering both and as described in Appendix 2, which involves the relative abundances of all atoms with ionisation potential below the Lyman cutoff (assumed to be fixed) and the UV radiation field Footnote4. With these derivations of and , the system formed by Equations (Equation14(14) (14) ) and (Equation15(15) (15) ) depends on five astrophysical parameters: Footnote5, T, , , and ζ, where we introduced the molecular fraction, .
4.2. Numerical solution of the coupled equations
We developed ExcitH3p, a FORTRAN code that solves the coupled system of equations (Equation (Equation14(14) (14) )) for any number of levels, and then computes the two-level excitation temperature (Equation (Equation4(4) (4) )) (18) (18) Here and denote the relative populations of the ground state and the first excited state , respectively.
In practice, the 24 energetically lowest rotational levels of are included in the model. Those are the levels for which radiative emission rates as well as collisional excitation/de-excitation rates are known or have been estimated. The highest level in this framework is the ortho level, located above the ground state. We find that the solution to the set of equations (Equation14(14) (14) ) recovers very nicely the full results obtained with the Meudon PDR code for diffuse line of sight conditions. The advantage of the master equation approach is that it is much less computationally expensive, and allows for a rapid exploration of the parameter space and the various possible hypotheses concerning the less well-known physical processes.
Results derived by the FORTRAN routine are presented in Figure for a range of typical values of the cosmic ionisation rate ζ and gas kinetic temperature T. The total gas density is set to and the molecular fraction to . We find that in the entire parameter space, , the excitation temperature given by the two lowest levels, is systematically lower than the kinetic temperature. The overall magnitude of the discrepancy between the two temperatures is in agreement with the observations [Citation10,Citation17]. This outcome is a sole property of the microscopic excitation and de-excitation mechanisms of the individual quantum levels of , which are detailed above. Figure also shows that a maximum of about is reached for at kinetic temperatures around . The occurrence of such a maximum is remarkable. We have verified that the maximum is still there, although somewhat different, when the dissociative recombination rate coefficients of para and ortho- are the same: = 38 K for T = 150 K under the same physical conditions. We come back to that point in Section 4.2.1 when discussing the influence of the number of levels included in the model.
Figure shows the influence of density, showing that is not sensitive to this parameter up to densities of more than if the gas temperature is below ∼ 80K.
4.2.1. Influence of the number of levels included in the model
Interesting conclusions can be drawn from the examination of the dependence of our results on the number N of levels included in the model. Except for the Galactic centre, only the two lowest rotational levels of H have been detected in the ISM so far. Consequently, the interpretation of astrophysical observations is usually restricted to these two levels. Figure displays the computed value as a function of the kinetic temperature T, for different values of N. We perform this comparison for the following parameters: , , and . We find that the curve corresponding to N = 2 is far from the values obtained with 24 levels, except for a small range at very low temperatures. Nevertheless, even with N = 2, is systematically below the kinetic temperature. N = 5 is the minimum needed for acceptable results, and convergence is reached for N = 10 for the physical conditions considered here. It is important to realise that N = 5 involves the metastable level, representing a possible sink for population.
The decrease of with increasing gas kinetic temperature above seems puzzling at first. However, it can be understood by considering the various possible decay paths from the two excited para levels (corresponding to relative population in our enumeration) and (corresponding to ), and the lowest excited ortho level (corresponding to ). It is important to note that the level is metastable with an infinite radiative lifetime. Hence it can only be depopulated in collisional processes. Examination of the collision rates with shows that its branching ratios towards ortho and para are approximately equal in the 30–300 temperature range. Both para levels and , on the other hand, can decay both radiatively and collisionally. The critical densities of these levels, given by the ratio of the radiative emission rate and the total collisional de-excitation rate coefficient, are a few thousand cubic centimeters. Thus, at densities of a few hundred cubic centimeters, relevant for the astrophysical observations discussed here, radiative decay of para levels to the ground state is much more likely than collisional de-excitation, while the ground ortho level (1,0) is underpopulated compared to a thermal Boltzmann distribution, as population may be trapped in the level. This leads to a perceived overpopulation of the lowest para state, compared to the lowest ortho state. Chemical formation can populate efficiently all levels due to the large exothermicity of the reaction. So, increasing N leads to more open channels to populate levels that decay to and inhibits population of the lowest ortho state .
4.3. Two-level approximation
In a restricted range of kinetic temperatures around , the two-level case is a fair approximation to the full system, as seen in Figure . It allows to considerably simplify the system and to derive an analytic expression for that offers the opportunity to highlight the key microscopic mechanisms at work.
We restrict our study to molecular hydrogen collisions and introduce , the total collisional excitation rate due to from level 1 to level 2 (and the equivalent for ). Then, the coupled equations system (Equation (Equation14(14) (14) )) reduces to: (19a) (19a) (19b) (19b)
Introducing , we can compute the ratio: (20) (20) The factor R cancels out and the ratio does not depend on . We introduce the electronic fraction and the configuration-specific formation rates of H with and . Furthermore, we apply detailed balance to the H–H collisional rates, yielding . Finally, we get the following expression (21) (21) Apart from the kinetic temperature T, the ratio depends on the collisional and dissociative recombination rate coefficients, the molecular fraction of , and the electronic fraction . It is independent of the cosmic ray ionisation rate ζ and of the total formation reaction rate coefficient of H, but – crucially – not of the branching ratios, which effectively enter the equation through the H para fraction .
We can interpret the term in square brackets as a correction factor to the Boltzmann value for the kinetic temperature. At the low temperatures considered here (25–50) the value of ranges from 0.6 to 1 (see Figure ), and numerical analysis reveals that the and pre-factors in the nominator and the denominator are sufficient to keep the overall correction factor smaller than unity in the entire temperature range, resulting in a lowering of the ratio compared to the thermal value (in agreement with astronomical observations). The extent of the deviation from thermal equilibrium depends critically on the ratios and of the electron recombination rates to the collisional rate coefficients. If we, for the sake of argument, extrapolate the formula using an artificially large value for – corresponding to very efficient collisional thermalisation – the entire correction factor tends toward unity, and we will recover the ratios given by the gas kinetic temperature. The same is obviously true for very small DR rate coefficients and .
The picture that emerges is that the excitation temperature of appears lower than the nominal gas temperature because of incomplete thermalisation. The formation process strongly favours the formation of at low temperatures, as the para-fraction of is large. The electron recombination process then removes H before it can reach thermal equilibrium in collisions with . While we stress that these conclusions based on the two-level approximation are only valid in a limited temperature range, and that for quantitative results more levels need to be considered, we can reproduce the general trend very well using the master equation approach described above. For artificially enlarged collisional rate coefficients (or sufficiently reduced electron recombination rates) the calculations reach complete thermal equilibrium between the excitation temperature and the kinetic temperature T.
In essence, the nuclear spin restrictions of the formation reaction produce an over-proportional amount of in the para configuration, and thermalisation in collisions with H is too slow to reach equilibrium before the ions are destroyed by free electrons. This is to be contrasted to the situation for , where the destruction process is much slower compared to thermalising collisions (see Appendix 1 for details).
5. Comparison to observations
5.1. Master equation approach
To validate our approach, we compare the results of the master equation approach (Equation (Equation14(14) (14) )) with observations. Table presents the excitation temperature reported in the literature for a few local diffuse clouds as well as the corresponding excitation temperatures. Data for come from Copernicus [Citation19] and FUSE [Citation20,Citation21] satellite observations. The proton densities reported in this table are those published in the papers reporting the data. Determination of diffuse cloud density with observations of and is not straightforward, because a fraction of the hydrogen atoms observed on the line of sight may not be related to the diffuse clouds to which belongs. So, these densities are to be seen as order of magnitude estimates, and should be considered only as indicative.
Using the ExcitH3p program that solves Equation (Equation14(14) (14) ), we compute for all the lines of sight. We assume a cosmic ray ionisation rate of and a molecular fraction of . For the gas density, the values in Table are used, and for the gas temperature we use the values derived from observations .
5.1.1. Model results and sensitivity to the DR rate coefficients
Figure shows the results for three different hypothesis for , and the value provided by Equation (Equation12a(12a) (12a) ) for . We see a quasi-linear variation of with , which is well reproduced by the models. Using from Equation (Equation12b(12b) (12b) ) (red circles) leads to temperatures which are too high, while using the same rate for and (pink circles) leads to temperatures which are too low. An empirical adjustment using from Equation (Equation12a(12a) (12a) ) and (green circles) gives a very satisfying result, given that no attempt has been made to optimise other parameters.
There is a single noticeable exception: Per. However, of this line of sight suffers from a particularly large error bar as listed in Table . A high excitation temperature of can only be reached for a kinetic temperatures close to with our models, as can be seen in Figure .
We note that typically 10 to of ions are in excited stable or metastable levels above the respective lowest ortho or para levels, such as , , , and .
5.2. Full PDR model
For comparison and validation, we compare the results of Section 5.1 to those obtained with our full PDR code for the case of ζ Per. We do not seek a complete model of this line of sight, and do not try to optimise the free parameters to reproduce other observations than and excitation.
We use our previous study dedicated to ζ Per [Citation52] as a starting point, but we also account for the many updates made since then [Citation11,Citation18,Citation53]. For the present study, we have introduced in the Meudon PDR code the ortho/para dependence of the formation reaction of via + collisions and the nuclear spin dependence of the dissociative recombination rate coefficient, in addition to the other excitation mechanisms of . The collisional excitation of by , that has been revisited by [Citation54] for highly rovibrationally excited levels, has also been updated. Observations suggest a total visual extinction (, , ). In our models the cloud is illuminated from both sides by the standard ISRF ().
We consider two different scenarios
Constant density and constant temperature,
Constant density, with computation of the thermal balance.
Both models A and B are fairly standard. To account for the presence of purely atomic gas along the line of sight, we use an for both models A and B, allowing to account for the molecular hydrogen abundance.
The results are summarised in Table . The examples shown here have been selected from an evaluation of the following , where σ are the observational uncertainties (22) (22)
The proposed ionisation rates are rather high, but in line with the other recent evaluations based on and abundances [Citation9,Citation10,Citation55]. We note that our previous estimate of the cosmic ionisation rate towards ζ Per was somewhat lower, as we included in that study additional constraints provided by and column densities. The constraints imposed by the excitation temperature of and are not sensitive to the cosmic ionisation rate of , as shown in Figure . We then recover the predictions based on molecular ion observations.
Using the recombination rate from Equation (Equation12b(12b) (12b) ), the excitation temperature of is slightly over-estimated, as found in the previous Section. Applying the same empirical correction to the recombination rate leads to a lower excitation temperature, without any impact on . The total amount of is lower due to the overall larger recombination rate. The resulting varies accordingly. This illustrates the very high sensitivity of to the exact value of the electron recombination rate. The situation now is much better than it used to be, but smaller uncertainties are still needed for quantitative analysis. Note that we do not claim to estimate these rates from observational data.
Overall, this comparison validates the excitation model presented in Section 4 and shows that the temperatures of and in the diffuse ISM can, in fact, be predicted considering a a strongly reduced set of reactions.
6. Discussion and conclusion
In this paper, we review all physical processes relevant to formation, excitation and destruction of in the diffuse interstellar medium. We provide references to the best data available to date.
We show that a statistical model of level populations, including formation/destruction terms and updated collisional excitation/deexcitation processes, allows to explain the low excitation temperature observed in diffuse clouds that has puzzled observers and modellers alike for twenty years. In particular, we show that it is mandatory to include state-to-state chemical formation and destruction processes to explain the departure from thermal equilibrium (Boltzmann ratio) observed for the two lowest levels of . Specifically, the formation of by at temperatures below is efficient. Considering the individual levels of , we find that reactive collisions with are generally too slow – when compared to spontaneous radiative decay and the fast destruction by electron recombination – to bring the populations into equilibrium with the gas kinetic temperature. These results are confirmed by an updated version of the Meudon PDR code that includes the specific ortho/para-dependence of the formation/destruction reactions of in addition to the radiative/collisional excitation balance of that molecular ion. While the formation process may be primarily responsible for the increased population of p-H3+ levels at low temperature, it is important to note that the consideration of different classes of processes – radiative transitions as well as chemical reactions with H2 – is required to achieve quantitative results. We find that all attempts to simplify our master equation further and remove more processes from our models lead to significant changes in the H3+ excitation temperature and impair the agreement with the observational data.
Moreover, we show that the inclusion of rotationally excited levels, besides the respective ortho and para ground states that are usually considered, has substantial implications for the population of the first two levels and thus for , and we find that at least 10 levels should be included in the coupled equations to get a converged result. Our models suggest that typically more than of ions are in metastable excited states, which may be observable in absorption towards bright stars or quasi-stellar objects.
Finally, we stress that some key processes are still either badly determined or completely unknown. In particular, precise quantitative computation of excitation will not be possible as long as we still lack accurate state specific recombination rates with electrons at temperatures between 20 and . Numerical manipulations of the rate coefficients shown in Section 5 reveal the sensitivity of to the ratio . However, it is well-known how dangerous it can be to try to infer reaction rates from observational results, and we do not claim that the rates used here are the final word. Another class of processes that are lacking accurate description are collision rates of with and . In particular, exchange of hydrogen atoms during collisions with may impact the ortho to para ratio of . Here we used estimated values for the rate coefficients scaled from reaction rates with in order to include these processes in the models. While our results seem not to depend strongly on the exact choice of the estimated rate coefficients, more accurate values for these reactions are clearly desirable.
Despite the remaining limitations, our models for the first time are able to account for the observed excitation temperature in diffuse cloud sightlines. This marks a major step in our understanding of interstellar hydrogen chemistry, providing a framework of state-selective chemistry for two of the most important and fundamental molecular gas phase species.
Disclosure statement
No potential conflict of interest was reported by the author(s).
Correction Statement
This article has been corrected with minor changes. These changes do not impact the academic content of the article.
Additional information
Funding
Notes
1 In fact, only the lowest para and ortho levels are considered in these studies.
2 Available at http://www.udfa.net.
3 Available at https://kida.astrochem-tools.org.
4 The UV radiation field strength is defined in [Citation51]. It controls the interstellar grain charge, which impacts the recombination of and .
5 , the proton density is expressed as = .
References
- T. Oka, Proc. Natl. Acad. Sci. 103 (33), 12235–12242 (2006). doi:10.1073/pnas.0601242103
- W.D. Watson, Astrophys. J. Lett. 183, L17 (1973). doi:10.1086/181242
- E. Herbst and W. Klemperer, Astrophys. J. 185, 505–534 (1973). doi:10.1086/152436
- T. Oka, Phys. Rev. Lett. 45 (7), 531–534 (1980). doi:10.1103/PhysRevLett.45.531
- T.R. Geballe and T. Oka, Nature 384 (6607), 334–335 (1996). doi:10.1038/384334a0
- B.J. McCall, T.R. Geballe, K.H. Hinkle and T. Oka, Science 279, 1910 (1998). doi:10.1126/science.279.5358.1910
- S. Miller, J. Tennyson, T.R. Geballe and T. Stallard, Rev. Mod. Phys. 92, 035003 (2020). doi:10.1103/RevModPhys.92.035003
- B.J. McCall, A.J. Huneycutt, R.J. Saykally, T.R. Geballe, N. Djuric, G.H. Dunn, J. Semaniak, O. Novotny, A. Al-Khalili, A. Ehlerding, F. Hellberg, S. Kalhori, A. Neau, R. Thomas, F. Österdahl and M. Larsson, Nature 422 (6931), 500–502 (2003). doi:10.1038/nature01498
- N. Indriolo, G.A. Blake, M. Goto, T. Usuda, T. Oka, T.R. Geballe, B.D. Fields and B.J. McCall, Astrophys. J. 724 (2), 1357–1365 (2010). doi:10.1088/0004-637X/724/2/1357
- N. Indriolo and B.J. McCall, Astrophys. J. 745 (1), 91 (2012). doi:10.1088/0004-637X/745/1/91
- F. Le Petit, M. Ruaud, E. Bron, B. Godard, E. Roueff, D. Languignon and J. Le Bourlot, Astron. Astrophys. 585, A105 (2016). doi:10.1051/0004-6361/201526658
- T. Oka, T.R. Geballe, M. Goto, T. Usuda, J. M. Benjamin and N. Indriolo, Astrophys. J. 883 (1), 54 (2019). doi:10.3847/1538-4357/ab3647
- S. Brünken, O. Sipilä, E.T. Chambers, J. Harju, P. Caselli, O. Asvany, C.E. Honingh, T. Kamiński, K.M. Menten, J. Stutzki and S. Schlemmer, Nature 516 (7530), 219–221 (2014). doi:10.1038/nature13924
- J. Spitzer, L. W.D. Cochran and A. Hirshfeld, Astrophys. J. Supp. 28, 373–389 (1974). doi:10.1086/190323
- J.M. Shull, C.W. Danforth and K.L. Anderson, Astrophys. J. 911 (1), 55 (2021). doi:10.3847/1538-4357/abe707
- K.N. Crabtree, N. Indriolo, H. Kreckel, B.A. Tom and B.J. McCall, Astrophys. J. 729 (1), 15 (2011). doi:10.1088/0004-637X/729/1/15
- T. Albertsson, N. Indriolo, H. Kreckel, D. Semenov, K.N. Crabtree and T. Henning, Astrophys. J. 787 (1), 44 (2014). doi:10.1088/0004-637X/787/1/44
- F. Le Petit, C. Nehmé, J. Le Bourlot and E. Roueff, Astrophys. J. 164 (2), 506–529 (2006). doi:10.1086/apjs.2006.164.issue-2
- B.D. Savage, R.C. Bohlin, J.F. Drake and W. Budich, Astrophys. J. 216, 291–307 (1977). doi:10.1086/155471
- B.L. Rachford, T.P. Snow, J. Tumlinson, J.M. Shull, W.P. Blair, R. Ferlet, S.D. Friedman, C. Gry, E.B. Jenkins, D.C. Morton, B.D. Savage, P. Sonnentrucker, A. Vidal-Madjar, D.E. Welty and D.G. York, Astrophys. J. 577 (1), 221–244 (2002). doi:10.1086/apj.2002.577.issue-1
- B.L. Rachford, T.P. Snow, J.D. Destree, T.L. Ross, R. Ferlet, S.D. Friedman, C. Gry, E.B. Jenkins, D.C. Morton, B.D. Savage, J.M. Shull, P. Sonnentrucker, J. Tumlinson, A. Vidal-Madjar, D.E. Welty and D.G. York, Astrophys. J. Supp. 180 (1), 125–137 (2009). doi:10.1088/0067-0049/180/1/125
- C.M. Lindsay and B.J. McCall, J. Mol. Spectr. 210 (1), 60–83 (2001). doi:10.1006/jmsp.2001.8444
- M. Goto, B.J. McCall, T.R. Geballe, T. Usuda, N. Kobayashi, H. Terada and T. Oka, Publ. Astron. Soc. Jpn. 54, 951–961 (2002). doi:10.1093/pasj/54.6.951
- T. Oka and T.R. Geballe, Astrophys. J. 927 (1), 97 (2022). doi:10.3847/1538-4357/ac3912
- I. Savić, S. Schlemmer and D. Gerlich, Chem. Phys. Chem. 21 (13), 1429–1435 (2020). doi:10.1002/cphc.v21.13
- P. Allmendinger, J. Deiglmayr, O. Schullian, K. Höveler, J.A. Agner, H. Schmutz and F. Merkt, Chem. Phys. Chem. 17 (22), 3596–3608 (2016). doi:10.1002/cphc.v17.22
- T. Glenewinkel-Meyer and D. Gerlich, Isr. J. Chem. 37 (4), 343–352 (1997). doi:10.1002/ijch.v37.4
- L.P. Theard and W.T. Huntress, J. Chem. Phys. 60 (7), 2840–2848 (1974). doi:10.1063/1.1681453
- T. Oka, J. Mol. Spectr. 228 (2), 635–639 (2004). doi:10.1016/j.jms.2004.08.015
- F. Merkt, K. Höveler and J. Deiglmayr, J. Phys. Chem. Letters 13, 864 (2022). doi:10.1021/acs.jpclett.1c03374
- T. Oka and E. Epp, Astrophys. J. 613 (1), 349–354 (2004). doi:10.1086/apj.2004.613.issue-1
- K. Park and J.C. Light, J. Chem. Phys. 126 (4), 044305 (2007). doi:10.1063/1.2430711
- E. Hugo, O. Asvany and S. Schlemmer, J. Chem. Phys. 130 (16), 164302 (2009). doi:10.1063/1.3089422
- M. Quack, Mol. Phys. 34 (2), 477–504 (1977). doi:10.1080/00268977700101861
- F. Grussie, M.H. Berg, K.N. Crabtree, S. Gärtner, B.J. McCall, S. Schlemmer, A. Wolf and H. Kreckel, Astrophys. J. 759 (1), 21 (2012). doi:10.1088/0004-637X/759/1/21
- S. Gómez-Carrasco, L. González-Sánchez, A. Aguado, C. Sanz-Sanz, A. Zanchet and O. Roncero, J. Chem. Phys. 137 (9), 094303 (2012). doi:10.1063/1.4747548
- A. Aguado, P. Barragán, R. Prosmiti, G. Delgado-Barrio, P. Villarreal and O. Roncero, J. Chem. Phys. 133 (2), 024306 (2010). doi:10.1063/1.3454658
- E. Roueff and F. Lique, Chem. Rev. 113 (12), 8906–8938 (2013). doi:10.1021/cr400145a
- V. Kokoouline, A. Faure, J. Tennyson and C.H. Greene, Mon. Not. R. Astron. Soc. 405 (2), 1195–1202 (2010). doi:10.1111/j.1365-2966.2010.16522.x
- A. Kálosi, M. Grieser, R. von Hahn, U. Hechtfischer, C. Krantz, H. Kreckel, D. Müll, D. Paul, D.W. Savin, P. Wilhelm, A. Wolf and O. Novotný, Phys. Rev. Lett. 128, 183402 (2022). doi:10.1103/PhysRevLett.128.183402
- M. Larsson, B. McCall and A. Orel, Chem. Phys. Lett. 462 (4), 145–151 (2008). doi:10.1016/j.cplett.2008.06.069
- M. Larsson, Phil. Trans. R. Soc. A 370 (1978), 5118–5129 (2012). doi:10.1098/rsta.2012.0020
- H. Kreckel, A. Petrignani, O. Novotný, K. Crabtree, H. Buhr, B.J. McCall and A. Wolf, Phil. Trans. R. Soc. A 370 (1978), 5088–5100 (2012). doi:10.1098/rsta.2012.0019
- V. Kokoouline, C.H. Greene and B.D. Esry, Nature 412 (6850), 891–894 (2001). doi:10.1038/35091025
- V. Kokoouline and C.H. Greene, Phys. Rev. Lett. 90 (13), 133201 (2003). doi:10.1103/PhysRevLett.90.133201
- S.F. dos Santos, V. Kokoouline and C.H. Greene, J. Chem. Phys. 127 (12), 124309 (2007). doi:10.1063/1.2784275
- H. Kreckel, M. Motsch, J. Mikosch, J. Glosík, R. Plašil, S. Altevogt, V. Andrianarijaona, H. Buhr, J. Hoffmann, L. Lammich, M. Lestinsky, I. Nevo, S. Novotny, D.A. Orlov, H.B. Pedersen, F. Sprenger, A.S. Terekhov, J. Toker, R. Wester, D. Gerlich, D. Schwalm, A. Wolf and D. Zajfman, Phys. Rev. Lett. 95 (26), 263201 (2005). doi:10.1103/PhysRevLett.95.263201
- H. Kreckel, O. Novotný, K.N. Crabtree, H. Buhr, A. Petrignani, B.A. Tom, R.D. Thomas, M.H. Berg, D. Bing, M. Grieser, C. Krantz, M. Lestinsky, M.B. Mendes, C. Nordhorn, R. Repnow, J. Stützel, A. Wolf and B.J. McCall, Phys. Rev. A 82 (4), 042715 (2010). doi:10.1103/PhysRevA.82.042715
- J. Varju, M. Hejduk, P. Dohnal, M. Jílek, T. Kotrík, R. Plašil, D. Gerlich and J. Glosík, Phys. Rev. Lett. 106 (20), 203201 (2011). doi:10.1103/PhysRevLett.106.203201
- L. Pagani, C. Vastel, E. Hugo, V. Kokoouline, C.H. Greene, A. Bacmann, E. Bayet, C. Ceccarelli, R. Peng and S. Schlemmer, Astron. Astrophys. 494 (2), 623–636 (2009). doi:10.1051/0004-6361:200810587
- B.T. Draine, Astrophys. J. Supplt. 36, 595–619 (1978). doi:10.1086/190513
- F. Le Petit, E. Roueff and E. Herbst, Astron. Astrophys. 417, 993–1002 (2004). doi:10.1051/0004-6361:20035629
- J. Le Bourlot, F. Le Petit, C. Pinto, E. Roueff and F. Roy, Astron. Astrophys. 541, A76 (2012). doi:10.1051/0004-6361/201118126
- T. González-Lezana, P. Hily-Blant and A. Faure, J. Chem. Phys. 154 (5), 054310 (2021). doi:10.1063/5.0039629
- X.L. Bacalla, H. Linnartz, N.L.J. Cox, J. Cami, E. Roueff, J.V. Smoker, A. Farhang, J. Bouwman and D. Zhao, Astron. Astrophys. 622, A31 (2019). doi:10.1051/0004-6361/201833039
- A. Dalgarno, J.H. Black and J.C. Weisheit, Astrophys. Lett. 14, 77 (1973).
- J.L. Spitzer and W.D. Cochran, Astrophys. J. Lett. 186, L23 (1973). doi:10.1086/181349
- A. Faure, P. Hily-Blant, C. Rist, G. Pineau des Forêts, A. Matthews and D.R. Flower, Monthl. Not. R. Astr. Soc. 487 (3), 3392–3403 (2019). doi:10.1093/mnras/stz1531
- A. Sternberg, F. Le Petit, E. Roueff and J. Le Bourlot, Astrophys. J. 790 (1), 10 (2014). doi:10.1088/0004-637X/790/1/10
- B.T. Draine, Physics of the Interstellar and Intergalactic Medium (Princeton University Press, 2011).
- N.R. Badnell, Astrophys. J. Supp. 167 (2), 334–342 (2006). doi:10.1086/apjs.2006.167.issue-2
- J.C. Weingartner and B.T. Draine, J. Astrophys. 563 (2), 842–852 (2001). doi:10.1086/apj.2001.563.issue-2
Appendices
Appendix 1. Two-level approximation for
As first recognised by [Citation56,Citation57] from Copernicus observations, of is an excellent proxy for the kinetic temperature T in diffuse clouds, where collisions with protons allow the two rotational levels to reach thermal equilibrium, in absence of any radiative transition. This is not true anymore in dense cloud conditions, where the ortho-to-para ratio is expected to be very far from thermal equilibrium [Citation58], and where collisions with modify the excitation balance. We discuss rapidly the conditions of validity of this feature through a simple 2-level approximation. Transitions between J = 0 and J = 1 occur only through reactive collisions with with rates and (reactive collisions with other species (, ) are negligible here). Besides collisions, these two levels are populated by direct formation on grains with rates and or depopulated by photodissociation with rates and . Other chemical reactions have only a minor impact. The resulting balance equations are This system can be solved for and and leads to In the temperature range from 50 to appropriate do diffuse and translucent cloud conditions, the Boltzmann factor varies from 0.5 to 1.6. So, and remain close to one another. Furthermore, the formation rates and are close to one another and simplify. The ratio can be arranged using detailed balance as: In regions of low radiation field where most diffuse clouds are found the transition is very close to the edge of the cloud, as shown in [Citation59]. Hence, in most of the region that builds column density the dissociation rates and are about 4 orders of magnitude lower than the products and . Thus, the correction factor coming from the chemistry is very close to 1 and the ratio gives a very good measure of the kinetic temperature.
Appendix 2
Computation of the electronic fraction
Ionization balance can be solved analytically for diffuse cloud conditions. Due to ultraviolet photons, all metals with an ionisation threshold below are ionised, providing a minimal electronic abundance of , where is the fraction of relevant metals (mostly and , with traces of and other heavier species). In the following, we take .
Additional electrons come from and resulting from the balance between ionisation via cosmic rays and recombination with electrons and grains. We follow here for the most part the presentation of [Citation60], Section 13.6, extended to include . The balance equations are Where is the radiative recombination rate, and the rate of recombination on grains. and are the cosmic ray ionisation rates of and , respectively. With respect to , we use , including secondary ionisation, and . The total abundance of electrons is For all abundances, we write . We take , with and compute abundance from the molecular fraction The resulting system of two equations can be written in a more compact form by using X for hydrogen and Y for helium This system is easily solved using a Newton-Raphson scheme once the rates are known.
Radiative electronic recombination rates are taken from [Citation61], the relevant coefficients are given in Table with
Electronic recombination on grains, comes from [Citation62] with Here, is the Inter Stellar Radiation Field (ISRF) intensity in units of Draine's ISRF. Coefficients to are given in Table ().
Appendix 3
Fortran code
The code solving for populations from Equation (Equation14(14) (14) ) is available at Meudon ISM Services Platform. Compilation requires a modern Fortran 90 compiler (gfortran will do) and access to the LAPACK library. The later is usually provided with all standard compilers. Otherwhile, it is self contained.
The code takes very few input parameters: and the number of levels used, from the command line or redirection of a small input file. It uses the latest data available, as described in this paper. Comments in the source file, coupled to this paper, should be enough for easy use and adaptation.