Rovibrational relaxation model for H + 3

. A rovibrational relaxation model based on existing levels and transition probabilities for H +3 is presented. In this model all H +3 levels below 12 000 cm − 1 and their respective transitions are included. The results are compared to Coulomb explosion measurements of the relaxation time of the ﬁrst vibrational breathing mode and found to yield excellent agreement. Furthermore, the existence of long-lived rotational states with energies up to 1 eV, as detected in DR imaging experiments at the TSR storage ring, is backed by the model. A reduced set of differential equations—including only states below 8000 cm − 1 — was set up, which incorporates also induced transitions. This model was used to investigate the heating of subthermal ensembles of H +3 ions exposed to 300 K blackbody radiation.

is the main reaction agent responsible for the formation of complex molecules in the reaction network of the interstellar medium.
The spectroscopy of H + 3 is hindered by a vanishing electric dipole moment and the lack of stable electronically excited states, leaving only vibrational transitions in the infrared region for examination. Due to rather weak transition strengths and the experimental difficulties that were encountered while working in the infrared region, it took several campaigns and many years of work until the first laboratory spectrum could be measured by Oka in 1980 [3]. Nowadays more than 800 vibrational transitions are documented (a comprehensive listing can be found in [4]) and the deuterated isotopomers D + 3 , D 2 H + and H 2 D + have also been subject to spectroscopic studies.
It took another 16 years until H + 3 could be detected in dense interstellar clouds [5] and it came as a surprise that only two years later a signal was seen also in a diffuse interstellar cloud revealing a similar H + 3 column density [6] as compared to those of dense clouds, in apparent contradiction to contemporary models of cloud chemistry. This discrepancy is still unresolved; however, if one considers the formation and destruction processes of H + 3 in dense and diffuse clouds, in both cases the crucial parameter controlling the formation is the cosmic ionization rate ζ, while the destruction proceeds along two different pathways. In dense clouds the metallicity is high and the fractional electron density is rather low, thus chemical reactions (e.g. with CO) provide the main destruction mechanism. Conversely, in diffuse clouds, the fractional electron density is much higher rendering the electron recombination of H + 3 , followed by the dissociation into neutral particles (dissociative recombination, DR), the dominant destruction path. This interrelation implies that the H + 3 dissociative recombination rate is one of the key parameters for cloud models, in fact, it was shown that relatively small changes in the DR rate coefficient α have the potential of severely altering the whole reaction scheme [7]. During the last few decades, several measurements of the DR rate coefficient have been performed (for reviews see e.g. [8,9]) utilizing various techniques; unfortunately, the results scatter by several orders of magnitude, despite small individual error bars. As one of the reasons for these conflicting results it has been suggested that different populations of internal molecular excitations might be present in these measurements; while the situation for rotational excitations is less clear, it is known that vibrational excitation can indeed drastically change the DR rate coefficient.
One of the more recent experimental approaches to determine DR rate coefficients is the storage ring technique, where an H + 3 beam is recirculated in an ultra high vacuum ring structure and a cold electron beam of variable energy can be superimposed to initiate the recombination process. The neutral fragments can be easily detected by suitable detector systems and permit to determine the rate coefficient and to examine the dissociation dynamics as well as the energy sharing between the fragments. Moreover, the long storage times of seconds up to minutes allow the cooling of internal excitations due to spontaneous decay. In the case of H + 3 not all vibrational excitations are expected to cool on the same time scale: H + 3 , which is an equilateral triangle in its ground state, hosts three normal-mode oscillations, one symmetric stretch motion often referred to as the breathing mode, and two degenerate asymmetric oscillations called the bending modes. While the latter create a small temporal dipole moment and thus decay spontaneously with lifetimes of the order of milliseconds, the breathing mode preserves the equilateral symmetry and consequently is metastable in character. Since from three ion-storage-ring experiments [10]- [12] the H + 3 DR rate coefficient was reported to be high as compared to some previous experimental and theoretical findings, it was suggested that long-lived vibrational states-which could survive in the essentially collision-free storage environment-might be responsible for the enhanced recombination, even though calculations [13] predicted lifetimes of distinctly less than one second for all vibrationally excited states except for the first breathing mode, for which a lifetime of τ = 1.18 s was calculated in the absence of any rotational excitation. According to these calculations, all H + 3 should have decayed to the vibrational ground state after a few seconds of storage, that is at the time the DR rate measurements took place.
To verify the vibrational cooling of H + 3 , a Coulomb Explosion Imaging (CEI) [14] measurement was carried out at the TSR storage ring of the Max-Planck-Institut für Kernphysik in Heidelberg, that demonstrated [15] the presence of vibrationally excited H + 3 states during the beginning of storage which, however, were shown to decay to the vibrational ground state on time-scales even shorter than predicted. On the other hand, complementary measurements at the same facility with a DR fragment imaging technique revealed that rotationally highly excited states are populated at all storage times of up to 60 s [16,17] with a distribution that could be modelled by a Boltzmann distribution with an average energy of E ≈ 0.3 eV. To sum up the experimental results, one was confronted with vibrational relaxations that were faster than the rotationless calculation predicted, while at the same time evidence was found for long-lived rotations with energies reaching up to 1 eV.
That trapping of energy in rotationally excited states of H + 3 can occur, which could even lead to maser action in astronomical environments, has been pointed out already by Black [18], and Goto et al [19] have found first observational evidence for such a population trapping by detecting very non-thermal populations of H + 3 in its (J, K) = (3, 3) rotational state. These results have recently been modelled by Oka and Epp [20].  [21] (only states with E f 5000 cm −1 are shown). In this work we present a comprehensive description of a relaxation model for rovibrationally excited H + 3 ions, which is based on the transition line list calculated by Neale et al [21] at University College London (UCL), and which has been employed already in [15] to explain the results of the CEI experiment. Emphasis will be laid on the discussion of the remaining rotational excitations after long cooling times of initially hot H + 3 ensembles, and on the extension of the model to investigate the influence of the blackbody radiation on cold ensembles.

3
The UCL line list (available via FTP [22]) comprises all the transitions between energy states in H + 3 below 15 000 cm −1 . The energy levels and wave functions were acquired variationally using the DVR3D program suite of Tennyson et al [23] and utilizing the discrete variable representation method to solve the nuclear motion problem. Special care was taken to adopt the basis set to cope with the triangular structure at low energies as well as possible linear arrangements above the so-called barrier to linearity at ∼10 000 cm −1 . The calculations employed the spectroscopically determined H + 3 surface of Dinelli et al [24] and the ab initio dipole surface of Röhse et al [25]. The convergence error for the energy levels was estimated to be ∼0.05 cm −1 in the low energy region <5000 cm −1 , ∼0.2 cm −1 at 10 000 cm −1 and ∼1 cm −1 at 15 000 cm −1 .
On the basis of the computed levels, the corresponding wave functions and all possible dipole transitions were determined taking into account only the strict selection rules for angular momentum J = 0, ±1 and parity change ± ↔ ∓. From the complete calculation of 3 187 125 transitions, those with an Einstein A-coefficient of less than 10 −7 s −1 were removed, reducing the set to 3 070 571 lines.
In the list (table 1 shows an excerpt) the following parameters are given to characterize each particular transition: the angular momentum quantum number of the initial (J i ) and final state (J f ), the respective energies (E i , E f ), the transition frequencyν if , the Einstein coefficient of spontaneous decay A if , and the nuclear spin degeneracy g i , if known; otherwise, the average value of 8 3 [26] was chosen. (The nuclear spin assignment was not done automatically by the program code, thus for most of the states below 5000 cm −1 the nuclear spin degeneracy g i was retrieved by an algorithm starting from a few hand-assigned states.) For more details on the line list see the original publication [21] and references therein. Note that the excitation energy scale is chosen such that the absolute ground state of H + 3 is at 64.126 cm −1 (with respect to the non-existent J = 0 state).
The vibrational quantum numbers (ν 1 , ν |l| 2 ) as well as the quantum number K, representing the projection of the angular momentum vector on to the molecular axis, are not given in the UCL line list. However, a second publication by Dinelli et al [27] where all the levels up to 9000 cm −1 and J 9 are tabulated together with all their quantum numbers, allowed a manual identification of the levels in the UCL list up to this energy.

Relaxation and initial populations
In order to trace the relaxation of internal excitations in H + 3 , a set of differential equations has to be generated that connects all the relevant rovibrationally excited levels according to their spontaneous decay rates. In a first step a PERL script was used to browse through the complete UCL line list in order to extract all the energy levels. The search was restricted to levels below 12 000 cm −1 (∼17 000 K), firstly to reduce the number of transitions to ∼245 000 (and the number of relevant states to ∼2500) and secondly because the levels above 9000 cm −1 could not be uniquely assigned (see section 2). In addition, the lifetime T i of each state i was determined by where the sum was taken over all states with E f < E i . For states where no decay path was found, the lifetime was artificially set to 1 × 10 7 s. Denoting the population of the state i at time t by P i (t), the differential equation for P i (t) is given by where the sum runs over all states with E k > E i . The coupled set of equations was solved using the DRKSTP routine [28] of the CERN program library, which employs the Runge-Kutta method for numerical integration of the first-order differential equations [29]. The initial population P i (0) of each of the ∼2500 states was assumed to be determined by a Boltzmann distribution of temperature kT , i.e.
where the normalization N 0 is chosen such that i P i (0) = 1. It should be noted that the nuclear spin degeneracy g i could not be retrieved for most of the high-J states (e.g. all states above 5000 cm −1 ); in these cases the average value of 8 3 was adopted. Although the transitions within the respective ortho-or para-manifold are not affected, this procedure might lead to a slight overpopulation of para-states.

Incorporating radiative heating
So far the 300 K radiation, present in ion storage rings operating at room temperature, is neglected in the relaxation model described above. This is justified as long as the population distribution is characterized by an effective temperature much larger than 300 K. However, in more recent storage ring experiments, where subthermal ion ensembles are being stored (created either by buffer gas cooling in the ion source itself [30] or by electron cooling in the storage ring [31]), it is also important to estimate typical heating times caused by the exposure to the 300 K radiation and to find out which states are involved in the heating process. A reduced model covering all states below 8000 cm −1 (11 500 K) was therefore set up, which includes the coupling to the ambient radiation field. This reduced model still contains 693 states and ∼15 600 transitions.
The Einstein coefficient A 21 , describing the spontaneous decay from level 2 to level 1, is connected to the Einstein coefficient of induced emission, B 21 , and induced absorption, B 12 , by where ν 12 is the transition frequency. Denoting the spectral energy density of the radiation field by ρ(ν), the differential equation for the time-dependent population of state i is now given by where j runs over all f and k states connected to state i. Again, the DRKSTP routine [28] of the CERN program library was used to solve the set of coupled equations. The spectral densitygiven by Planck's radiation law-was calculated by an external function for each transition frequency ν ij and transferred as a vector to the differential equation solver (for details see [32]).

Relaxation of hot H + 3 ensembles
The calculations were performed on a Compaq Alpha Es40 processor running at 500 MHz. Time steps of 1 ms were chosen for the Runge-Kutta routine after convergence checks had been done with both smaller and larger step sizes. A typical calculation, starting with a Boltzmann distribution for kT = 0.23 eV and following the relaxation for 60 s, took approximately 1 h. The longest computation that was done with the full set of equations-for a relaxation time of 1 hour-took about three days of CPU time.
In figure 1 a stick diagram of the initial population resulting from a Boltzmann distribution with kT = 0.23 eV is plotted. Initial temperatures of this size were inferred from DR fragment imaging measurements performed at the TSR storage ring [16] when producing the H + 3 ions in a 'hot filament' ion source. Tests with different starting conditions between kT = 0.1 and 0.3 eV showed, however, that the general behaviour of the rovibrational decay as presented in this subsection is rather insensitive to the precise value of the initial temperature, since the higher states decay rapidly resulting in similar distributions after a few ms. For the same reason the cut of the initial distribution at E i > 1.5 eV, which is caused by the limitation of the relaxation model to level energies of <12 000 cm −1 , is expected to be unimportant with respect to the later relaxation behaviour. Figure 2 shows the relaxation of the kT = 0.23 eV hot H + 3 ensemble in six time slices (an animation of the cooling within the first minute is available with the online version of this text). To give a more realistic picture of the energy stored in the H + 3 ions, the stick diagrams displaying P i (t) as a function of E i have been converted into histograms P(E, t) by summing up the population within 6.5 meV wide energy bins. The displayed time slices of 0.1, 0.3, 1, 10 and 60 s span the typical time range of a storage ring experiment and reveal that a considerable amount of energy is present in the H + 3 ensemble even after 1 min of storage: the average energy E carried by the H + 3 ions (which is different from the temperature because the density of states increases dramatically with energy) drops rapidly from the initial value of 0.74 to 0.35 eV in the first second, while it stays as high as 0.31 eV after 10 s and 0.29 eV after 1 min, in agreement with the average excess energy of ≈0.3 eV observed experimentally at storage times greater than a few seconds [16,17]. The explanation for these two different time constants can be readily traced to the fast decay of vibrationally excited states and the slow relaxation of a group of only rotationally excited levels.
As already discussed above, H + 3 can host two principally different modes of vibrational excitation, the doubly degenerate bending motion and the symmetric-stretch or breathing mode. The creation of a temporary dipole moment renders short lifetimes of the bending modes of a few ms, while the decay of the breathing mode, being infrared inactive, is considerably slower. The decay of the longest lived vibrational excitations, which are based on the first and second excited breathing level, could be measured at the TSR storage ring using the CEI technique and was discussed already in [15] within the framework of the present relaxation model. For this comparison all the states (ν 1 , ν l 2 ) carrying one or two breathing mode quanta ν 1 were manually identified using the list of Dinelli et al [27]. In this way all levels below 9000 cm −1 were assigned, leading to 248 states with ν 1 = 1 and 73 states with ν 1 = 2. The result of the relaxation model for the absolute size and the time dependence of the fractional population of the first breathing mode is compared in figure 3 to the results of the CEI experiment [15]. Moreover, a single-exponential-representing the rotationless calculation of Dinelli et al [13] for the decay constant of the (1, 0 0 ) level-is shown, which was adjusted to the experimental values at storage times greater than 500 ms. The comparison clearly shows that at time <500 ms the observed decay of the first breathing mode proceeds on a time scale which is considerably faster than predicted by the rotationless calculation, but which is reasonably well accounted for by the present relaxation model involving also rotational excitations. The remaining slight underestimation of the absolute size of the fractional population by our model might well be caused by our choice of a Boltzmann distribution for the initial population, which is motivated by simplicity rather than its ability to reflect exactly-not necessarily thermal-the initial excitation distribution produced by the ion source.
There are two mechanisms whereby the metastable ν 1 = 1 vibrational state can decay. The first is by the symmetry allowed but weak vibrational difference band to the ν 2 = 1 state [33]; this is a purely vibrational effect. As to the second mechanism, introducing rotational motion allows coupling between different vibrational states and hence intensity stealing. For ν 1 = 1 this effect is particularly strong for J values of about 7. Experimentally, the effect has led to the observation of infrared hyforbidden transitions in the laboratory [34]. We consider this rotational coupling to be the main driving force for the observed speeding up of the ν 1 = 1 decay. The most striking feature of the distributions shown in figure 2 are the seemingly equidistant structures appearing at longer relaxation times. Since even the most long-lived breathing mode excitation decays within the first few seconds, it is clear that the patterns are caused exclusively by rotational excitations. To focus on the long-lived states, even longer relaxation times than accessible in the storage ring experiments were computed. Figure 4(a) shows a stick diagram of the states remaining after a relaxation time of 1800 s, again starting from a Boltzmann distribution with kT = 0.23 eV to fix the initial populations P i (0). Qualitatively, the picture is similar to that after 60 s, but at the longer storage time the doublet structures at higher energies are even more clearly born out. Although the published assignments end at 9000 cm −1 , from the level systematics it is simple to assign the K and J values to all the levels in this picture. In table 2, all states with a population 0.01% after 1800 s are tabulated, accounting for 99.98% of the total population in the model. In the table, the energy E, the population P after 1800 s, the lifetime T , the angular momentum quantum number J together with its projection K on to the molecular axis as well as the nuclear spin degeneracy g are given.
The origin of the doublet structures at higher energies manifests itself in the fact that for high J values always two states survive, namely those with the two highest possible K-values, K = J and K = J − 1. Figure 4(b) shows a zoom at the high-energy region of the distribution together with the level assignment. It is interesting to note that for states with J > 7 and K = J − 2 the decay almost exclusively proceeds via (J, J − 2) → (J + 1, J + 1) transitions, which are 'allowed' by the selection rules for rotationally forbidden transitions studied here [35]. Thus only a doublet of metastable states persists for all angular momentum states J > 7.
Furthermore, table 2 reveals the general trend that the levels with K approaching J are longlived, while those with K J decay fast. This again can be explained by simple considerations of the molecular symmetry: Since H + 3 is an equilateral triangle in its ground state, rotations around the symmetry axis, i.e. states with J = K, may experience a stretch but with the symmetry undisturbed, whereas with decreasing K the triangular shape can be increasingly distorted, resulting in the formation of a small dipole moment.  + 3 , where E denotes the energy, P the population after 1800 s, T the lifetime, J the angular momentum quantum number, K the projection of J on to the molecular symmetry axis and g the nuclear spin degeneracy. All states with population 0.01% after a relaxation time of 1800 s (starting from a Boltzmann distribution for kT = 0.23 eV) are listed. Note that all the states are purely rotational, e.g. ν 1 = ν 2 = 0 and especially K = G. The states listed account for 99.98% of the total population. For the states with J > 12, the K quantum number and the nuclear spin degeneracy g were hand-assigned according to the level systematics. In light of the relaxation model, the evidence for excess energy of ≈0.3 eV-compatible with a Boltzmann distribution of ≈0.23 eV temperature-that was found in the DR imaging experiments [16,17] is exothermic by 1.7 eV, hence enough energy is available to populate even high-lying rotational states. A simple estimate shows that if a collision such as formulated in equation (6) takes place in a hot ion source, easily enough angular momentum can be picked up to populate all the high-J states that are listed in table 2. Test beamtimes at the TSR, employing three different types of ion Stick diagram displaying the population of H + 3 levels assuming an initial thermal equilibrium population for a temperature of 100 K, and after exposure to blackbody radiation of 300 K for 10, 100 and 3600 s, respectively. sources under various pressure and temperature conditions, revealed very similar H + 3 excitation energies, leading us to the conjecture that the energy release in the formation reaction dominates the initial excitation pattern rather than the thermal conditions in the ion source.

Radiative heating of cold H + 3 ensembles
To analyse the radiative heating of a cold H + 3 ensemble, an initial population of states in H + 3 corresponding to a subthermal Boltzmann distribution of kT = 100 K temperature was processed with the reduced model described in section 3.2. To check for consistency, the result of an exposure to a 300 K blackbody radiation for 10 h was compared to a 300 K Boltzmann distribution [32]. Good agreement was obtained, although small deviations in the relative populations revealed that even longer times are required for complete thermalization.
To illustrate the reheating time scales, figure 5 shows the development of the initial 100 K distribution under the influence of the room temperature radiation field in three time slices of 10, 100 s and 1 h. While the first three distributions are practically identical at first glance, the last plot shows that higher states slowly start to be populated. To quantify the changes, one may consider the five lowest-lying states which account for 99.0% of the total population in the initial Table 3. List of the 20 strongest transition rates for the five lowest-lying states of H + 3 in a 300 K radiation field. E i , J i , G i , U i and E f , J f , G f , U f denote the energies and rotational quantum numbers of the initial and final state, respectively. The vibrational quantum numbers ν 1 ν |l| 2 of the final state are also given, while all initial states are in the vibrational ground state (ν 1 = ν 2 = |l| = 0 and K = G). Although these results indicate that there is no noteworthy radiative heating to be expected in the present storage-ring experiments, it is interesting to examine which transitions drive the slow reheating process visible in figure 5. In table 3 the 20 strongest heating transitions in a 300 K radiation field, heating up the five lowest excited H + 3 states, are listed together with the quantum number of the lower and upper states and the excitation rate in s −1 . Remarkably, only one of these transitions is purely rotational, whereas all others are populating the first vibrational bending mode (0, 1 1 ) which opens at 2521 cm −1 (equivalent to 3600 K). Obviously, the transition strengths to states involving the first bending mode are so much stronger than rotational transitions that they even compensate the exponentially decreasing spectral density of the 300 K blackbody radiation at the corresponding transition energies. This behaviour reflects the strong suppression of pure rotational transitions which are forbidden in the picture of separable wavefunctions and gain signal strength only by intensity stealing through level mixing with vibrationally excited levels [36].
This finding provoked considerations whether stray light of higher temperatures might be much more effective in heating up internally cold H + 3 in a storage ring than the room temperature radiation field. As possible light sources, the hot filament ion gauges used, e.g., in the TSR storage ring were spotted. Similar to light bulbs, these ion gauges illuminate their surroundings with filaments powered by 8-12 W. To scrutinize their effect, a worst case estimate of the spectral density created by these gauges was added to the 300 K field and the calculations were repeated. As a result, the heating effects were considerably faster [32] (the population of the five lowest states dropped to 92.1% after 100 s) but still too slow to play an important role in a typical storage ring experiment. On the other hand, the calculation demonstrated that high temperature stray light can easily be the dominant source of radiative heating for a system like H + 3 .

Conclusions and outlook
Several conclusions can be drawn from the relaxation model as described above. Firstly, it resolves the discrepancy between the experimentally measured relaxation times of the vibrational breathing modes and the rotationless calculation, which is clearly not adequate if the ions are produced in a 'hot' ion source. Explicitly, it is shown that the presence of rotationally excited states speeds up the vibrational decay substantially. Secondly, the existence of highly excited rotational states in storage-ring experiments even after several ten seconds of storage is reproduced by the model and the origin and properties of the long-lived states are readily understood. Considering the lifetimes of these states, is is apparent that the preparation of cold H + 3 ensembles cannot be accomplished by longer storage times, since many of the J = K and J = K + 1 levels have no decay route and might survive for many months in a collision-free environment. This finding supports the idea that these states may mase in appropriate astronomical environments [18].
Furthermore, the exothermic formation reaction is likely to be responsible for most of the internal energy and angular momentum of the molecular ions, so that an active mechanism is necessary for successful cooling. The most promising approach is to use collisions with a cold inert buffer gas to thermalize the ions. At the Cryring facility in Stockholm a supersonic expansion source has recently been used to produce a cold H + 3 beam for injection into the storage ring [30]. The measurement revealed a thermal DR rate coefficient that is about a factor of 2 lower than observed previously at the same setup. A different concept is under development at the TSR in Heidelberg; here a cryogenic 22-pole trap was constructed [37] in collaboration with the Technische Universität Chemnitz in order to set up a temperature-variable ion source for the temperature range found in cold interstellar clouds (10-100 K). Another option for the production of subthermal ion beams in storage rings is to use the electron cooler to cool also the internal degrees of freedom of the molecular ion. In experiments with D 2 H + -the doubly deuterated sibling of H + 3 -the preparation of subthermal temperatures by electron bombardment could be demonstrated at the TSR [31].
While the heating of cold ensembles of infrared-active ions like H + 3 in a thermal 300 K radiation field is rather small even at very long storage times, a 'heating up' of the ions due to residual gas collisions may occur. If the time-scale for spontaneous relaxation is slower than the production of excited states by collisions with the diffuse ambient gas, the effective temperature of the ions circulating in the storage ring may rise with increasing storage times. Evidence for a slow heating process due to residual gas scattering was recently found at the TSR in experiments with H + 3 [38] as well as 4 He 4 He + ions [39]. Coming back to the dissociative recombination rate of H + 3 , it follows from the above discussion that the population of individual states in H + 3 at the time the recombination takes place depends strongly on the history of the H + 3 ion and on its environment. An additional issue that has not yet been thoroughly addressed experimentally is the influence of nuclear spin modifications on the DR process. While radiative transitions between ortho-and para-H + 3 are forbidden, collisions on the other hand may lead to a change of the nuclear spin. Cordonnier and coworkers [40,41] could show that especially the chemical reaction of H + 3 with the neutral hydrogen molecule is an effective pathway for ortho-para conversions through the temporary formation of an H + 5 collision complex that may lead to a scrambling of all five protons. At low temperatures (<20 K) the (J = 1, K = 1) para-state being about 30 K lower in energy than the next lowest (ortho) state (J = 1, K = 0) is supposed to be dominant, and the most recent theoretical calculations predict a difference of a factor of ∼2 in the DR rate coefficients of the two states at 10 K [42]. Taken together, one is tempted to conclude that the electron recombination of a polyatomic molecule like H + 3 is simply too complex to be characterized by a thermal rate coefficient alone. In fact, it is possible that a profound knowledge of all the states and rate coefficients involved is necessary to eventually understand the various laboratory experiments and to unravel the role of H + 3 DR in interstellar clouds.