Understanding the temperatures of H3+ and H2 in diffuse interstellar sightlines

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


Introduction
The triatomic hydrogen ion 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 Owing to the comparatively low proton affinity of 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 included in the PDR code to account for the ortho-para character of H + 3 .The paper concludes with a brief discussion in Section 6.

Nuclear spin and rotational states of H + 3 and H 2
Both H 2 and H + 3 exist in two different nuclear spin configurations.For the lowest rotational state of H 2 , with J = 0, the proton spins are anti-parallel and add up to I = 0.This configuration is denoted as para-H 2 (or p-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 H 2 states with odd rotational quantum numbers have parallel nuclear spin configurations with I = 1, which is denoted as ortho-H 2 (or o-H 2 ).Likewise, all even rotational quantum numbers in H 2 can be attributed to p-H 2 with I = 0. Figure 1 shows the three lowest levels of H 2 (with J ≤ 2) and their respective energies expressed as k b T (in units of kelvin).
The excitation temperature of the two lowest states of H 2 is defined as where ∆E 01 /k b = 170.476K 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 ∆E 02 /k b = 509.864K 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 K [15,[19][20][21].
The nuclear spin configurations of H + 3 are also denoted by para and ortho for I = 1/2 and I = 3/2, respectively.As with 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 H + 3 , while all other levels (effectively fulfilling G = 3n ± 1) are of the para configuration, with I = 1/2.The right-hand side of Figure 1 shows all H + 3 rotational levels with J ≤ 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 H + 3 , we will from now on refer to all 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 center [12,23,24], only the lowest (1, 1) para-state of H + 3 and the lowest ortho-state (1, 0) have ever been detected in the interstellar medium.Consequently, the H + 3 excitation temperature, T 12 (H + 3 ), derived from observations is usually calculated using the column density ratio of these two states ) (1,0) Level scheme of all rotational levels with J ≤ 2 for H 2 and with J ≤ 3 for H + 3 .The level energy (in cm −1 ) is given on the y-axis at the left-hand-side, while the corresponding values in 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.
where ∆E/k b = 32.86K, 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 (H + 3 ) between 20 and 40 K [16,17], systematically lower than the H 2 excitation temperature T 01 .As para-and ortho-states of H 2 and H + 3 can not be inter-converted by radiative transitions, we now describe the various collisional and reactive processes allowing ortho/para exchange.

Processes controlling the para fraction of H
3 formation via the H + 2 + H 2 reaction (Eq. 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 − 60 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 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 × 10 9 cm 3 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 H + 3 ions created by reaction (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 H + 2 and 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-H 2 ) and n(o-H 2 ) of the para-and ortho-species of H 2 , and n(p-H + 3 ) and n(o-H + 3 ) of H + 3 , respectively.The para-fraction for both species is then denoted by and Following the argumentation by Crabtree et al. [16], we assume that the cosmic ray ionization of H 2 (constituting the main source of ionization that initiates the H + 3 formation) does not affect the nuclear spin of the molecule, and therefore the para-fraction of H + 2 has the same value p 2 .With these assumptions the p-H + 3 fraction at formation, p f 3 , 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 The values of p 2 and p f 3 are displayed as a function of temperature in Figure 2 for H 2 in thermal equilibrium.
We further assume that the rotational levels of p-H + 3 and o-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.

Thermalizing collisions of H + 3
Collisions with electrons, He, H and H 2 contribute to the energy exchange between the rovibrational levels of H + 3 , but only collisions with H 2 and H may change the nuclear spin of the H + 3 ions.To the best of our knowledge, of the above collision partners, only for collisions with H 2 and electrons detailed information is available in the literature.

Collisions between H +
3 and H 2 Before any detailed study of H + 3 collision rates with H 2 was available, Oka and Epp [31] suggested to use the Langevin expressions to describe the excitation of H + 3 detected towards the Galactic center.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 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 H + 3 and 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 H + 3 -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 H + 5 potential energy surface [37].We employ a set of rate coefficients that covers collisions of the lowest 24 rotational states of H + 3 with ortho-H 2 and para-H 2 (in their lowest rotational states) and temperatures up to 500 K, which were kindly provided by O. Roncero.Those rate coefficients are also currently used in our PDR model calculations [18].

Collisions between H +
3 , He and H He and H are additional collision partners that should be taken into account when describing the excitation equilibrium of H + 3 , but, to the best of our knowledge, no detailed studies are presently available on these systems.Collisions of H + 3 with He can not modify the nuclear spin configuration of H + 3 , and we approximate the corresponding collision rates by taking the rates for p-H 2 and scaling them with the reduced mass factor (while vetoing channels that might change the nuclear spin of H + 3 ), as described previously [38].
Collisions of H with H + 3 may lead to proton exchange or inelastic collisions, similar to the channels discussed for the H + 3 + H 2 reaction.However, since we are not aware of more detailed information on this process, we will use the collisional rates with p-H 2 as a proxy for collisions with H.We checked that this choice has no major influence on our model calculations for the conditions studied here.

Inelastic collisions between H +
3 and electrons Electronic collisions with H + 3 preserve the ortho-para character of 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 CH + allowed for a comparison between experiment and theory for the lowest rotational states [40], which revealed very good agreement.

Dissociative recombination of H + 3
In diffuse gas, the main destruction reaction of H + 3 is the dissociative recombination (DR) with electrons.
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 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][42][43]).In summary, the present consensus is that the absolute rate of the 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][45][46].
Both astrochemistry databases, UMIST2 and KIDA3 , report rate coefficients with branching ratios of 2/3 for reaction (10a) and 1/3 for reaction (10b), based on storage ring experiments, and the total DR reaction rate coefficient is given as However, there are recent suggestions concerning a possible difference between the DR rate coefficients of the two nuclear spin modifications of 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-H + 3 recombining much faster than o-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 Table 1 summarizes the values of the principal reaction rate coefficients introduced in the present study, where we have assumed the same branching ratios for reactions (10a) and (10b) as described above.

Master equation for H + 3 state populations
In this section, we describe the master equation for the level populations of H + 3 , and we show that 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 (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] 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-H + 3 by p-H 2 followed by incomplete thermalization is responsible for the deviation from thermal equilibrium.

Differential equations
The variation over time of the density of a level i (i ∈ [1, N ]) of H + 3 , n i , can be described by (formation and destruction) (13) (collisions with other species X) with k form,i and α DR,i denoting the state-dependent chemical formation and destruction rates, k p,o ij the collisional excitation/de-excitation rates with p-H 2 and o-H 2 , and k X ij denote collision rates where X stands for H, 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 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 H + 3 , so that n i = x i n H + 3 , the total destruction rate is α DR = i x i α DR,i .The total formation rate is k form = i k 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 α DR,i differ for individual states, and thus n H + 3 cannot be computed independently.At steady state, dni dt = 0 and the system of equations to solve is . The value of R is initially unknown, as the total destruction rate of H + 3 requires the knowledge of the relative level populations of the molecular ion.So, we explicitly add the conservation equation We get a system of N + 1 equations with N + 1 unknowns, which is easily solved if the densities of H 2 , 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.
Derivation of H + 2 density and electronic fraction It is possible to reduce the set of equations further by estimating n H + 2 and n (e − ).The main formation and destruction reactions of H + 2 are: where ) 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 C + .However, for high cosmic ray ionization 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 ionization fraction.
Here we compute the electronic density considering both H + and He + as described in App.B, which involves the relative abundances of all atoms with ionization potential below the Lyman cutoff (assumed to be fixed) and the UV radiation field 4 .With these derivations of n H + 2 and n (e − ), the system formed by Eq. ( 14) and Eq. ( 15) depends on five astrophysical parameters: n H5 , T , f m , G 0 , and ζ, where we introduced f m the molecular fraction, f m = 2 n(H 2 )/n H .

Numerical solution of the coupled equations
We developed ExcitH3p, a FORTRAN code that solves the H + 3 coupled system of equations (Eq.14) for any number of levels, and then computes the two-level excitation temperature (Eq.4) 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 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 K above the ground state.We find that the solution to the set of equations ( 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 ionization rate ζ and gas kinetic temperature T .The total gas density is set to n H = 10 2 cm −2 and the molecular fraction to f m = 0.8.We find that in the entire parameter space, T 12 (H + 3 ), the excitation temperature given by the two lowest 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 H + 3 , which are detailed above.Figure 3 also shows that a maximum of about 44 K is reached for T 12 (H + 3 ) at kinetic temperatures around 100 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-H + 3 are the same: T 12 (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 H + 3 levels included in the model.Figure 4 shows the influence of density, showing that T 12 H + 3 is not sensitive to this parameter up to densities of more than 10 3 cm −3 if the gas temperature is below ∼ 80K.

Influence of the number of 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 H + 3 levels included in the model.Except for the Galactic center, 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 (H + 3 ) value as a function of the kinetic temperature T , for different values of N .We perform this comparison for the following parameters: ζ = 5 10 −16 s −1 , n H = 100 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 (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)

population.
The decrease of T 12 (H + 3 ) with increasing gas kinetic temperature above 100 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 H 2 shows that its branching ratios towards ortho and para are approximately equal in the 30 − 300 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).

Two-level approximation
In a restricted range of kinetic temperatures around T 25 − 50 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 o 12 (1 − p 2 ) + k p 12 p 2 , the total collisional excitation rate due to H 2 from level 1 to level 2 (and the equivalent for k 21 ).Then, the coupled equations system (Eq.14) reduces to: Introducing k form = k form,1 + k form,2 , we can compute the x 2 /x 1 ratio: The factor R cancels out and the ratio x 2 /x 1 does not depend on n H + 3 .We introduce the electronic fraction x e = n (e − ) /n H and the configuration-specific formation rates of H + 3 with k form,1 = p f 3 k form and k form,2 = (1 − p f 3 ) k form . Furthermore, we apply detailed balance to the H + 3 -H 2 collisional rates, yielding k 12 = g2 g1 exp − E21 T k 21 .Finally, we get the following expression 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 H 2 , and the electronic fraction x e .It is independent of the cosmic ray ionization 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 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+2 p 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 x2 x1 ratio compared to the thermal value (in agreement with astronomical observations).The extent of the deviation from thermal equilibrium depends critically on the ratios αDR,1 k21 and αDR,2 k21 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 thermalization -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 α DR,1 and α DR,2 .
The picture that emerges is that the excitation temperature of H + 3 appears lower than the nominal gas temperature because of incomplete thermalization.The formation process strongly favors the formation of p-H + 3 at low temperatures, as the para-fraction p 2 of H 2 is large.The electron recombination process then removes H + 3 before it can reach thermal equilibrium in collisions with 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 H + 3 −H 2 collisional rate coefficients (or sufficiently reduced electron recombination rates) the calculations reach complete thermal equilibrium between the excitation temperature T 12 H + 3 and the kinetic temperature T .In essence, the nuclear spin restrictions of the H + 3 formation reaction produce an overproportional amount of H + 3 in the para configuration, and thermalization 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 H 2 , where the destruction process is much slower compared to thermalizing collisions (see Appendix A for details).

Master equation approach
To validate our approach, we compare the results of the master equation approach (Eq.14) with observations.Table 2 presents the H + 3 excitation temperature reported in the literature for a few local diffuse clouds as well as the corresponding H 2 excitation temperatures.Data for 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 H 2 data.Determination of diffuse cloud density with observations of H and 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 H 2 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 Eq. ( 14), we compute T 12 H + 3 for all the lines of sight.We assume a cosmic ray ionization rate of ζ = 10 −16 s −1 and a molecular fraction of f m = 0.8.For the gas density, the values in Tab. 2 are used, and for the gas temperature we use the values derived from H 2 observations T obs 01 (H 2 ).There is a single noticeable exception: X Per.However, T 12 H + 3 of this line of sight suffers from a particularly large error bar as listed in Table 2 .A high excitation temperature of H + 3 can only be reached for a kinetic temperatures close to 100 K with our models, as can be seen in Figure 3.

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 optimize the free parameters to reproduce other observations than H 2 and 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 H + 3 via H 2 + H + 2 collisions and the nuclear spin dependence of the dissociative recombination rate coefficient, in addition to the other excitation mechanisms of H + 3 .The collisional excitation of H 2 by H + , that has been revisited by [54] for highly rovibrationally excited levels, has also been updated.Observations suggest a total visual extinction A V = 0.9 mag (R V = 2.8, E B−V = 0.32, N H /E B−V = 5.2 10 21 cm −2 ).In our models the cloud is illuminated from both sides by the standard ISRF (G 0 = 1).
We consider two different scenarios Model A Constant density and constant temperature, Model B Constant density, with computation of the thermal balance.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 H + 3 and H 2 in the diffuse ISM can, in fact, be predicted considering a a strongly reduced set of reactions.

Discussion and conclusion
In this paper, we review all physical processes relevant to formation, excitation and destruction of 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 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 modelers 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 H + 3 .Specifically, the formation of p-H + 3 by p-H 2 at temperatures below 70 K is efficient.Considering the individual levels of H + 3 , we find that reactive collisions with 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 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 + 3 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 2 -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 + 3 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 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 % of 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 un-Additional electrons come from H + and He + resulting from the balance between ionization via cosmic rays and recombination with electrons and grains.We follow here for the most part the presentation of [60] For all abundances, we write x (X) = n (X) n H .We take x (He) = δ He , with δ He = 0.1 compute H abundance from the molecular fraction f m n (H) = 1 − f m − x H + n H .
The resulting system of two equations can be written in a more compact form by using X for hydrogen and Y for helium α X rr n H X 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 B1 Electronic recombination on grains, comes from [62] α gr X + , G 0 n (e − ) , T = 10 −14 C 0 1 + C 1 ψ C2 (1 + C 3 T C4 ψ −C5−C6 ln T )

Figure 2 .
Figure 2. Fraction p 2 of p-H 2 and p f 3 of nascent p-H + 3 (Eq.8).The typical range of diffuse cloud kinetic temperature is outlined in yellow.

Figure 3 .
Figure 3.The excitation temperature T 12 (H +3 ) as calculated using ExcitH3p, as a function of ζ and T for a density n H = 10 2 cm −3 and a molecular fraction fm = 0.8.Note that the lowest 24 levels of H + 3 are included in the model, while only the lowest two states are used to determine T 12 (H + 3 ).

Figure 4 .Figure 5 .
Figure 4.The excitation temperature T 12 (H + 3 ) as a function of the gas kinetic temperature T and density n H for ζ = 10 −16 s −1 , a molecular fraction fm = 0.8 and considering the first 24 levels of H + 3 .
, Section 13.6, extended to include He.The balance equations are ζ H n (H) = α rr H + n H + n e − + α gr H + n H n H + , ζ He n (He) = α rr He + n He + n e − + α gr He + n H n He + .Where α rr is the radiative recombination rate, and α gr the rate of recombination on grains.ζ H and ζ He are the cosmic ray ionization rates of H and He, respectively.With respect to H 2 , we use ζ H = 0.77 ζ, including secondary ionization, and ζ He = 0.5 ζ.The total abundance of electrons is n e − = δ M n H + n H + + n He + .

Table 1 .
Principal chemical reactions for H + 3 formation and destruction.
CRP represents cosmic ray particles and ζ the corresponding ionization rate of H 2 .At steady state, this leads to (with α H + 2 the dissociative recombination rate of H + 2

Table 2 .
Selection of sightlines with both H 2 and H + 3 observations.valueunavailable.Educated guess only.Model results and sensitivity to the DR rate coefficientsFigure 6 shows the results for three different hypothesis for α DR (o-H + 3 ), and the value provided by equation (12a) for α DR (p-H + 3 ).We see a quasi-linear variation of T 12 H + 3 with T 01 , which is well reproduced by the models.Using α DR from equation.(12b)(red circles) leads to temperatures which are too high, while using the same rate for o-H + 3 and p-H + 3 (pink circles) leads to temperatures which are too low.An empirical adjustment using α DR (p-H + 3 ) from Eqs. (12a) and α DR (o-H + 3 ) = α DR (p-H +3 )/1.5 (green circles) gives a very satisfying result, given that no attempt has been made to optimize other parameters. *

Table 3 .
PDR model results.Column densities are in cm −2 and excitation temperatures in K. Numbers in parenthesis are power of 10.
2 + α X rr n H X Y + ζ H + α X rr δ M n H + α X gr n H X = ζ H (1 − f m ) , α Y rr n H X Y + α Y rr n H Y 2 + α Y rr δ M + α Y gr n H Y = ζ He δ He .