Understanding the temperatures of H₃⁺ and H₂ in diffuse interstellar sightlines

Abstract

The triatomic hydrogen ion ${\mathrm{H}}_3^{+}$ is one of the most important species for the gas phase chemistry of the interstellar medium. Observations of ${\mathrm{H}}_3^{+}$ are used to constrain important physical and chemical parameters of interstellar environments. However, the temperatures inferred from the two lowest rotational states of ${\mathrm{H}}_3^{+}$ in diffuse lines of sight – typically the only ones observable – appear consistently lower than the temperatures derived from H${}_2$ observations in the same sightlines. All previous attempts at modelling the temperatures of ${\mathrm{H}}_3^{+}$ in the diffuse interstellar medium failed to reproduce the observational results. Here we present new studies, comparing an independent master equation for H${}_3^{+}$ level populations to results from the Meudon PDR code for photon dominated regions. We show that the populations of the lowest rotational states of ${\mathrm{H}}_3^{+}$ are strongly affected by the formation reaction and that H${}_3^{+}$ ions experience incomplete thermalisation before their destruction by free electrons. Furthermore, we find that for quantitative analysis more than two levels of H${}_3^{+}$ have to be considered and that it is crucial to include radiative transitions as well as collisions with H₂. 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.

GRAPHICAL ABSTRACT

Keywords:

• Astrochemistry
• interstellar medium
• triatomic hydrogen
• molecular clouds

1. Introduction

The triatomic hydrogen ion ${\mathrm{H}}_3^{+}$ is one of the main drivers of interstellar chemistry in the gas phase [1]. It is formed very efficiently in the interstellar medium by collisions between hydrogen molecules and hydrogen molecular ions (1) ${\mathrm{H}}_2+{\mathrm{H}}_2^{+}⟶{\mathrm{H}}_3^{+}+\mathrm{H}.$(1) Owing to the comparatively low proton affinity of ${\mathrm{H}}_2$, 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) ${\mathrm{H}}_3^{+}+\mathrm{X}⟶{\mathrm{XH}}^{+}+{\mathrm{H}}_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 ${\mathrm{H}}_3^{+}$ ions and neutral $\mathrm{O}$, and $\mathrm{C}$ atoms will lead to the formation of ${\mathrm{OH}}^{+}$ and ${\mathrm{CH}}^{+}$, respectively, and thus facilitate the introduction of the heavier atomic species into the chemical networks.

The relevance of ${\mathrm{H}}_3^{+}$ for interstellar chemistry was recognised already in early quantitative models of interstellar clouds [2,3]. After the breakthrough work of Oka [4], who identified the infrared spectrum of the ${\mathrm{H}}_3^{+}$ fundamental vibrational band in the laboratory, ${\mathrm{H}}_3^{+}$ was found in both dense [5] and diffuse lines of sight [6], 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 [7] for a recent review of ${\mathrm{H}}_3^{+}$ astronomy). Moreover, owing to its seemingly simple formation and destruction mechanisms, ${\mathrm{H}}_3^{+}$ observations have been used to constrain important astrophysical parameters like, e.g. the cosmic ray ionisation rate [8–11] and the temperatures and densities of the molecular gas in the vicinity of the Galactic centre [12]. Submillimeter observations of H${}_2$D${}^{+}$, an isotopic variant of triatomic hydrogen, have been used to infer the minimum age of a star-forming molecular cloud [13].

While the hydrogen molecule ${\mathrm{H}}_2$ is by far the most abundant molecule in space, the bulk of ${\mathrm{H}}_2$ molecules in colder environments is difficult to observe with ground-based telescopes. Direct observation of ${\mathrm{H}}_2$ in diffuse and translucent interstellar clouds are obtained either from satellite absorption observations of electronic transitions in the ultraviolet regime (see, e.g. [14,15]) 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 $N_0$ and $N_1$ of the two lowest rotational states of ${\mathrm{H}}_2$ (with J = 0 and J = 1, respectively) and the column densities $N_{(1,1)}$ and $N_{(1,0)}$ of the two lowest states of ${\mathrm{H}}_3^{+}$ (with $(J,G)=(1,1)$ and $(J,G)=(1,0)$, respectively) in the same lines of sight led to a surprise [16]. The temperature derived from the lowest states of ${\mathrm{H}}_3^{+}$, $T_{12}({\mathrm{H}}_3^{+})=32.9/\ln (2N_{(1,1)}/N_{(1,0)})\mathrm{K}$ appeared systematically lower than the temperature of the lowest states of ${\mathrm{H}}_2$, $T_{01}=170.5/\ln (9N_0/N_1)\mathrm{K}$, which is usually found to be in equilibrium with the gas kinetic temperature [15,16]. While the initial studies used very simple model calculations for the ${\mathrm{H}}_3^{+}$ abundances and thermalisation processes, later models, employing large chemical networks, focussed principally on the chemical evolution of the ortho/para forms of ${\mathrm{H}}_2$ and ${\mathrm{H}}_3^{+}$ abundances¹ without considering detailed collisional excitation mechanisms nor introducing thermal balance considerations [17].

Here, we present a different approach and introduce first a master equation describing the evolution of rotational levels of ${\mathrm{H}}_3^{+}$, 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 [18], which solves both the chemical and thermal equilibrium of the cloud and add their contribution to the equilibrium state of ${\mathrm{H}}_2$, 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 ${\mathrm{H}}_3^{+}$ and ${\mathrm{H}}_2$. In Section 3, we describe the most important processes driving the ortho-para ratio of ${\mathrm{H}}_3^{+}$. A master equation for H${}_3^{+}$ state populations is presented in Section 4 together with the corresponding results on the excitation temperature of ${\mathrm{H}}_3^{+}$. 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 ${\mathrm{H}}_3^{+}$. The paper concludes with a brief discussion in Section 6.

2. Nuclear spin and rotational states of ${\mathrm{H}}_3^{+}$ and ${\mathrm{H}}_2$

Both ${\mathrm{H}}_2$ and ${\mathrm{H}}_3^{+}$ exist in two different nuclear spin configurations. For the lowest rotational state of ${\mathrm{H}}_2$, with J = 0, the proton spins are anti-parallel and add up to I = 0. This configuration is denoted as $\mathrm{para}-{\mathrm{H}}_2$ (or $p-{\mathrm{H}}_2$). 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 ${\mathrm{H}}_2$ states with odd rotational quantum numbers have parallel nuclear spin configurations with I = 1, which is denoted as $\mathrm{ortho}-{\mathrm{H}}_2$ (or $o-{\mathrm{H}}_2$). Likewise, all even rotational quantum numbers in H${}_2$ can be attributed to $p-{\mathrm{H}}_2$ with I = 0. Figure 1 shows the three lowest levels of ${\mathrm{H}}_2$ (with $J\le 2$) and their respective energies expressed as $k_{b}T$ (in units of kelvin).

Figure 1. Level scheme of all rotational levels with $J\le 2$ for ${\mathrm{H}}_2$ and with $J\le 3$ for ${\mathrm{H}}_3^{+}$. The level energy (in ${\mathrm{cm}}^{-1}$) is given on the y-axis at the left-hand-side, while the corresponding values in $\mathrm{K}$ are given on the right-hand-side of the graph. Para levels are in blue and ortho levels in red. Black arrows show the possible radiative transitions. Stable and metastable levels are shown with a heavier line. The symmetry-forbidden levels $(0,0)$ and $(2,0)$ are shown by dotted lines for completeness. The origin of the energies is fixed at the lowest permitted level.

The excitation temperature of the two lowest states of ${\mathrm{H}}_2$ is defined as (3) $T_{01}=\frac{\mathrm{\Delta }{E}_{01}/k_b}{\ln (g_1/g_0\cdot N_0/N_1)},$(3) where $\mathrm{\Delta }{E}_{01}/k_b=170.476\mathrm{K}$ stands for the energy difference between the states, $g_1/g_0=9$ is the ratio of the multiplicities of the states with J = 1 and J = 0, respectively, and $N_1$ and $N_0$ denote their populations. An analogous equation can be given for the excitation temperature $T_{02}$ between the states with J = 2 and J = 0, for which $\mathrm{\Delta }{E}_{02}/k_b=509.864$ K and $g_2/g_0=5$, and both temperatures are usually found to be consistent with one another [15], giving a good proxy of the kinetic gas temperature. Typical values for $T_{01}$ in these diffuse sightlines range from 50 to $70\mathrm{K}$ [15,19–21].

The nuclear spin configurations of ${\mathrm{H}}_3^{+}$ are also denoted by para and ortho for $I=1/2$ and $I=3/2$, respectively. As with ${\mathrm{H}}_2$, 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 [22] 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 $I=3/2$ or ortho levels of ${\mathrm{H}}_3^{+}$, while all other levels (effectively fulfilling $G=3n\pm 1$) are of the para configuration, with $I=1/2$. The right-hand side of Figure 1 shows all ${\mathrm{H}}_3^{+}$ rotational levels with $J\le 3$ 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 $(J,G)=(0,0)$ ground state and the $(J,G)=(2,0)$ 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 ${\mathrm{H}}_3^{+}$, we will from now on refer to all ${\mathrm{H}}_3^{+}$ states by giving their $(J,G)$ quantum numbers only, e.g. $(1,1)$ and $(1,0)$ for the lowest para and ortho states, respectively.

With the notable exception of the Galactic centre [12,23,24], only the lowest $(1,1)$ para-state of ${\mathrm{H}}_3^{+}$ and the lowest ortho-state $(1,0)$ have ever been detected in the interstellar medium. Consequently, the ${\mathrm{H}}_3^{+}$ excitation temperature, $T_{12}({\mathrm{H}}_3^{+})$, derived from observations is usually calculated using the column density ratio of these two states (4) $T_{12}({\mathrm{H}}_3^{+})=\frac{\mathrm{\Delta }E/k_b}{\ln ((g_{1,0}/g_{1,1})\cdot N_{(1,1)}/N_{(1,0)})},$(4) where $\mathrm{\Delta }E/k_b=32.86\mathrm{K}$, and $g_{1,0}=12$ and $g_{1,1}=6$ denote the total degeneracies of the ortho- and para-state, respectively. Most observations yield values for $T_{12}({\mathrm{H}}_3^{+})$ between 20 and $40\mathrm{K}$ [16,17], systematically lower than the ${\mathrm{H}}_2$ excitation temperature $T_{01}$. As para- and ortho-states of ${\mathrm{H}}_2$ and ${\mathrm{H}}_3^{+}$ 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 ${\mathrm{H}}_3^{+}$

3.1. Formation of ${\mathrm{H}}_3^{+}$

${\mathrm{H}}_3^{+}$ formation via the ${\mathrm{H}}_2^{+}$ + ${\mathrm{H}}_2$ reaction (Equation   (1(1) ${\mathrm{H}}_2+{\mathrm{H}}_2^{+}⟶{\mathrm{H}}_3^{+}+\mathrm{H}.$(1) )) has been studied for many years by various experimental techniques. A recent compilation of the results can be found in [25]. Particularly noteworthy is a recent study that employs excited Rydberg molecules in a merged beams approach to reach collision energies between 5 and 60$\mathrm{K}$ [26]. The results are in very good agreement with the previous measurements [27] conducted at somewhat higher energies. A recommended fit to the overall cross section can be found in [25], 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) $k_{\mathrm{form}}=2.27{10}^{-9}\times {\left(\frac{T}{300}\right)}^{-0.06}{\mathrm{cm}}^3{\mathrm{s}}^{-1},$(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 $2.08\times {10}^9{\mathrm{cm}}^3{\mathrm{s}}^{-1}$ that is reported in most astrochemistry databases (and based on a previous experimental study [28], which was conducted at room temperature).

To determine the nuclear spin of the ${\mathrm{H}}_3^{+}$ ions created by reaction (1(1) ${\mathrm{H}}_2+{\mathrm{H}}_2^{+}⟶{\mathrm{H}}_3^{+}+\mathrm{H}.$(1) ), we refer to the selection rules outlined by Oka [29]. 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 ${\mathrm{H}}_2^{+}$ and ${\mathrm{H}}_2$. 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 $n(p-{\mathrm{H}}_2)$ and $n(o-{\mathrm{H}}_2)$ of the para- and ortho-species of ${\mathrm{H}}_2$, and $n(p-{\mathrm{H}}_3^{+})$ and $n(o-{\mathrm{H}}_3^{+})$ of ${\mathrm{H}}_3^{+}$, respectively. The para-fraction for both species is then denoted by (6) $p_2=\frac{n(p-{\mathrm{H}}_2)}{n(p-{\mathrm{H}}_2)+n(o-{\mathrm{H}}_2)}$(6) and (7) $p_3=\frac{n(p-{\mathrm{H}}_3^{+})}{n(p-{\mathrm{H}}_3^{+})+n(o-{\mathrm{H}}_3^{+})}.$(7) Following the argumentation by Crabtree et al. [16], we assume that the cosmic ray ionisation of ${\mathrm{H}}_2$ (constituting the main source of ionisation that initiates the ${\mathrm{H}}_3^{+}$ formation) does not affect the nuclear spin of the molecule, and therefore the para-fraction of ${\mathrm{H}}_2^{+}$ has the same value $p_2$.

With these assumptions the $p-{\mathrm{H}}_3^{+}$ fraction at formation, $p_3^f$, can be derived as a function of $p_2$, using the nuclear spin branching fractions given in [29] (for details see Table 4 in [16]), resulting in a simple linear dependence of the form (8) $p_3^f=\frac{1+2p_2}{3}.$(8) The values of $p_2$ and $p_3^f$ are displayed as a function of temperature in Figure 2 for ${\mathrm{H}}_2$ in thermal equilibrium.

Figure 2. Fraction $p_2$ of $p-{\mathrm{H}}_2$ and $p_3^f$ of nascent $p-{\mathrm{H}}_3^{+}$ (Equation (8(8) $p_3^f=\frac{1+2p_2}{3}.$(8) )). The typical range of diffuse cloud kinetic temperature is outlined in yellow.

We further assume that the rotational levels of $p-{\mathrm{H}}_3^{+}$ and $o-{\mathrm{H}}_3^{+}$ are each populated at formation according to their Boltzmann distribution at a temperature corresponding to 2/3 of the reaction exothermicity [30]. 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 ${\mathrm{H}}_3^{+}$

Collisions with electrons, $\mathrm{He}$, $\mathrm{H}$ and ${\mathrm{H}}_2$ contribute to the energy exchange between the rovibrational levels of ${\mathrm{H}}_3^{+}$, but only collisions with ${\mathrm{H}}_2$ and $\mathrm{H}$ may change the nuclear spin of the ${\mathrm{H}}_3^{+}$ ions. To the best of our knowledge, of the above collision partners, only for collisions with ${\mathrm{H}}_2$ and electrons detailed information is available in the literature.

3.2.1. Collisions between ${\mathrm{H}}_3^{+}$ and ${\mathrm{H}}_2$

Before any detailed study of ${\mathrm{H}}_3^{+}$ collision rates with ${\mathrm{H}}_2$ was available, Oka and Epp [31] suggested to use the Langevin expressions to describe the excitation of ${\mathrm{H}}_3^{+}$ detected towards the galacatic centre. However, as a result of the five identical Fermion nuclei involved in these collisions, [29,32,33] pointed out that one should consider the nuclear spin dependence of the total wavefunction. Three different possibilities have to be accounted for (9a) $\begin{array}{rl}& {\mathrm{H}}_3^{+}+{\tilde{\mathrm{H}}}_2⟶{\mathrm{H}}_3^{+}+{\tilde{\mathrm{H}}}_2\mathrm{identity}\\ & \sim \mathrm{inelastic}\mathrm{collision},\end{array}$(9a) (9b) $\begin{array}{rl}& ⟶{\mathrm{H}}_2+(\mathrm{H}{\tilde{\mathrm{H}}}_2{)}^{+}\mathrm{proton}-\mathrm{hop}\\ & \sim \mathrm{reactive}\mathrm{collision},\end{array}$(9b) (9c) $\begin{array}{rl}& ⟶\mathrm{H}\tilde{\mathrm{H}}+(\tilde{\mathrm{H}}{\mathrm{H}}_2{)}^{+}\mathrm{exchange}\\ & \sim \mathrm{reactive}\mathrm{collision}.\end{array}$(9c) Specific selection rules on the nuclear spins of the products can be derived, following [29,34], and they have been used to model hydrogen plasma experiments [16]. The calculations of [32] compared favourably to ion trap measurements of the nuclear spin equilibrium of ${\mathrm{H}}_3^{+}$ and ${\mathrm{H}}_2$ at low temperature [35], 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 ${\mathrm{H}}_3^{+}$ – ${\mathrm{H}}_2$ reaction so far, which was presented by [36]. Their approach is similar to the method of [32], 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 ${\mathrm{H}}_5^{+}$ potential energy surface [37]. We employ a set of rate coefficients that covers collisions of the lowest 24 rotational states of ${\mathrm{H}}_3^{+}$ with ortho-${\mathrm{H}}_2$ and para-${\mathrm{H}}_2$ (in their lowest rotational states) and temperatures up to $500\mathrm{K}$, which were kindly provided by O. Roncero. Those rate coefficients are also currently used in our PDR model calculations [18].

3.2.2. Collisions between ${\mathrm{H}}_3^{+}$, $\mathrm{He}$ and $\mathrm{H}$

$\mathrm{He}$ and $\mathrm{H}$ are additional collision partners that should be taken into account when describing the excitation equilibrium of ${\mathrm{H}}_3^{+}$, but, to the best of our knowledge, no detailed studies are presently available on these systems. Collisions of ${\mathrm{H}}_3^{+}$ with $\mathrm{He}$ can not modify the nuclear spin configuration of ${\mathrm{H}}_3^{+}$, and we approximate the corresponding collision rates by taking the rates for $p-{\mathrm{H}}_2$ and scaling them with the reduced mass factor (while vetoing channels that might change the nuclear spin of ${\mathrm{H}}_3^{+}$), as described previously [38].

Collisions of $\mathrm{H}$ with ${\mathrm{H}}_3^{+}$ may lead to proton exchange or inelastic collisions, similar to the channels discussed for the ${\mathrm{H}}_3^{+}+{\mathrm{H}}_2$ reaction. However, since we are not aware of more detailed information on this process, we will use the collisional rates with $p-{\mathrm{H}}_2$ as a proxy for collisions with $\mathrm{H}$. We checked that this choice has no major influence on our model calculations for the conditions studied here.

3.2.3. Inelastic collisions between ${\mathrm{H}}_3^{+}$ and electrons

Electronic collisions with ${\mathrm{H}}_3^{+}$ preserve the ortho-para character of ${\mathrm{H}}_3^{+}$. They have been computed by [39] 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 ${\mathrm{CH}}^{+}$ allowed for a comparison between experiment and theory for the lowest rotational states [40], which revealed very good agreement.

3.3. Dissociative recombination of ${\mathrm{H}}_3^{+}$

In diffuse gas, the main destruction reaction of ${\mathrm{H}}_3^{+}$ is the dissociative recombination (DR) with electrons. (10a) $\begin{array}{rl}{\mathrm{H}}_3^{+}+e^{-}& ⟶\mathrm{H}+\mathrm{H}+\mathrm{H},\end{array}$(10a) (10b) $\begin{array}{rl}& ⟶{\mathrm{H}}_2+\mathrm{H}.\end{array}$(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 ${\mathrm{H}}_3^{+}$ have been published, with outcomes that differ by orders of magnitude (there are a number of reviews on this topic, see, e.g, [41–43]). In summary, the present consensus is that the absolute rate of the ${\mathrm{H}}_3^{+}$ DR rate coefficient is correctly derived from the storage ring merged beams measurements [43]. Theoretical studies identified the Jahn-Teller effect as a driver for the recombination process at low temperatures [44–46].

Both astrochemistry databases, UMIST ² and KIDA³, report rate coefficients with branching ratios of 2/3 for reaction (10a(10a) $\begin{array}{rl}{\mathrm{H}}_3^{+}+e^{-}& ⟶\mathrm{H}+\mathrm{H}+\mathrm{H},\end{array}$(10a) ) and 1/3 for reaction (10b(10b) $\begin{array}{rl}& ⟶{\mathrm{H}}_2+\mathrm{H}.\end{array}$(10b) ), based on storage ring experiments, and the total DR reaction rate coefficient is given as (11) ${\alpha }_{\mathrm{DR}}^{tot}=6.70\times {10}^{-8}{\left(\frac{T}{300}\right)}^{-0.52}{\mathrm{cm}}^3{\mathrm{s}}^{-1}.$(11) However, there are recent suggestions concerning a possible difference between the DR rate coefficients of the two nuclear spin modifications of ${\mathrm{H}}_3^{+}$. These were first pointed out in the storage ring experiments of [47,48], 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 [49], and updated theoretical studies [46] predict that at low temperature the difference between the rate coefficients for the two nuclear spin modifications may exceed an order of magnitude, with $p-{\mathrm{H}}_3^{+}$ recombining much faster than $o-{\mathrm{H}}_3^{+}$. Two different theoretical values are explicitly reported in [50], supporting plasma studies. We have derived an analytic expression from these results, which we can describe by the formulae (12a) $\begin{array}{rl}{\alpha }_{\mathrm{DR}}^{tot}(p-{\mathrm{H}}_3^{+})& =5.25\times {10}^{-8}(T/300{)}^{-0.75}{\mathrm{cm}}^3{\mathrm{s}}^{-1},\end{array}$(12a) (12b) $\begin{array}{rl}{\alpha }_{\mathrm{DR}}^{tot}(o-{\mathrm{H}}_3^{+})& =6\times {10}^{-8}{\mathrm{cm}}^3{\mathrm{s}}^{-1}T\le 250\mathrm{K}\\ & ={\alpha }_{\mathrm{DR}}(p-{\mathrm{H}}_3^{+})T\ge 250\mathrm{K}.\end{array}$(12b)

Table 1 summarises the values of the principal reaction rate coefficients introduced in the present study, where we have assumed the same branching ratios for reactions (10a(10a) $\begin{array}{rl}{\mathrm{H}}_3^{+}+e^{-}& ⟶\mathrm{H}+\mathrm{H}+\mathrm{H},\end{array}$(10a) ) and (10b(10b) $\begin{array}{rl}& ⟶{\mathrm{H}}_2+\mathrm{H}.\end{array}$(10b) ) as described above.

Table 1. Principal chemical reactions for ${\mathrm{H}}_3^{+}$ formation and destruction.

				A ${}^{\ast }$
Refs	Reaction			(cm${}^3$ s${}^{-1}$)	β	Comment
[25]	${\mathrm{H}}_2^{+}+{\mathrm{H}}_2$	→	${\mathrm{H}}_3^{+}+\mathrm{H}$	2.27 × 10${}^{-9}$	–0.06	computed from x-section
[50]	$p-{\mathrm{H}}_3^{+}+e^{-}$	→	$\mathrm{H}+\mathrm{H}+\mathrm{H}$	3.5 × 10${}^{-8}$	–0.75	present fit
[50]	$p-{\mathrm{H}}_3^{+}+e^{-}$	→	${\mathrm{H}}_2+\mathrm{H}$	1.75 × 10${}^{-8}$	–0.75	present fit
[50]	$o-{\mathrm{H}}_3^{+}+e^{-}$	→	$\mathrm{H}+\mathrm{H}+\mathrm{H}$	4.0 × 10${}^{-8}$	–	present fit, T ≤ 250K
[50]	$o-{\mathrm{H}}_3^{+}+e^{-}$	→	$\mathrm{H}+\mathrm{H}+\mathrm{H}$	3.5 × 10${}^{-8}$	–0.75	present fit, T ≥ 250K
[50]	$o-{\mathrm{H}}_3^{+}+e^{-}$	→	${\mathrm{H}}_2+\mathrm{H}$	2.0 × 10${}^{-8}$	–	present fit, T ≤ 250K
[50]	$o-{\mathrm{H}}_3^{+}+e^{-}$	→	${\mathrm{H}}_2+\mathrm{H}$	1.75 × 10${}^{-8}$	–0.75	present fit, T ≥ 250K

${}^{\ast }$ The reaction rate coefficient, k, is expressed as $A\times (T/300{)}^{\beta }$.

4. Master equation for ${\mathrm{H}}_3^{+}$ state populations

In this section, we describe the master equation for the level populations of ${\mathrm{H}}_3^{+}$, and we show that ${\mathrm{H}}_3^{+}$ 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 $T_{12}({\mathrm{H}}_3^{+})$ 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 $[25:50]\mathrm{K}$ 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-${\mathrm{H}}_3^{+}$ by p-${\mathrm{H}}_2$ 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 ($i\in [1,N]$) of ${\mathrm{H}}_3^{+}$, $n_i$, can be described by (13) $\begin{array}{rl}\frac{\mathrm{d}{n}_i}{\mathrm{d}t}& =k_{\mathrm{form},i}n\left({\mathrm{H}}_2^{+}\right)n\left({\mathrm{H}}_2\right)\\ & -{\alpha }_{\mathrm{DR},i}{n}_{i}n\left(e^{-}\right)& (\mathrm{formation}\mathrm{and}\mathrm{destruction})\\ & +\sum _{j\ne i}{k}_{ji}^{p}n\left({\mathrm{H}}_2\right)p_2{n}_j\\ & -\sum _{j\ne i}{k}_{ij}^{p}n\left({\mathrm{H}}_2\right)p_2{n}_i& (\mathrm{collisions}\mathrm{with}p-{\mathrm{H}}_2)\\ & +\sum _{j\ne i}{k}_{ji}^{o}n\left({\mathrm{H}}_2\right)(1-p_2)n_j\\ & -\sum _{j\ne i}{k}_{ij}^{o}n\left({\mathrm{H}}_2\right)(1-p_2)n_i& (\mathrm{collisions}\mathrm{with}o{\mathrm{H}}_2)\\ & +\sum _{j\ne i,X}{k}_{ji}^{X}n\left(\mathrm{X}\right)n_j\\ & -\sum _{j\ne i,X}{k}_{ij}^{X}n\left(\mathrm{X}\right)n_i& (\mathrm{collisions}\mathrm{with}\mathrm{other}\mathrm{species}X)\\ & +\sum _{i<j}{A}_{ji}{n}_j-\sum _{i>j}{A}_{ij}{n}_i& (\mathrm{radiative}\mathrm{transitions})\end{array}$(13) with $k_{\mathrm{form},i}$ and ${\alpha }_{\mathrm{DR},i}$ denoting the state-dependent chemical formation and destruction rates, $k_{ij}^{p,o}$ the collisional excitation/de-excitation rates with $p-{\mathrm{H}}_2$ and $o-{\mathrm{H}}_2$, and $k_{ij}^X$ denote collision rates where X stands for $\mathrm{H}$, $\mathrm{He}$, $e^{-}$. $A_{ij}$ 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 ${\mathrm{H}}_3^{+}$ is a highly exothermic reaction, which preferentially populates high energy levels. Thus, chemical formation can be regarded as an excitation mechanism.

Introducing $x_i$, the relative populations of ${\mathrm{H}}_3^{+}$, so that $n_i=x_{i}n({\mathrm{H}}_3^{+})$, the total destruction rate is ${\alpha }_{\mathrm{DR}}=\sum _i{x}_i{\alpha }_{\mathrm{DR},i}$. The total formation rate is $k_{\mathrm{form}}=\sum _i{k}_{\mathrm{form},i}$. Since the relative populations $x_i$ are unknown until the differential equations are solved, it is not possible to compute the destruction rate beforehand, as long as the ${\alpha }_{\mathrm{DR},i}$ differ for individual states, and thus $n({\mathrm{H}}_3^{+})$ cannot be computed independently. At steady state,$\frac{d{n}_i}{dt}=0$ and the system of equations to solve is (14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(14) with $R=\frac{n({\mathrm{H}}_2^{+})n({\mathrm{H}}_2)}{n({\mathrm{H}}_3^{+})}$. The value of R is initially unknown, as the total destruction rate of ${\mathrm{H}}_3^{+}$ requires the knowledge of the relative level populations of the molecular ion. So, we explicitly add the conservation equation (15) $\sum _i{x}_i=1.$(15) We get a system of N + 1 equations with N + 1 unknowns, which is easily solved if the densities of ${\mathrm{H}}_2$, ${\mathrm{H}}_2^{+}$, $e^{-}$ 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 ${\mathrm{H}}_2^{+}$ density and electronic fraction

It is possible to reduce the set of equations further by estimating $n({\mathrm{H}}_2^{+})$ and $n(e^{-})$. The main formation and destruction reactions of ${\mathrm{H}}_2^{+}$ are: (16a) $\begin{array}{rl}& {\mathrm{H}}_2+\mathrm{CRP}⟶{\mathrm{H}}_2^{+}+e^{-}\zeta ({\mathrm{s}}^{-1}),\end{array}$(16a) (16b) $\begin{array}{rl}& {\mathrm{H}}_2^{+}+{\mathrm{H}}_2⟶{\mathrm{H}}_3^{+}+\mathrm{H}{k}_{\mathrm{form}}({\mathrm{cm}}^3{\mathrm{s}}^{-1}),\end{array}$(16b) (16c) $\begin{array}{rl}& {\mathrm{H}}_2^{+}+e^{-}⟶\mathrm{H}+\mathrm{H}{\alpha }^{{\mathrm{H}}_2^{+}}({\mathrm{cm}}^3{\mathrm{s}}^{-1}),\end{array}$(16c) where $\mathrm{CRP}$ represents cosmic ray particles and ζ the corresponding ionisation rate of ${\mathrm{H}}_2$. At steady state, this leads to (with ${\alpha }^{{\mathrm{H}}_2^{+}}$ the dissociative recombination rate of ${\mathrm{H}}_2^{+}$) (17) $n\left({\mathrm{H}}_2^{+}\right)=\frac{\zeta n\left({\mathrm{H}}_2\right)}{k_{\mathrm{form}}n\left({\mathrm{H}}_2\right)+{\alpha }^{{\mathrm{H}}_2^{+}}n\left(e^{-}\right)}.$(17) This removes one parameter from the system.

As for $n(e^{-})$, in diffuse gas it is often assumed to be equal to the density of ${\mathrm{C}}^{+}$. However, for high cosmic ray ionisation rates, as, e.g. in the Central Molecular Zone (CMZ) of our galaxy, Le Petit et al. [11] have shown that protons may contribute significantly to the ionisation fraction. Here we compute the electronic density considering both ${\mathrm{H}}^{+}$ and ${\mathrm{He}}^{+}$ 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 ⁴. With these derivations of $n({\mathrm{H}}_2^{+})$ and $n(e^{-})$, the system formed by Equations (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(14) ) and  (15(15) $\sum _i{x}_i=1.$(15) ) depends on five astrophysical parameters: $n_{\mathrm{H}}$⁵, T, $f_m$, $G_0$, and ζ, where we introduced $f_m$ the molecular fraction, $f_m=2n({\mathrm{H}}_2)/n_{\mathrm{H}}$.

4.2. Numerical solution of the coupled equations

We developed ExcitH3p, a FORTRAN code that solves the ${\mathrm{H}}_3^{+}$ coupled system of equations (Equation  (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(14) )) for any number of levels, and then computes the two-level excitation temperature (Equation  (4(4) $T_{12}({\mathrm{H}}_3^{+})=\frac{\mathrm{\Delta }E/k_b}{\ln ((g_{1,0}/g_{1,1})\cdot N_{(1,1)}/N_{(1,0)})},$(4) )) (18) $T_{12}\left({\mathrm{H}}_3^{+}\right)=32.86\mathrm{K}/\ln \left(\frac{2x_1}{x_2}\right).$(18) Here $x_1$ and $x_2$ denote the relative populations of the $(1,1)$ ground state and the first excited state $(1,0)$, respectively.

In practice, the 24 energetically lowest rotational levels of ${\mathrm{H}}_3^{+}$ 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 $(7,6)$ ortho level, located $2190\mathrm{K}$ above the ground state. We find that the solution to the set of equations (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(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 3 for a range of typical values of the cosmic ionisation rate ζ and gas kinetic temperature T. The total gas density is set to $n_{\mathrm{H}}={10}^2{\mathrm{cm}}^{-2}$ and the molecular fraction to $f_m=0.8$. We find that in the entire parameter space, $T_{12}({\mathrm{H}}_3^{+})$, the excitation temperature given by the two lowest ${\mathrm{H}}_3^{+}$ levels, is systematically lower than the kinetic temperature. The overall magnitude of the discrepancy between the two temperatures is in agreement with the observations [10,17]. This outcome is a sole property of the microscopic excitation and de-excitation mechanisms of the individual quantum levels of ${\mathrm{H}}_3^{+}$, which are detailed above. Figure 3 also shows that a maximum of about $44\mathrm{K}$ is reached for $T_{12}({\mathrm{H}}_3^{+})$ at kinetic temperatures around $100\mathrm{K}$. 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-${\mathrm{H}}_3^{+}$ are the same: $T_{12}({\mathrm{H}}_3^{+}{)}^{max}$ = 38 K for T${}_{gas}$ = 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 ${\mathrm{H}}_3^{+}$ levels included in the model.

Figure 3. The excitation temperature $T_{12}({\mathrm{H}}_3^{+})$ as calculated using ExcitH3p, as a function of ζ and T for a density $n_{\mathrm{H}}={10}^2{\mathrm{cm}}^{-3}$ and a molecular fraction $f_m=0.8$. Note that the lowest 24 levels of ${\mathrm{H}}_3^{+}$ are included in the model, while only the lowest two states are used to determine $T_{12}({\mathrm{H}}_3^{+})$.

Figure 4 shows the influence of density, showing that $T_{12}({\mathrm{H}}_3^{+})$ is not sensitive to this parameter up to densities of more than ${10}^3{\mathrm{cm}}^{-3}$ if the gas temperature is below ∼ 80K.

Figure 4. The excitation temperature $T_{12}({\mathrm{H}}_3^{+})$ as a function of the gas kinetic temperature T and density $n_{\mathrm{H}}$ for $\zeta ={10}^{-16}{\mathrm{s}}^{-1}$, a molecular fraction $f_m=0.8$ and considering the first 24 levels of ${\mathrm{H}}_3^{+}$.

4.2.1. Influence of the number of ${\mathrm{H}}_3^{+}$ levels included in the model

Interesting conclusions can be drawn from the examination of the dependence of our results on the number N of ${\mathrm{H}}_3^{+}$ levels included in the model. Except for the Galactic centre, only the two lowest rotational levels of H${}_3^{+}$ have been detected in the ISM so far. Consequently, the interpretation of astrophysical observations is usually restricted to these two levels. Figure 5 displays the computed $T_{12}({\mathrm{H}}_3^{+})$ value as a function of the kinetic temperature T, for different values of N. We perform this comparison for the following parameters: $\zeta =5{10}^{-16}{\mathrm{s}}^{-1}$, $n_{\mathrm{H}}=100{\mathrm{cm}}^{-3}$, and $f_m=0.8$. 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, $T_{12}({\mathrm{H}}_3^{+})$ 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 $(3,3)$ level, representing a possible sink for ${\mathrm{H}}_3^{+}$ population.

Figure 5. $T_{12}({\mathrm{H}}_3^{+})$ as a function of the number of levels taken into account in the coupled set of equations (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(14) ). We assume values of $n_{\mathrm{H}}={10}^2{\mathrm{cm}}^{-3}$, $\zeta =5{10}^{-16}{\mathrm{s}}^{-1}$, $f_m=0.8$. The straight line in grey is the first diagonal showing complete thermalisation.

The decrease of $T_{12}({\mathrm{H}}_3^{+})$ with increasing gas kinetic temperature above $100\mathrm{K}$ seems puzzling at first. However, it can be understood by considering the various possible decay paths from the two excited para levels $(2,2)$ (corresponding to relative population $x_3$ in our enumeration) and $(2,1)$ (corresponding to $x_4$), and the lowest excited ortho level $(3,3)$ (corresponding to $x_5$). It is important to note that the $(3,3)$ level is metastable with an infinite radiative lifetime. Hence it can only be depopulated in collisional processes. Examination of the collision rates with ${\mathrm{H}}_2$ shows that its branching ratios towards ortho and para are approximately equal in the 30–300$\mathrm{K}$ temperature range. Both para levels $(2,2)$ and $(2,1)$, 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 $(3,3)$ 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 $(3,3)$ and inhibits population of the lowest ortho state $(1,0)$.

4.3. Two-level approximation

In a restricted range of kinetic temperatures around $T\simeq 25-50\mathrm{K}$, the two-level case is a fair approximation to the full system, as seen in Figure 5. It allows to considerably simplify the system and to derive an analytic expression for $x_2/x_1$ that offers the opportunity to highlight the key microscopic mechanisms at work.

We restrict our study to molecular hydrogen collisions and introduce $k_{12}=k_{12}^o(1-p_2)+k_{12}^p{p}_2$, the total collisional excitation rate due to ${\mathrm{H}}_2$ from level 1 to level 2 (and the equivalent for $k_{21}$). Then, the coupled equations system (Equation  (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(14) )) reduces to: (19a) $\begin{array}{rl}& (n\left({\mathrm{H}}_2\right)k_{12}+{\alpha }_{\mathrm{DR},1}n\left(e^{-}\right))x_1-n\left({\mathrm{H}}_2\right)k_{21}{x}_2\\ & =k_{\mathrm{form},1}R\end{array}$(19a) (19b) $\begin{array}{rl}& -n\left({\mathrm{H}}_2\right)k_{12}{x}_1+(n\left({\mathrm{H}}_2\right)k_{21}+{\alpha }_{\mathrm{DR},2}n\left(e^{-}\right))x_2\\ & =k_{\mathrm{form},2}R\end{array}$(19b)

Introducing $k_{\mathrm{form}}=k_{\mathrm{form},1}+k_{\mathrm{form},2}$, we can compute the $x_2/x_1$ ratio: (20) $\frac{x_2}{x_1}=\frac{n\left({\mathrm{H}}_2\right)k_{12}{k}_{\mathrm{form}}+{\alpha }_{\mathrm{DR},1}n\left(e^{-}\right)k_{\mathrm{form},2}}{n\left({\mathrm{H}}_2\right)k_{21}{k}_{\mathrm{form}}+{\alpha }_{\mathrm{DR},2}n\left(e^{-}\right)k_{\mathrm{form},1}},$(20) The factor R cancels out and the ratio $x_2/x_1$ does not depend on $n({\mathrm{H}}_3^{+})$. We introduce the electronic fraction $x_e=n({\mathrm{e}}^{-})/n_{\mathrm{H}}$ and the configuration-specific formation rates of H${}_3^{+}$ with $k_{\mathrm{form},1}=p_3^f{k}_{\mathrm{form}}$ and $k_{\mathrm{form},2}=(1-p_3^f)k_{\mathrm{form}}$. Furthermore, we apply detailed balance to the H${}_3^{+}$–H${}_2$ collisional rates, yielding $k_{12}=\frac{g_2}{g_1}\exp (-\frac{E_{21}}{T})k_{21}$. Finally, we get the following expression (21) $\begin{array}{rl}\frac{x_2}{x_1}& =\frac{g_2}{g_1}\exp (-\frac{E_{21}}{T})\\ & \times \left[\frac{1+(1-p_2)\frac{4}{3}\frac{{\alpha }_{\mathrm{DR},1}}{k_{21}}\frac{x_e}{f_m}\frac{g_1}{g_2}\exp \left(\frac{E_{21}}{T}\right)}{1+(1+2p_2)\frac{2}{3}\frac{{\alpha }_{\mathrm{DR},2}}{k_{21}}\frac{x_e}{f_m}}\right]\\ & =\frac{g_2}{g_1}\exp (-\frac{E_{21}}{T_{12}\left({\mathrm{H}}_3^{+}\right)}).\end{array}$(21) Apart from the kinetic temperature T, the $x_2/x_1$ ratio depends on the collisional and dissociative recombination rate coefficients, the molecular fraction $f_m$ of ${\mathrm{H}}_2$, and the electronic fraction $x_e$. It is independent of the cosmic ray ionisation rate ζ and of the total formation reaction rate coefficient of H${}_3^{+}$, but – crucially – not of the branching ratios, which effectively enter the equation through the H${}_2$ para fraction $p_2$.

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$\mathrm{K}$) the value of $p_2$ ranges from 0.6 to 1 (see Figure 2), and numerical analysis reveals that the $(1-p_2)$ and $(1+2p_2)$ 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 $\frac{x_2}{x_1}$ ratio compared to the thermal value (in agreement with astronomical observations). The extent of the deviation from thermal equilibrium depends critically on the ratios $\frac{{\alpha }_{\mathrm{DR},1}}{k_{21}}$ and $\frac{{\alpha }_{\mathrm{DR},2}}{k_{21}}$ 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 $k_{21}$ – 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 ${\alpha }_{\mathrm{DR},1}$ and ${\alpha }_{\mathrm{DR},2}$.

The picture that emerges is that the excitation temperature of ${\mathrm{H}}_3^{+}$ appears lower than the nominal gas temperature because of incomplete thermalisation. The formation process strongly favours the formation of $p-{\mathrm{H}}_3^{+}$ at low temperatures, as the para-fraction $p_2$ of ${\mathrm{H}}_2$ is large. The electron recombination process then removes H${}_3^{+}$ before it can reach thermal equilibrium in collisions with ${\mathrm{H}}_2$. 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 ${\mathrm{H}}_3^{+}-{\mathrm{H}}_2$ collisional rate coefficients (or sufficiently reduced electron recombination rates) the calculations reach complete thermal equilibrium between the excitation temperature $T_{12}({\mathrm{H}}_3^{+})$ and the kinetic temperature T.

In essence, the nuclear spin restrictions of the ${\mathrm{H}}_3^{+}$ formation reaction produce an over-proportional amount of ${\mathrm{H}}_3^{+}$ in the para configuration, and thermalisation in collisions with H${}_2$ is too slow to reach equilibrium before the ions are destroyed by free electrons. This is to be contrasted to the situation for ${\mathrm{H}}_2$, 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  (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(14) )) with observations. Table 2 presents the ${\mathrm{H}}_3^{+}$ excitation temperature reported in the literature for a few local diffuse clouds as well as the corresponding ${\mathrm{H}}_2$ excitation temperatures. Data for ${\mathrm{H}}_2$ come from Copernicus [19] and FUSE [20,21] satellite observations. The proton densities reported in this table are those published in the papers reporting the ${\mathrm{H}}_2$ data. Determination of diffuse cloud density with observations of $\mathrm{H}$ and ${\mathrm{H}}_2$ 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 ${\mathrm{H}}_2$ belongs. So, these densities are to be seen as order of magnitude estimates, and should be considered only as indicative.

Table 2. Selection of sightlines with both H${}_2$ and H${}_3^{+}$ observations.

		$n_{\mathrm{H}}$	$T_{01}^{obs}({\mathrm{H}}_2)$	$T_{12}^{obs}({\mathrm{H}}_3^{+})$	$N({\mathrm{H}}_3^{+})$	$N({\mathrm{H}}_2)$
		$({\mathrm{cm}}^{-3})$	$(\mathrm{K})$	$(\mathrm{K})$	$({10}^{13}{\mathrm{cm}}^{-2})$	$({10}^{20}{\mathrm{cm}}^{-2})$	References
HD154368		240	$51\pm 8$	$20\pm 4$	9.37	14.4	[10,16,20]
HD73882		520	$51\pm 6$	$23\pm 3$	9.02	12.9	[10,16,20]
HD27778	62 Tau	280	$55\pm 7$	$29\pm 4$	6.49	6.23	[17,20]
HD24398	ζ Per	215	$57\pm 6$	$28\pm 4$	6.26	4.75	[10,16,19]
HD24534	X Per	325	$57\pm 4$	${46}_{-13}^{+21}$	7.34	8.38	[10,16,20]
HD41117	${\chi }^2$ Ori	200	$60\pm 7$	$29\pm 13$	5.29	4.90	[17,21]
HD110432		140	$68\pm 5$	$30\pm 2$	5.22	4.37	[10,16]
HD210839	λ Cep	115	$72\pm 6$	$34\pm 2$	7.58	6.88	[10,20]
HD43384	9 Gem	120${}^{\ast }$	$74\pm 15$	$38\pm 11$	4.07	7.36	[17,21]

${}^{\ast }$ Value unavailable. Educated guess only.

Using the ExcitH3p program that solves Equation (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(14) ), we compute $T_{12}({\mathrm{H}}_3^{+})$ for all the lines of sight. We assume a cosmic ray ionisation rate of $\zeta ={10}^{-16}{\mathrm{s}}^{-1}$ and a molecular fraction of $f_m=0.8$. For the gas density, the values in Table 2 are used, and for the gas temperature we use the values derived from ${\mathrm{H}}_2$ observations $T_{01}^{obs}({\mathrm{H}}_2)$.

5.1.1. Model results and sensitivity to the DR rate coefficients

Figure 6 shows the results for three different hypothesis for ${\alpha }_{DR}(o-{\mathrm{H}}_3^{+})$, and the value provided by Equation (12a(12a) $\begin{array}{rl}{\alpha }_{\mathrm{DR}}^{tot}(p-{\mathrm{H}}_3^{+})& =5.25\times {10}^{-8}(T/300{)}^{-0.75}{\mathrm{cm}}^3{\mathrm{s}}^{-1},\end{array}$(12a) ) for ${\alpha }_{DR}(p-{\mathrm{H}}_3^{+})$. We see a quasi-linear variation of $T_{12}({\mathrm{H}}_3^{+})$ with $T_{01}$, which is well reproduced by the models. Using ${\alpha }_{DR}$ from Equation (12b(12b) $\begin{array}{rl}{\alpha }_{\mathrm{DR}}^{tot}(o-{\mathrm{H}}_3^{+})& =6\times {10}^{-8}{\mathrm{cm}}^3{\mathrm{s}}^{-1}T\le 250\mathrm{K}\\ & ={\alpha }_{\mathrm{DR}}(p-{\mathrm{H}}_3^{+})T\ge 250\mathrm{K}.\end{array}$(12b) ) (red circles) leads to temperatures which are too high, while using the same rate for $o-{\mathrm{H}}_3^{+}$ and $p-{\mathrm{H}}_3^{+}$ (pink circles) leads to temperatures which are too low. An empirical adjustment using ${\alpha }_{\mathrm{DR}}(p-{\mathrm{H}}_3^{+})$ from Equation (12a(12a) $\begin{array}{rl}{\alpha }_{\mathrm{DR}}^{tot}(p-{\mathrm{H}}_3^{+})& =5.25\times {10}^{-8}(T/300{)}^{-0.75}{\mathrm{cm}}^3{\mathrm{s}}^{-1},\end{array}$(12a) ) and ${\alpha }_{\mathrm{DR}}(o-{\mathrm{H}}_3^{+})={\alpha }_{\mathrm{DR}}(p-{\mathrm{H}}_3^{+})/1.5$ (green circles) gives a very satisfying result, given that no attempt has been made to optimise other parameters.

Figure 6. Computed H${}_3^{+}$ excitation temperature $T_{12}({\mathrm{H}}_3^{+})$ as a function of the observed ${\mathrm{H}}_2$ excitation temperature for the 9 lines of sight given in Table 2. Sightlines with the same $T_{01}$ have been shifted by $0.1\mathrm{K}$ for clarity. Three different hypothesis are used for the dissociative recombination rate of $o-{\mathrm{H}}_3^{+}$ with electrons (see text).

There is a single noticeable exception: $\mathrm{X}$ Per. However, $T_{12}({\mathrm{H}}_3^{+})$ of this line of sight suffers from a particularly large error bar as listed in Table 2. A high excitation temperature of ${\mathrm{H}}_3^{+}$ can only be reached for a kinetic temperatures close to $100\mathrm{K}$ with our models, as can be seen in Figure 3.

We note that typically 10 to $12\mathrm{\%}$ of ${\mathrm{H}}_3^{+}$ ions are in excited stable or metastable levels above the respective lowest ortho or para levels, such as $(3,3)$, $(4,4)$, $(5,5)$, $(6,6)$ and $(7,7)$.

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 ${\mathrm{H}}_2$ and ${\mathrm{H}}_3^{+}$ excitation.

We use our previous study dedicated to ζ Per [52] as a starting point, but we also account for the many updates made since then [11,18,53]. For the present study, we have introduced in the Meudon PDR code the ortho/para dependence of the formation reaction of ${\mathrm{H}}_3^{+}$ via ${\mathrm{H}}_2$ + ${\mathrm{H}}_2^{+}$ collisions and the nuclear spin dependence of the dissociative recombination rate coefficient, in addition to the other excitation mechanisms of ${\mathrm{H}}_3^{+}$. The collisional excitation of ${\mathrm{H}}_2$ by ${\mathrm{H}}^{+}$, that has been revisited by [54] for highly rovibrationally excited levels, has also been updated. Observations suggest a total visual extinction $A_{\mathrm{V}}=0.9\mathrm{mag}$ ($R_{\mathrm{V}}=2.8$, $E_{B-V}=0.32$, $N_{\mathrm{H}}/E_{B-V}=5.2{10}^{21}{\mathrm{cm}}^{-2}$). In our models the cloud is illuminated from both sides by the standard ISRF ($G_0=1$).

We consider two different scenarios

1. Constant density and constant temperature,

2. 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 $A_{\mathrm{V}}=0.7$ for both models A and B, allowing to account for the molecular hydrogen abundance.

The results are summarised in Table 3. The examples shown here have been selected from an evaluation of the following ${\chi }^2$, where σ are the observational uncertainties (22) $\begin{array}{rl}{\chi }^2& =\frac{1}{4}(\frac{(N_{Obs}({\mathrm{H}}_2)-N_{Mod}({\mathrm{H}}_2){)}^2}{{\sigma }_{N({\mathrm{H}}_2)}^2}\\ & +\frac{(N_{Obs}({\mathrm{H}}_3^{+})-N_{Mod}({\mathrm{H}}_3^{+}){)}^2}{{\sigma }_{N({\mathrm{H}}_3^{+})}^2}\\ & +\frac{(T_{01,Obs}({\mathrm{H}}_2)-T_{01,Mod}({\mathrm{H}}_2){)}^2}{{\sigma }_{T_{01}({\mathrm{H}}_2)}^2}\\ & +\frac{(T_{12,Obs}({\mathrm{H}}_3^{+})-T_{12,Mod}({\mathrm{H}}_3^{+}){)}^2}{{\sigma }_{T_{12}({\mathrm{H}}_3^{+})}^2}).\end{array}$(22)

Table 3. PDR model results. Column densities are in ${\mathrm{cm}}^{-2}$ and excitation temperatures in $\mathrm{K}$.

Obs.	$N({\mathrm{H}}_2)$	$T_{01}$	$N({\mathrm{H}}_3^{+})$	$T_{12}({\mathrm{H}}_3^{+})$	${\chi }^2$
ζ Per	4.75(20)	$57\pm 6$	$6.26(13)$	28 ± 4
[10]	$\pm 0.95(20)$		$\pm 0.52(13)$
Models	$n_{\mathrm{H}}=150{\mathrm{cm}}^{-3}$, ($\zeta =8{10}^{-16}{\mathrm{s}}^{-1}$
A	$5.00(20)$	56.9	$5.62(13)$	31.0	0.54
B	$5.00(20)$	62.3	$6.30(13)$	33.7	0.73
	Same parameters, but ${\alpha }_{DR}^o={\alpha }_{DR}^p/1.5$
A	$5.00(20)$	56.9	$4.58(13)$	23.0	3.0
B	$5.00(20)$	62.3	$5.18(13)$	26.0	1.4

Note: Numbers in parenthesis are power of 10.

The proposed ionisation rates are rather high, but in line with the other recent evaluations based on ${\mathrm{H}}_3^{+}$ and ${\mathrm{OH}}^{+}$abundances [9,10,55]. 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 $\mathrm{OH}$ and $\mathrm{HD}$ column densities. The constraints imposed by the excitation temperature of ${\mathrm{H}}_2$ and ${\mathrm{H}}_3^{+}$ are not sensitive to the cosmic ionisation rate of ${\mathrm{H}}_2$, as shown in Figure 3. We then recover the predictions based on molecular ion observations.

Using the recombination rate from Equation (12b(12b) $\begin{array}{rl}{\alpha }_{\mathrm{DR}}^{tot}(o-{\mathrm{H}}_3^{+})& =6\times {10}^{-8}{\mathrm{cm}}^3{\mathrm{s}}^{-1}T\le 250\mathrm{K}\\ & ={\alpha }_{\mathrm{DR}}(p-{\mathrm{H}}_3^{+})T\ge 250\mathrm{K}.\end{array}$(12b) ), the excitation temperature of ${\mathrm{H}}_3^{+}$ 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 ${\mathrm{H}}_2$. The total amount of ${\mathrm{H}}_3^{+}$ is lower due to the overall larger recombination rate. The resulting ${\chi }^2$ varies accordingly. This illustrates the very high sensitivity of ${\mathrm{H}}_3^{+}$ 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 ${\mathrm{H}}_3^{+}$ and ${\mathrm{H}}_2$ 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 ${\mathrm{H}}_3^{+}$ in the diffuse interstellar medium. We provide references to the best data available to date.

We show that a $0D$ statistical model of ${\mathrm{H}}_3^{+}$ 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 ${\mathrm{H}}_3^{+}$. Specifically, the formation of $p-{\mathrm{H}}_3^{+}$ by $p-{\mathrm{H}}_2$ at temperatures below $70\mathrm{K}$ is efficient. Considering the individual levels of ${\mathrm{H}}_3^{+}$, we find that reactive collisions with ${\mathrm{H}}_2$ 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 ${\mathrm{H}}_3^{+}$ 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-H₃⁺ 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 H₂ – 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 H₃⁺ 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 $T_{12}({\mathrm{H}}_3^{+})$, 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 $10\mathrm{\%}$ of ${\mathrm{H}}_3^{+}$ 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 ${\mathrm{H}}_3^{+}$ excitation will not be possible as long as we still lack accurate state specific recombination rates with electrons at temperatures between 20 and $200\mathrm{K}$. Numerical manipulations of the rate coefficients shown in Section 5 reveal the sensitivity of $T({\mathrm{H}}_3^{+})$ to the ratio ${\alpha }_{DR}^o/{\alpha }_{DR}^p$. 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 ${\mathrm{H}}_3^{+}$ with $\mathrm{He}$ and $\mathrm{H}$. In particular, exchange of hydrogen atoms during collisions with $\mathrm{H}$ may impact the ortho to para ratio of ${\mathrm{H}}_3^{+}$. Here we used estimated values for the rate coefficients scaled from reaction rates with $p-{\mathrm{H}}_2$ 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 ${\mathrm{H}}_3^{+}$ 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

ER, JLB & FLP were supported in part by the Programme National ‘Physique et Chimie du Milieu Interstellaire’ (PCMI) of Centre National de la Recherche Scientifique (CNRS)/INSU with INC/INP, co-funded by CEA and CNES. FK, AO & HK acknowledge financial support by the Max Planck Society (Max-Planck-Gesellschaft).

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 $G_0$ is defined in [51]. It controls the interstellar grain charge, which impacts the recombination of ${\mathrm{H}}^{+}$ and ${\mathrm{He}}^{+}$.

5 $n_{\mathrm{H}}$, the proton density is expressed as = $n(\mathrm{H})+2n({\mathrm{H}}_2)+n({\mathrm{H}}^{+})$.

Appendices

Appendix 1. Two-level approximation for ${\mathrm{H}}_2$

As first recognised by [56,57] from Copernicus observations, $T_{01}$ of ${\mathrm{H}}_2$ 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 [58], and where collisions with ${\mathrm{H}}_3^{+}$ 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 ${\mathrm{H}}^{+}$ with rates $k_{01}$ and $k_{10}$ (reactive collisions with other species ($\mathrm{H}$, ${\mathrm{H}}_3^{+}$) are negligible here). Besides collisions, these two levels are populated by direct formation on grains with rates $k_{f0}$ and $k_{f1}$ or depopulated by photodissociation with rates $d_0$ and $d_1$. Other chemical reactions have only a minor impact. The resulting balance equations are $\begin{array}{rl}& (k_{01}n\left({\mathrm{H}}^{+}\right)+d_0)x_0-k_{10}n\left({\mathrm{H}}^{+}\right)x_1=k_{f0}n\left(\mathrm{H}\right)\frac{n_{\mathrm{H}}}{n\left({\mathrm{H}}_2\right)},\\ & (k_{10}n\left({\mathrm{H}}^{+}\right)+d_1)x_1-k_{01}n\left({\mathrm{H}}^{+}\right)x_0=k_{f1}n\left(\mathrm{H}\right)\frac{n_{\mathrm{H}}}{n\left({\mathrm{H}}_2\right)}.\end{array}$This system can be solved for $x_0$ and $x_1$ and leads to $\frac{x_1}{x_0}=\frac{k_{01}n\left({\mathrm{H}}^{+}\right)k_{f0}+(k_{01}n\left({\mathrm{H}}^{+}\right)+d_0)k_{f1}}{(k_{10}n\left({\mathrm{H}}^{+}\right)+d_1)k_{f0}+k_{10}n\left({\mathrm{H}}^{+}\right)k_{f1}}.$In the temperature range from 50 to $100\mathrm{K}$ appropriate do diffuse and translucent cloud conditions, the Boltzmann factor $\frac{g_1}{g_0}\exp (-\frac{E_{10}}{T})=9\exp (-\frac{170.5}{T})$ varies from 0.5 to 1.6. So, $k_{01}$ and $k_{10}$ remain close to one another. Furthermore, the formation rates $k_{f0}$ and $k_{f1}$ are close to one another and simplify. The ratio can be arranged using detailed balance as: $\frac{x_1}{x_0}=\frac{g_1}{g_0}\exp (-\frac{E_{10}}{T})\frac{1+\frac{d_0}{2k_{01}n\left({\mathrm{H}}^{+}\right)}}{1+\frac{d_1}{2k_{10}n\left({\mathrm{H}}^{+}\right)}}$In regions of low radiation field where most diffuse clouds are found the $\mathrm{H}/{\mathrm{H}}_2$ transition is very close to the edge of the cloud, as shown in [59]. Hence, in most of the region that builds ${\mathrm{H}}_2$ column density the dissociation rates $d_0$ and $d_1$ are about 4 orders of magnitude lower than the products $k_{10}n({\mathrm{H}}^{+})$ and $k_{01}n({\mathrm{H}}^{+})$. Thus, the correction factor coming from the chemistry is very close to 1 and the ratio $\frac{x_1}{x_0}$ 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 $13.6\mathrm{eV}$ are ionised, providing a minimal electronic abundance of ${\delta }_M{n}_{\mathrm{H}}$, where ${\delta }_M$ is the fraction of relevant metals (mostly $\mathrm{C}$ and $\mathrm{S}$, with traces of $\mathrm{Si}$ and other heavier species). In the following, we take ${\delta }_M=1.55{10}^{-4}$.

Additional electrons come from ${\mathrm{H}}^{+}$ and ${\mathrm{He}}^{+}$ 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 [60], Section 13.6, extended to include $\mathrm{He}$. The balance equations are $\begin{array}{rl}& {\zeta }_{\mathrm{H}}n\left(\mathrm{H}\right)={\alpha }_{rr}\left({\mathrm{H}}^{+}\right)n\left({\mathrm{H}}^{+}\right)n\left(e^{-}\right)+{\alpha }_{gr}\left({\mathrm{H}}^{+}\right)n_{\mathrm{H}}n\left({\mathrm{H}}^{+}\right),\\ & {\zeta }_{\mathrm{He}}n\left(\mathrm{He}\right)={\alpha }_{rr}\left({\mathrm{He}}^{+}\right)n\left({\mathrm{He}}^{+}\right)n\left(e^{-}\right)\\ & +{\alpha }_{gr}\left({\mathrm{He}}^{+}\right)n_{\mathrm{H}}n\left({\mathrm{He}}^{+}\right).\end{array}$Where ${\alpha }_{rr}$ is the radiative recombination rate, and ${\alpha }_{gr}$ the rate of recombination on grains. ${\zeta }_{\mathrm{H}}$ and ${\zeta }_{\mathrm{He}}$ are the cosmic ray ionisation rates of $\mathrm{H}$ and $\mathrm{He}$, respectively. With respect to ${\mathrm{H}}_2$, we use ${\zeta }_{\mathrm{H}}=0.77\zeta$, including secondary ionisation, and ${\zeta }_{\mathrm{He}}=0.5\zeta$. The total abundance of electrons is $n\left(e^{-}\right)={\delta }_{\mathrm{M}}{n}_{\mathrm{H}}+n\left({\mathrm{H}}^{+}\right)+n\left({\mathrm{He}}^{+}\right).$For all abundances, we write $x(\mathrm{X})=n(\mathrm{X})n_{\mathrm{H}}$. We take $x(\mathrm{He})={\delta }_{\mathrm{He}}$, with ${\delta }_{\mathrm{He}}=0.1$ and compute $\mathrm{H}$ abundance from the molecular fraction $f_m$ $n\left(\mathrm{H}\right)=(1-f_m-x\left({\mathrm{H}}^{+}\right))n_{\mathrm{H}}.$The resulting system of two equations can be written in a more compact form by using X for hydrogen and Y for helium $\begin{array}{rl}& {\alpha }_{rr}^X{n}_{\mathrm{H}}{X}^2+{\alpha }_{rr}^X{n}_{\mathrm{H}}XY+({\zeta }_{\mathrm{H}}+{\alpha }_{rr}^X{\delta }_{\mathrm{M}}{n}_{\mathrm{H}}+{\alpha }_{gr}^X{n}_{\mathrm{H}})X\\ & ={\zeta }_{\mathrm{H}}(1-f_m),\\ & {\alpha }_{rr}^Y{n}_{\mathrm{H}}XY+{\alpha }_{rr}^Y{n}_{\mathrm{H}}{Y}^2+({\alpha }_{rr}^Y{\delta }_{\mathrm{M}}+{\alpha }_{gr}^Y)n_{\mathrm{H}}Y={\zeta }_{\mathrm{He}}{\delta }_{\mathrm{He}}.\end{array}$This system is easily solved using a Newton-Raphson scheme once the rates are known.

Radiative electronic recombination rates are taken from [61], the relevant coefficients are given in Table A1 $\begin{array}{r}{\alpha }_{rr}=A\times [\sqrt{\frac{T}{T_0}}{(1+\sqrt{\frac{T}{T_0}})}^{1-BB}\times {(1+\sqrt{\frac{T}{T_1}})}^{1+BB}]\end{array}$with $BB=B+C\exp (-\frac{T_2}{T}).$

Table A1. Electronic recombination rates.

Coefficient	A	B	$T_0$	$T_1$	C	$T_2$
${\mathrm{H}}^{+}$	$8.32{10}^{-11}$	0.7472	2.96	$7.0{10}^5$	0	0
${\mathrm{He}}^{+}$	$5.23{10}^{-11}$	0.6988	7.3	$4.48{10}^6$	0.0829	$1.68{10}^5$

Electronic recombination on grains, comes from [62] ${\alpha }_{gr}(X^{+},\frac{G_0}{n\left(e^{-}\right)},T)=\frac{{10}^{-14}{C}_0}{\begin{array}{c}1+C_1{\psi }^{C_2}\\ (1+C_3{T}^{C_4}{\psi }^{-C_5-C_6\ln T})\end{array}}$with $\psi =\frac{G_0\sqrt{T}}{n\left(e^{-}\right)}.$Here, $G_0$ is the Inter Stellar Radiation Field (ISRF) intensity in units of Draine's ISRF. Coefficients $C_0$ to $C_6$ are given in Table (A2).

Table A2. Grain recombination rates.

Coefficient	$C_0$	$C_1$	$C_2$	$C_3$	$C_4$	$C_5$	$C_6$
${\mathrm{H}}^{+}$	12.25	$8.074{10}^{-6}$	1.378	$5.087{10}^2$	$1.586{10}^{-2}$	0.4723	$1.102{10}^{-5}$
${\mathrm{He}}^{+}$	5.572	$3.185{10}^{-7}$	1.512	$5.115{10}^3$	$3.902{10}^{-7}$	0.4956	$5.494{10}^{-7}$

Appendix 3

Fortran code

The code solving for ${\mathrm{H}}_3^{+}$ populations from Equation  (14(14) $\begin{array}{rl}& (\sum _{i>j}{A}_{ij}+n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ij}^o(1-p_2)+k_{ij}^p{p}_2)\\ & +\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X+{\alpha }_{\mathrm{DR},i}n\left(e^{-}\right))x_i\\ & -\sum _{i<j}{A}_{ji}{x}_j-n\left({\mathrm{H}}_2\right)\sum _{j\ne i}(k_{ji}^o(1-p_2)+k_{ji}^p{p}_2)x_j\\ & -\sum _{j\ne i,X}n\left(\mathrm{X}\right)k_{ij}^X{x}_j-k_{\mathrm{form},i}R=0,\end{array}$(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: $(n_{\mathrm{H}},T,f_m,I_m,\zeta )$ 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.