Direct dissociative recombination of HCO⁺via the core-excited doublet and the lowest quartet states

Abstract

The dissociative recombination of HCO⁺ up to 1 eV collision energy is studied. New calculations for several core-excited HCO states provide improved potential energy surfaces crossing the HCO⁺ ground state surface in the vicinity of its equilibrium geometry. Wave packet analysis leads to significantly higher contributions of the direct mechanism to the cross section for electron energy ε < 0.7 eV than according to earlier studies [Larson et al., Phys. Rev. A, 2012, 85, 042702]. The limit H + CO(a³Π) is found to be the most probable exit channel. We discuss the improved agreement of theory with the latest experiments [Hamberg et al., J. Phys. Chem., 2014, 118, 6034] resulting from combination with the most recent calculations of the indirect process [Fonseca dos Santos et al., J. Chem. Phys., 2014, 140, 164308]. For the lowest quartet surfaces some vibrational states and their population and depopulation (mediated by spin–orbit coupling) are examined.

---

1. Introduction

The formyl cation HCO⁺ is one of the most abundant polyatomic ions in interstellar clouds and planetary atmospheres. The observation of its rotational spectrum may serve as a sensitive indicator for density and hence of the star formation rate in molecular clouds. In cold clouds, where the H₂ concentration is not directly measurable, the most important production and destruction reactions of HCO⁺ are H₃⁺ + CO → HCO⁺ + H₂ and HCO⁺ + e⁻ → H + CO respectively. Both the H₃⁺ and the electron concentration depend on cosmic ray ionization of H₂, so the (linear) HCO⁺ ion tends to be a rather reliable hydrogen density marker.¹

Dissociative recombination (DR) mechanisms of molecular ions with free electrons are generally classified as direct² or indirect.³ The direct mechanism is usually dominant, if the potential energy surface (PES) of a (doubly) excited dissociative state crosses the ground state surface of the ion not too far from its equilibrium geometry. The captured electron may then either be re-emitted (autoionization) or the neutral molecule dissociates along the PES. The characteristic for the indirect process is the temporary population of rovibrationally excited Rydberg states followed by predissociation or re-ionization. Both mechanisms may proceed under the same conditions; the total DR cross section can then be evaluated as the incoherent sum of the two processes.

Most theoretical studies agree that in the low electron energy limit for HCO⁺ (ε ≲ 0.1 eV) the indirect DR mechanism, driven by the Renner–Teller coupling of electronic and nuclear motions, is dominant.^4–10 The theory of the indirect mechanism predicts some structure between 0.1 eV and 0.5 eV, drops of σ(ε) (the DR cross section as a function of the colliding energy) connected with the vibrational modes of the ion.^6,9,11 However, such declines are only barely visible in the experiments,^12,13 if at all. Therefore, based on an analytical model, Schneider et al.⁸ supposed that these drops are concealed due to a simultaneous direct DR mechanism.

A possible contribution of the direct mechanism to the DR effect has been investigated in a series of papers that have Å. Larson as a common author;^14–17 the most comprehensive one is ref. 17 by Larson, Stenrup and Orel, henceforth sometimes abbreviated (LStO). In the work of Nordhorn et al.¹⁶ experimental branching ratios of DCO are compared with two-dimensional wave packet model computations for HCO. The results of the (to our knowledge) latest calculation of the total DR cross section,⁹ combining indirect effects with the direct contributions of (LStO), are systematically about 10% to 20% lower than the CRYRING experiment¹³ for 0.02 eV ≤ ε ≤ 0.2 eV.

According to the results of the quantum chemical (multi-reference configuration interaction MRCI) calculations of Larson et al.,^14,15 no resonant state PES crosses the surface of HCO⁺(X¹Σ⁺) near its minimum. Larson and Orel¹⁵ find the lowest vertical excitation energy for a ²Π resonant state to be 2.21 eV relative to their ionic minimum (at r_CH = 2.046a_o, r_CO = 2.099a_o). On the contrary our new MC-CEPA¹⁸ calculations show a ⁴Π and an ion-core-excited (‘resonant’) ²Π surface of HCO with vertical excitation energies only 0.91 eV and 1.72 eV resp. above the same (theoretical) HCO⁺ minimum, whose crossing seams with the HCO⁺ PES are for both multiplicities and both spatial symmetries (A′, A′′) partially located within the range of the ionic ground state vibrational wave function. Motivated by this observation we reexamine the possible significance of the direct mechanism to the DR cross section of HCO⁺. According to our findings the contributions of all open channels add up to a clearly higher cross section of direct DR than in the results of (LStO), who predicted that the energetically high exit channel O + CH(X) is absolutely dominant up to an electron energy of at least 4 eV. We find the limit H + CO(a³Π) (with the lowest possible electronic energy) to be the most probable. Our intention is to determine the dissociation rates for the direct DR processes via the lower core-excited doublet states and to explore the possible contributions of the lowest quartet states (1⁴A′ or 1⁴A′′).

This paper is organized as follows: In Section 2 our quantum chemical calculations of energies and couplings are explained. Some additional details will be given in a separate publication.¹⁹ The applied computational methods of wave packet dynamics (on single as well as coupled doublet and quartet surfaces) are introduced in Section 3, and the resulting DR cross sections are presented and discussed in Section 4. A discussion of the differences of our results compared to ref. 9, 15, and 17 and measurements¹³ is incorporated in Sections 2.4, 4.2, and 5. In Section 6 we introduce population and depopulation rates for long-lived bound vibrational states on the quartet PES, and a concise summary is found in Section 7.

The most recent theoretical papers on HCO, not dealing with DR, focus on the 1²A′ and 1²A′′ PES including vibrational analysis²⁰ and photodissociation.^21,22 The first excited state 1²A′′ of HCO may be populated by absorption of a photon from the bent ground state 1²A′. As the two states are degenerate at linear geometry (Renner–Teller coupling), the subsequent dissociation proceeds on the 1²A′ ground state PES showing an oscillatory pattern of the CO rotational distribution, which has been studied experimentally and theoretically.²¹ In addition to the vibrational bending excitation of the 1²A′′ intermediate state, its rotational excitation and C–H vibrational mode population turns out to be of great importance for the final state CO(ν,j) distribution.²² Usually a molecule (or ion) with the (ground state) PES minimum at non-bent geometry is classified as linear. Winterhoff et al.²³ have calculated the bond angle expectation value for HCO⁺ and DCO⁺. Their theoretical results are compared with Coulomb Explosion Imaging experiments agreeing on non-linear average bond angles for these undoubtedly linear ions.

2. Electronic states and coupling

2.1 Potential energy surfaces

Here, and in the following sections, we presuppose that the electron in an initial (unbound) ²Π state with continuum energy ε recombines with the (linear) molecular ion HCO⁺ in its ground state, where (according to C_s classification) the seven orbitals 1a′ to 6a′ and 1a′′ are doubly occupied. In HCO⁺ 6a′ and 1a′′ are the degenerate components of 1π, whereas in linear HOC⁺ 1π corresponds to 5a′ and 1a′′, and the highest occupied orbital (HOMO) 6a′ is 5σ. The LUMO of the ion is always 2π (or π*; 7a′ and 2a′′). The a′′ component of 2π remains qualitatively unchanged upon bending while the a′ component is increasingly mixed with the valence orbitals of the ion core. But also in bent geometry these two components form a pair of C–O antibonding valence orbitals.

In Table 1 we have collected the experimental energy data of the various adiabatic dissociation limits of some excited HCO surfaces (thereby printing in bold face the numbers of the electronic states discussed here), and also the ground levels of the neutral molecule, the ion, and of the well documented B̃ state. Note that these energies are not thermochemical data, in contrast to e.g. ref. 13 and 24. The main problem here was the position of the ion within the table. In order to calibrate HCO⁺ we combined the value of 8.14 eV for the adiabatic ionization potential IP_a predicted early by Staemmler²⁵ and observed by Dyke,²⁶ with the value of 0.629 eV for the dissociation energy of HCO, converted by Werner et al.²⁷ from data of ref. 28. The table shows that the 2²A′, 2²A′′, 1⁴A′ and 1⁴A′′ surfaces all adiabatically dissociate to the limit C(³P) + OH(X²Π) as well as to O(³P) + CH(X²Π), and three of them to H(1s) + CO(a³Π), whereas 2²A′′ converges adiabatically to the higher limit H(1s) + CO(d³Δ). None of the four doublet or quartet surfaces considered here dissociates adiabatically to the intermediate limit H(1s) + CO(a′³Σ⁺) or to the lowest limit H(1s) + CO(X¹Σ⁺). We see that all these limits are either below the minimum energy of the HCO⁺ ion recombining with a zero energy (ε = 0) free electron or very close above it. The CO interatomic distance of all possible CO fragment states is longer than in HCO⁺ (by about 10% for a³Π and nearly 25% for a′³Σ⁺ or d³Δ). Therefore we expect vibrational excitation upon dissociation to H + CO.

Table 1 Experimental dissociation limits for HCO doublet and quartet states (including zero point energies), and state numbers (in C_s -symmetry) of the converging adiabatic surfaces. The numbers of the states considered in this paper are printed in bold face

Dissociation limits	E (eV)	^2A′	^2A′′	^4A′	^4A′′	E el (eV)
a Energy origin: the separated ground state atoms H, C, and O. b Based on the data of the Huber–Herzberg tables.^29 c Electronic energies with respect to the HCO^+ minimum (zero point energies ZPE subtracted). d Adiabatic ionisation potential IPa = 8.14 eV.^25,26 e ZPE = 3466.1 cm^−1 (0.430 eV);^30 minimum of the HCO^+ PES at −4.01 eV. f T o = 38695.5 cm^−1 (4.798 eV), ZPE = 2521.6 cm^−1 (0.313 eV).^30,31 g ZPE = 2691.7 cm^−1 (0.334 eV),^30D^oo = 0.629 eV.^27,28
O(X^3Pg) + CH(X^2Πr)	−3.465	1, 2, 3	1, 2, 3	1, 2, 3	1, 2, 3	0.38
H(1s) + CO(d^3Δi)	−3.576	4	2	3	2	0.36
HCO^+(X^1Σ^+)^d^e	−3.58					0.00
H(1s) + CO(a′^3Σ^+)	−4.229	3		2		−0.30
C(X^3Pg) + OH(X^2Πi)	−4.392	1, 2, 3	1, 2, 3	1, 2, 3	1, 2, 3	−0.61
H(1s) + CO(a^3Πr)	−5.082	2	1	1	1	−1.18
HCO(B^2A′)^f	−6.923					−3.23
H(1s) + CO(X^1Σ^+)	−11.092	1				−7.21
HCO(X^2A′)^g	−11.721					−8.05

---

In this paper (and in the figures) we apply the appropriate system of Jacobi coordinates r,R,γ depending on the respective dissociation limit, where r denotes the inter-atomic distance in the diatom (CO, OH or CH), R is the distance from the centre of mass of the diatom to the third atom and γ = (r,R) means always a Jacobi angle, while (as in ref. 15 and 17) θ is used for bending angles (θ = 0° for linear arrangement). For the (linear) equilibrium geometry of the ion HCO⁺ we choose γ(H + CO) = γ(O + CH) = 0°, γ(C + OH) = 180°; the equilibrium of HOC⁺ is at γ(C + OH) = 0°. (Of course all possible geometries of the triatom HCO can be described by any of these Jacobi coordinate systems.)

The here studied core-excited states of (neutral) HCO are characterised by a single excitation out of the HOMO of the ion into the (now doubly occupied) LUMO components. In the linear molecule the resulting ⁴Π components are the lowest states of these configurations, followed by a ²Π and a ²Φ state within about 1 eV (and two further ²Π). Complementary to Table 1, Fig. 1 illustrates the paths on the lowest resonant doublet and quartet PES leading to the dissociation limits considered in this study. The lowest core-excited doublet and quartet state PES cross the zero vibration level of the ion close to the minimum of the HCO⁺ surface, as can be seen in Fig. 2a and b. The 1⁴A′ and 1⁴A′′ surfaces, with minima at bent geometries, are strictly adiabatic. Furthermore – as already conjectured by Bruna et al. (see Fig. 3 of ref. 32) – upon bending the molecule the ²A′ and ²A′′ components of the lowest resonant ²Π state are diabatically connected with the second lowest (adiabatic) PES of their symmetry (2²A′ and 2²A′′). With increasing bending angle they fall below the lowest quartet PES.

Fig. 1 Energy correlation and reaction coordinate diagram, energies in eV relative to the minimum of the ground state surface of HCO⁺ (red solid lines). Minimum potential energy paths for the lowest core-excited doublet (blue, purple) and the lowest quartet (green) states as functions of Jacobi coordinate R in the H + CO (blue, green), C + OH (blue) and O + CH (purple) system (solid A′, dashed A′′). All curves start at the equilibrium geometry of HCO⁺ (filled rectangles). The coordinate values r and γ are energy optimised with limited variation between adjacent sampling points in order to avoid irregularities. Vibrational zero point energies (dotted horizontal lines). Some calculated vibrational energies for the quartet (1⁴A′, 1⁴A′′) and the HCO⁺ (X) PES are given in Table 3, doublet vibration (2²A′, 2²A′′) will be discussed in a forthcoming paper¹⁹).

Fig. 2 Potential energy surfaces in H + CO Jacobi coordinates (R = 3.20a_o) for the lowest quartet states (green) and the core-excited doublet states (blue). Energies in eV relative to the minimum of the ground state surface of HCO⁺ (red), (r in a_o and γ in °). The brown planes mark the ZPE of HCO⁺. (a) 2²A′ and 1 ⁴A′′; (b) 2 ²A′′ and 1 ⁴A′. For small values of γ and r the adiabatic surfaces (2 ²A′, 2 ²A′′, yellow) split from the diabatic (core-excited) states. (c) Angular sections through the adiabatic potentials of 2 ²A′ (B̃ state), 2 ²A′′ (C̃), 1⁴A′ (b̃), and 1 ⁴A′′ (c̃) close to their minima, with R = 2.80a_o and r = 2.60a_o.

Fig. 3 Diabatic (blue) and adiabatic curves as functions of r in the linear molecule (a) and of γ for r = 2.06a_o (b) and r = 2.28a_o (c). Red is the potential of HCO⁺ with ZPE as a horizontal dashed line, and green indicates the quartet potentials. In the left part of (a) r_CH ≈ 2.0a_o; in (b) and (c) γ = 40° corresponds to r_CH ≈ 2.4a_o and ∢_HCO ≈ 120°. At γ = 0° r_CH ≈ 2.0a_o in (b), ≈1.9a_o in (c).

The doublet minima have been named B̃ (2²A′) and C̃ (2²A′′) by Tanaka and Davidson.³³ As in the region of the minima the quartets are very similar to the doublets, we propose – violating the rules of nomenclature – to call them b̃ (1⁴A′) and c̃ (1⁴A′′). The characteristics of the minima are given in Table 2, and a section through their almost coinciding potential wells is shown in Fig. 2c (for the quartets see also Fig. 4 and Table 3). The bond lengths in Table 2 are about 10⁻¹a_o shorter than the corresponding results of the CI-based results of ref. 33; our T_e values are up to 1 eV higher than those. Compared to more recent calculations (e.g. ref. 47) the differences of lengths are in the order of 10⁻²a_o, of energies typically a few times 10⁻¹ eV.

Table 2 Calculated characteristics of the minima of the lowest core-excited doublet and quartet surfaces. We obtain the minimum of the ground state surface of HCO at −113.6868 a.u

		r CO (ao)	r CH (ao)	∢HCO (°)	T e (eV)
HCO^+	1^1Σ^+	2.08	2.06	180	8.00
C̃	2^2A′′	2.69	2.08	108	5.04
b̃	1^4A′	2.65	2.05	122	4.94
B̃	2^2A′	2.57	2.08	115	4.90
c̃	1^4A′′	2.62	2.04	118	4.83
X̃	1^2A′	2.24	2.12	126	0.00

---

Fig. 4 Potential energy in eV relative to the HCO⁺ minimum for the lowest ⁴A′ and ⁴A′′ states of HCO as functions of γ and R (in a_o), whereby r is energy-optimized. Jacobi coordinates for the H + CO (a and b) and C + OH (c) dissociation channel. Left ⁴A′, middle and right ⁴A′′. Equilibrium geometry of HCO⁺ at γ = 0° (left and middle), γ = 180° (right). Potential wells b̃ (A′), c̃ (A′′) and c̃′ (A′′), see the text for details.

Table 3 Calculated vibrational energy levels in cm⁻¹ (J = 0) for the lowest quartet states 1⁴A′, 1⁴A′′ and HCO⁺ ground state. The equilibrium and ZPEs are given in eV relative to the HCO⁺ minimum. (Distances in a_o, H + CO Jacobi coordinates, except * C + OH Jacobi coordinates)

E min (eV),γ	−3.060, 34°	−3.177, 36°	0.000, 0°	−2.740, 68°	1.633, 0°
r, R, (label)	2.65, 3.12, b̃	2.62, 3.05, c̃	2.08, 3.25	1.82, 2.60, c̃′*	1.86 2.28*
ν 1 ν 2 ν 3	1^4A′	1^4A′′	HCO^+	1^4A′′	HOC^+
0 0 0	−2.728 eV	−2.840 eV	0.426 eV	−2.353 eV	1.987 eV
0 1 0	1077	1002	846	1133	250
0 0 1	1174	1226	2171	1261	1994
0 2 0	2138	2002	1622	2248	551
0 1 1	2234	2222	2944	2364	2223
1 0 0	2993	3081	3011	3673	3345
0 0 2	2341	2449	4111	2521	3962
0 3 0	3176	2994	2376	3337	880
0 0 3	3493	3646		3775

---

2.2 Method

All calculations of PES points are based on a version of the the multi-configuration reference coupled electron pair approximation method (MC-CEPA) by Fink and Staemmler,¹⁸ modified for applications including electronically excited states. This program processes the input of a multi-configuration reference function of Slater determinants. The orbitals – initially generated from self-consistent field (SCF) calculations (as a rule for the ionic ground state, sometimes directly for a desired configuration) – were improved by iterated medium size CI (configuration interaction) computations, separately for both A′- and A′′-symmetry. The orbitals had an expansion length of about 180 ANO (atomic natural orbitals³⁴) Gaussians, augmented in special cases up to more than 200 and including also elements for the description of low-lying Rydberg or continuum orbitals.³⁵ The CI consisted of roughly 10⁴ configurations, and using transformations to natural orbitals (by diagonalizing a single-particle density matrix) the length of the constructed reference function could in most cases be restricted iteratively to a number between one and approximately a dozen determinants of the most important configurations. In order to avoid stepping down in the CEPA procedure to an unwanted solution at lower energy we have removed certain single de-excitations in the list of single and double excitations generated by this program. As a rule this was sufficient for separating core-excited configurations from the singly excited Rydberg and continuum wave functions. At a few critical points we used an auxiliary wall potential as proposed by Guérout et al.³⁶ (details in Section 2.3) with the effect that the diffuse solutions were shifted to higher energies, while the valence states changed only negligibly. We have computed the PES for R values below 6a_o; for the needs of wave packet propagation the surfaces have then been extrapolated numerically towards their correct (experimental) dissociation potentials (Table 1) up to about R = 11a_o, thereby avoiding the CEPA energy error connected with the dissociation of a triple bond. For the linear minimum of the HCO⁺ PES we obtain −113.3927 a.u. (with frozen and non-correlated inner electron pairs 1s(C) and 1s(O)), to be compared with −113.3562 a.u. in ref. 14 and −113.1347 a.u. in ref. 15 and 17. Our calculated (electronic, adiabatic) ionization potential of HCO is 0.294 a.u. (8.00 eV, 0.05 eV less than the ZPE-corrected experimental value of Table 1).

For the construction of the diabatic core-excited ²A′ and ²A′′ surfaces the reference function was restricted to all Slater determinants with 13 electrons occupying the most stable 7 orbitals, two more electrons in the components of the π* LUMO or its a′ and a′′ C_s counterparts, and iterated as explained above. Rather to the surprise of the authors this ansatz led to a very smooth (quasi)-diabatic transition from the components of the linear core-excited ²Π state embedded in the continuum to the adiabatic 2 ²A′ and 2 ²A′′ PES within R < 4a_o, r < 2.7a_o and γ < 80°, without an appropriate diabatization procedure. Outside of this area we use the adiabatic doublet surfaces. This passage is documented with three one-dimensional sections through the surfaces in the plots of Fig. 3. The diabatic doublet state Ψ_d surfaces (blue in the Fig. 2a and b and 3) coincide everywhere with the strictly adiabatic PES (2 ²A′ and 2 ²A′′) except for a region not too far from the ion equilibrium geometry, where the adiabatic (yellow in Fig. 2, black in Fig. 3) and diabatic PES split. At linear geometry (Fig. 3a) the core-excited doublet and quartet states are degenerate Π states; at r ≈ 2.4a_o the adiabatic 2²Π curve separates from the diabatic state. The (blue) diabatic doublet curve falls below the ion minimum for r ≥ 2.28a_o (at R = 3.20a_o). In C_∞v symmetry the second lowest adiabatic ²A′′ corresponds to the A′′ component of 2²Π, whereas the second ²A′ can be a ²Σ or the A′ component of a ²Π, depending on the order of those states at different geometries (see the crossing of 2²Σ and 2²Π in plot Fig. 3a). The curves in Fig. 3b show the situation for r = 2.06a_o (and R = 3.20a_o): the adiabatic curves 2²A′, 3²A′ and 4²A′ turn into 1²Σ, 2²Σ and 2²Π for γ → 0. Remarkably enough the 3²A′ state features an increasing percentage of the core-excited configuration for γ ≥ 45° (and r = 2.06a_o). At r = 2.06a_o the core-excited A′ state stays above 2²A′, but coincides with the 3²A′ adiabatic state for sufficiently large γ. The seams, where the 2²A′ and 2²A′′ states reach the core-excited configurations, intersect the r = 2.28a_o coordinate plane at γ ≈ 38° and 48° respectively (see curves in Fig. 3c).

In other words the surfaces associated with the core-excited structure are merged into the second adiabatic doublet PES of both symmetries for γ around 40° to 60° (and not too small values of r). Thus the DR process may be described by wave packets starting on the core-excited doublet surfaces at the equilibrium geometry of the ion, and subsequently propagating ‘downhill’ on the diabatic surfaces towards larger r and γ, where the (adiabatic) dissociation limits of the 2²A′ and 2²A′′ PES are within reach. The lowest one is H + CO(a³Π), not accessible for ²A′′ electronic symmetry (see Table 1). To sum up, the wave packets describing the dissociation process on the doublet surfaces pass through three PES regions: Resonant state, diabatic transition zone, and second adiabatic state of C_s-symmetry. As mentioned above, the multi-configuration reference function as input for the MC-CEPA program was restricted to core-excited configurations in the initial region.

The doublet system and the lowest quartet states Ψ_q are coupled by spin–orbit interaction. The dominant one-particle contribution of the Breit–Pauli operator H_SO leads to non-vanishing matrix elements 〈⁴A′|H_SO|²A′′〉 and 〈⁴A′′|H_SO|²A′〉 between configurations of opposite spatial C_s-symmetry,^37,38 and the wave packet may partially pass from the doublet PES over to the quartet. As mentioned in Section 2.1 both doublet surfaces fall below the coupling quartet PES for sufficiently bent geometry. Part of the energetically accessible doublet–quartet transition area is to be seen in the central parts (around γ ≈ 40°) of Fig. 2 and 3. On the other hand the spin–orbit matrix elements coupling directly the initial (doublet) continuum state to one of the dissociative quartet surfaces are expected to be very small due to the fact that the involved Slater determinants differ in more than one spatial orbital.

On the lowest ⁴A′ and ⁴A′′ PES the dissociation channels to H(1s) + CO(a³Π) as well as to C(³P) + OH(X²Π) are open for the considered initial energy range, as is documented in Table 1 and Fig. 4. The dissociation path leading to the higher limit O(³P) + CH(X²Π) after spin–orbit transition to one of the quartet surfaces was not incorporated in the present study. The 1⁴A′ surface has only one minimum (b̃), located at bent geometry on the HCO side (see Table 3; the groove in the left part of plot 4a shows the exit channel to C + OH with r → ∞), the ⁴A′′ PES has a HCO (c̃) as well as a HOC minimum (c̃′, see Section 6 for further details). The saddle point between these two minima is near −1.0 eV and r_CO ≈ 2.60a_o, r_CH ≈ 2.35a_o, and r_OH ≈ 2.30a_o. The relative strength of the various coupling mechanisms is discussed in the following subsection.

All involved PES were calculated on a regular grid in Jacobi coordinates (15 angles γ, Δr = 0.12a_o, ΔR = 0.2a_o) and interpolated by three-dimensional spline functions. The density of the grid was enhanced in various critical regions (for example around the local minima).

2.3 Electronic coupling matrix elements

Ψ^Λ_c(ε) denotes the components of the initial antisymmetrized electronic wave function of the neutral unbound system (HCO⁺ + electron with energy ε in a continuum orbital χ^Λ_π), where Λ abbreviates A′ or A′′. Ψ^Λ_d means the symmetry components ( or ) of the core-excited ²Π state, and Φ_o is the vibrational ground state wave function of the ion HCO⁺. The transition from the initial state Ψ^Λ_c to Ψ^Λ_d is mediated through the electronic Hamiltonian H_el,³⁹ where the coupling matrix element

V^cdΛ_ε(r,R,γ) = 〈Ψ^Λ_c(ε)|H_el|Ψ^Λ_d〉_e, (1)

is used to map Φ_o onto the initial nuclear (energy normalized) wave packet Ω(t = 0):

Ω^Λ_ε(r, R, γ, t = 0) = V^cdΛ_ε(r, R, γ)·Φ_o(r, R, γ). (2)

(The subscripts〉_e or 〉_n at the closing brackets indicate whether electronic or nuclear coordinates are to be integrated over; vibrational energies of the ions HCO⁺ and HOC⁺ are documented in Table 3). V^cdΛ_ε(r, R, γ) measures the interaction of the continuum Ψ_c with the core-excited doublets. It also determines the (geometry dependent) autoionization width Γ of the Ψ^Λ_d states:⁴⁰

Γ^Λ_ε(r, R, γ) = 2π(V^cdΛ_ε(r, R, γ))². (3)

Due to energy normalization of the continuum orbital χ^Λ_π entering the calculation of V^cdΛ_ε, Γ carries the dimension of an energy.

As the continuum wave functions Ψ^Λ_c(ε) contribute to V^cdΛ in eqn (1) only in a restricted radial and angular range around the equilibrium conformation of the ion (defined by the size of its occupied orbitals), they may be described within this region by conventional quantum chemical methods, using an auxiliary wall potential (AWP) as proposed by Guérout et al.,³⁶ in order to obtain a locally correct description of χ^Λ_π. The applied cubic box, centred at the positive charge centre of the HCO⁺ ion, with half length ρ = 11a_o in the three Cartesian coordinates ξ_j, zero potential W = 0 inside and slope was chosen to be large enough to not affect the slender lower bound electronic states. The box parameter b (the ‘steepness’ of the wall) then only controls the shape of the (pseudo-)continuum orbital χ^Λ_π within the box and allows the variation of its energy. A linear combination of regular and irregular Coulomb functions f_l and g_l⁴¹ (with angular momentum l = 1) was fitted to the radial part of the quantum chemically calculated continuum orbitals in order to determine the energy ε of the captured electron and to energy normalize the wave functions. An additional angular momentum analysis at linear configurations confirmed the pπ character of the continuum orbitals.

V ^cdΛ_ε was then calculated numerically at the CI level on a three-dimensional grid of (nuclear) geometries around the equilibrium conformation of the ion, applying two different orbital sets for the ε-adjusted continuum state (b > 0) and the core-excited doublet state (b = 0). The selected p-orbitals were shifted to 0.02 < ε < 0.8 eV (and in some cases 1.0 eV) by box parameters b in the range 2 × 10⁻⁸ to 9 × 10⁻⁷ a.u. a_o⁻⁸. The evaluation of the matrix elements in eqn (1) involved the transformation between the two sets of orbitals and consideration of contributions from the one- and two-particle operators (i.e. one or two occupation differences of the contributing Slater determinants). As a rule increasing r or ε lets the interaction matrix element decline, whereas V^cdΛ_ε is slightly enlarged with increasing R. For low energies it is nearly ε independent; at the equilibrium of HCO⁺ and ε = 0.1 eV we obtain Γ = 564 cm⁻¹. But it was found that the characteristics of V^cdΛ_ε(r, R, γ) are symmetry dependent: bending the molecule significantly increases (decreases) its value for Λ = A′ (A′′). This becomes evident if alternatively the coupling strength of the core-excited doublet states to initial continuum is characterized by the (geometry independent) expectation values

V̄_ε,Λ = 〈Φ_o(r,R,γ)²·V^cdΛ_ε(r,R,γ)〉_n, (4)

averaged over the squared vibrational ground state wave function Φ_o and separating results for A′ and A′′ symmetry. For example at ε = 0.1 eV we obtain while is significantly smaller due to the opposed γ-dependence of the V^cdΛ.

As mentioned in Section 2.2 the Breit–Pauli spin–orbit matrix elements V^cqΛ = 〈Ψ_c(ε)|H_SO|Ψ^Λ_q〉_e, coupling electronic continuum directly to the lowest quartet states, are for both symmetries much smaller than V^cdΛ, thus their contribution to the DR cross section is neglected.

2.4 Comparison with previous calculations

Larson and Orel¹⁵ show doublet potential curves for the core-excited (resonant) states as functions of r_CO, r_CH and the bending angle θ. While their r_CO dependence agrees well with our curve for the lowest resonant ²Π state – compare the position of the core-excited doublet and quartet curves in Fig. 3b and c – the other plots exhibit distinct differences: The lowest ²A′ and ²A′′ potentials in ref. 15 are practically constant as functions of θ (for r_CO = 2.08a_o, r_CH = 2.0a_o), whereas our corresponding curves decrease monotonically up to γ≈50° and cross the HCO⁺ potential at γ ≈ 24° (θ ≈ 38°) and γ ≈ 33° (θ ≈ 52°) respectively (Fig. 3b). And as a function of r_CH our ²Π curve falls partly below the ion surface first for r_CO ≥ 2.22a_o (Fig. 3c), in contrast to in ref. 15 where the lowest resonant curve intersects the ion PES at r_CH ≈ 3.3a_o (for r_CO = 2.08a_o and θ = 0°). In ref. 17 a rather crude method for obtaining three-dimensional surfaces is introduced, neglecting the fact that r_CO increases on dissociation to H + CO. The dissociation limits and the electronic states of the fragments are not specified. Our results for the resonance state are lower in energy than those of (LStO), i.e. they cross the ion PES closer to its minimum.

The descriptions of the applied quasidiabatization methods here and in ref. 15 sound quite similar, but they lead to different results. It seems that the quasidiabatized surfaces of (LStO) are not connected with the 2²A′ (B̃) or 2²A′′ (C̃) states of HCO; in any case these core-excited states are never mentioned in their texts. The differences are probably due to the unlike angular comportment of the resonance state PES.

In contrast to our AWP approach (LStO) use the complex Kohn variational method^42,43 to calculate the autoionization width Γ. They detect a considerable increase of Γ at geometries, where the ionic and resonant state PES come close. For the lowest ²Π resonant state at linear geometry with r_CH = 2.0a_o, r_CO = 2.08a_o they obtain Γ = 290 cm⁻¹, about half our value mentioned above. But due to technical limitations concerning their program the θ-dependence of Γ was disregarded in ref. 17, which might lead to an underestimation of the (direct) DR cross section as result of too small a contribution of the core-excited ²A′ state (and by comparison too high for ²A′′).

3. Wave packet dynamics

This section is dedicated to the construction and propagation of wave packets of nuclear motion moving on the diabatized PES (introduced in Section 2.2), from the geometry of electron capture to the various dissociation limits of the adiabatic 2²A′ and 2²A′′ surfaces (in Section 3.1) or, after a doublet–quartet transition, to the dissociation limits of the lowest quartet surfaces (in Section 3.2). From the results DR cross sections as functions of the electron energy ε and also the (H + CO)/(C + OH) branching ratio (in Section 4.1), as well as the occupation of vibrational levels of the CO(a³Π) fragment (in Section 4.3) will be derived.

3.1 Wave packet propagation on a single core-excited doublet surface

We start the exposition of our wave packet calculations with the case of motion on a single PES U(r,R,γ). The propagating wave packet Ω^JM(r,R,γ,t) with rotational quantum numbers JM is expanded in an (orthonormal) Jacobi coordinate product basis set of the diatomic fragment (r) and the dissociation channel (R)with coupling Clebsch–Gordan coefficients in the angular part. The radial one-dimensional sets u_α(r) and v_β(R) solving the time-independent Schrödinger equation for suited potentials, are linear combinations of spatially distributed Gaussian functions in r and R respectively (i covers all occurring combinations of α, β, j and l). Time evolution of the wave packet on a calculated PES U(r,R,γ) is generated by the Hamiltonian H_U = T_N + U − iπV², where T_N denotes the kinetic energy operator of nuclear motion and the interaction potential V is defined in eqn (1). The effect of the autoionization term −iπV², which was included as first order perturbation, turned out be largely negligible. In other words recombination leads almost always to dissociation. All wave packet calculations in this paper were restricted to J = 0.

We may write the solution

of the time-dependent Schrödinger equation in the usual way, where the initial wave packet Ω^Λ_ε(t = 0) is given in eqn (2). E is the energy of the initial wave packet relative to the asymptotic dissociation limit, thus exp(−iEt′/ħ) in eqn (6) ensures conservation of energy.^44,45 The relaxation time τ controlling the propagation start was chosen to be τ = 8 fs to τ = 10 fs, whereas the typical time t_e of cross section evaluation was 26 fs ≤ t_e ≤ 33 fs (depending on ε and the final state). τ has to cover multiple periods (h/E) of oscillation but should be clearly smaller than time t_e of the final wave packet analysis. With the chosen (channel) basis sets v_β(R) we can treat initial energies ε from 0.02 eV up to about 0.8 or 1 eV. Expanding the angular dependence of the potential U(r,R,γ) in Legendre polynomials P_λ(cos(γ)) and applying the Spherical Harmonic Addition Theorem simplifies the calculation of the matrix elements H_ik = 〈η_i|H_U|η_k〉_n. In order to avoid expensive diagonalization of the Hamiltonian matrix H, typically of dimension 35 000, the Lanczos algorithm⁴⁶ was applied (as in our previous study on predissociation of H₃⁴⁵) to transform it into tridiagonal form in a new basis covering the range up to R ≈ 11a_o, in which wave packet dynamics could be carried out efficiently. The Lanczos transformation is the most time consuming step of the procedure, but it is performed only once and greatly accelerates the propagation step for any chosen energy E. Eqn (6) was integrated numerically, using a fourth order expansion of the integrand and Δt = 10⁻³ fs time steps.

Unfortunately the employed Jacobi coordinates (r,R,γ), where r asymptotically is the bond length of the respective diatomic product, are not well suited for the simultaneous description of multiple dissociation limits (H + CO, C + OH, O + CH). Therefore the branching ratios could only be estimated by carrying out separate wave packet calculations in every Jacobi coordinate system. As dissociation to the respective other fragmentation channels corresponds to r → ∞, it is important to use a sufficiently extended r-basis.

Optionally an alternating propagation technique was introduced as an attempt to improve interaction of the two most important ²A′ channels in the early stages of wave packet motion. During the first 16 fs propagation was alternated every fs between the two channels H + CO and C + OH, transforming the current wave packet to the respective other Jacobi coordinate system. This periodic superposition of the wave packets, calculated in different Jacobi coordinate systems, could be restricted to the (geometrically) inner part of non-vanishing overlap of the two channels. It was damped with the time constant τ (of eqn (6)), implicating the possibility of initial amplitude exchange between channels. This procedure also increases the number of necessary integration time steps (not, of course, the real time to dissociation). The effect of the alternating propagation approach on branching ratio and cross sections will be discussed in Section 4.1.

3.2 Wave packets on coupled (doublet and quartet) surfaces

As stated in Section 2.3 the direct recombination process typically leads to a population of the core-excited doublet states. At linear configuration the matrix elements of spin–orbit interaction 〈H^a_SO〉 = 〈⁴A′′|H_SO|²A′〉 and 〈H^b_SO〉 = 〈⁴A′|H_SO|²A′′〉 are both approximately 50 cm⁻¹; 〈H^a_SO〉 remains practically constant as a function of the Jacobi angle (see also ref. 47), whereas 〈H^b_SO〉 is distinctly smaller in the vicinity of the intersection of the corresponding PES (⁴A′, ²A′′). The pairs of surfaces to be coupled are labeled [U^a_d(²A′), U^a_q(⁴A′′), and U^b_d(²A′′), U^b_q(⁴A′)], respectively. In order to carry out coupled wave packet calculations we follow the method introduced by Mead and Truhlar.⁴⁸ For every point (r,R,γ) of the PES a 2 × 2 matrix of the elements (Uⁱ_d and Uⁱ_q) of the two interacting states on the diagonal and 〈Hⁱ_SO〉 as coupling element is diagonalized; the mixing ratio is interpreted as sin ϑ_i and cos ϑ_i of an additional (geometry-dependent) mixing angle ϑ_i(r,R,γ) (with i = a or b). In a second step the initial elements of Uⁱ_d, Uⁱ_q and 〈Hⁱ_SO〉 are replaced by

Uⁱ₁₁(r,R,γ) = Uⁱ_d cos²ϑ_i + Uⁱ_q sin²ϑ_i (7)

Uⁱ₂₂(r,R,γ) = Uⁱ_d sin²ϑ_i + Uⁱ_q cos²ϑ_i (8)

Uⁱ₁₂(r,R,γ) = (Uⁱ_d − Uⁱ_q) sin ϑ_i cos ϑ_i. (9)

These new coupled (double-layered) PES lead to Hamiltonian matrices H of double dimension including the possible doublet (d) to quartet (q) transition, and thus open for the wave packets the path to dissociation on one of the quartet PES in both cases (a,b). The operator of kinetic energy in eqn (6) does not change but the basis set (5) expansion length of the wave packet is doubled compared to propagation on a single surface. The Lanczos algorithm was applied to these full Hamiltonian matrices composed of the blocks ((7)–(9)) allowing wave packet flux between states of unlike multiplicity. This formalism will be used in Sections 4.1 and 6.

4. Direct DR cross sections

4.1 Method and results

Regardless of whether time evolution of the wave function Ω_ε(t) is restricted to a doublet PES or whether the quartet component of a coupled propagation is considered, the cross section for dissociative recombination may be calculated via projection of the asymptotic wave packet onto the (energy normalized) final statesHere ϕ_jν(r) describes the rovibrational states jν of the diatomic fragment, μ is the reduced mass, and k_E is a result of the final translational energy E. The cross section σ is given by ref. 49,50where propagation time t_e must be chosen large enough (see Section 3.1) so that the dissociating part of the wave packet has left the interaction zone, but is still within the range of the radial basis v_β(R). k is the wave number of the captured electron and g_i, g_f are the multiplicities of the initial and final state respectively. According to eqn ((1)–(3)) the autoionization width Γ^Λ_ε(r,R,γ) is included linearly in the cross section through the initial wave packet Ω^Λ_ε(r,R,γ,t = 0). In order to ensure that only the asymptotic part of the (complex) wave function enters evaluation, the part with R < 6a_o was suppressed by convolution with a function c(R) prior to integration over R in eqn (11). A (real) damping function at the outermost end of the R-basis was found to be more effective to prevent disturbing wave packet reflection at the end of the basis set v_β(R) than application of an imaginary absorbing potential.

Fig. 5 shows the computed DR cross sections σ as functions of the initial energy ε of the free electron. Our computations are limited to energies between 0.02 eV and 1 eV, because energy control of the incoming electron with too low energy is difficult (see Section 2.3) and because the use of a basis (5) suited for higher energies renders the Lanczos procedure clumsy and slow. For both symmetries (A′, A′′) the sums of the squared L²-norms of the asymptotic parts of the wave packets were calculated as a quality check. In the energy range 0 < ε < 0.9 eV the deviations of these sums from unity nowhere exceeded 0.1, and only slightly between 0.9 eV and 1.0 eV. Thus the applied technique of separately calculating the contributions to the DR cross section for every exit channel appears to be justified. In particular we may conclude that the percentage of the wave packet that remains in the vicinity of initial geometry, when treated within one of the three Jacobi coordinate systems, represents the share, which dissociates to the other two exit channels. Furthermore, the convolution function c(R) seems to be reasonably chosen. The remaining (minor) total normalization error was corrected for with energy-dependent factors (individually for both symmetries) resulting in the dashed curves in Fig. 5. As to be expected, propagation on the dissociative (diabatic) ²A′ PES is clearly the most important of all considered (direct) reaction mechanisms, and the limit H + CO (blue crosses) is more probable than C + OH (red stars) or O + CH (cyan disks). At ε = 0.15 eV the cross section ratio is σ_H+CO/σ_C+OH ≈ 2.7. Dissociation to H + CO(a) is also supported by the experiments by Adams and Babcock,⁵¹ and by Rosati et al.²⁴ who detected CO(a³Π_r) in a plasma producing HCO⁺ that subsequently reacts with free electrons.

Fig. 5 Calculated cross sections for DR of HCO⁺ as functions of electron energy ε. Upper plot: Direct dissociation via the core-excited ²A′ state (→H + CO(a³Π) blue crosses, → C + OH(X) red stars) and ²A′′ state (→ C + OH(X) green triangles). The black curves include all considered mechanisms. Dashed curves: Separate propagation in each channel. Solid curves: Improved propagation (alternating propagation, see the text). Lower plot: direct dissociation via the core-excited ²A′ state (→ O + CH(X) cyan disks) and ²A′′ state (→ H + CO(d³Δ) brown squares, → O + CH(X) orange circles). spin–orbit transition ²A′ → ⁴A′′ (filled diamonds) and ²A′′ → ⁴A′ (empty diamonds) leading to H + CO(a³Π).

Typically the wave packet fraction transferred from the ²A′ to the ⁴A′′ surface is of the order 10⁻⁴ and its contribution (marked with black diamonds in the lower part of Fig. 5) to the cross section is very small. As mentioned in Section 2.2 the spin–orbit matrix elements coupling the initial electronic doublet (continuum) state directly to the quartet surfaces is negligibly small and is not incorporated in Fig. 5.

The solid curves in the upper plot of Fig. 5 show the effect of the alternating propagation technique between the two channels H + CO and C + OH (on the ²A′ PES, as explained at the end of Section 3.1), which led to an enlargement of the branching ratio (σ_H+CO/σ_C+OH) to about 5.3 at ε = 0.15 eV and to a slightly larger total direct cross section. And the total number of time steps Δt to reach the asymptotic region of evaluation was increased to about 50 000, without affecting the final translational energies in both channels. The total normalization correction was carried out analogously as outlined above. The increase of the branching ratio is associated with a weaker (final state) vibrational CO excitation (see Section 4.3), hence channel interaction during the prolonged initial propagation phase might facilitate vibrational relaxation of the wave packet.

4.2 Discussion

For a number of reasons the recombination process populating the ²A′′ state yields for ε < 0.5 eV (green, brown and orange curves in Fig. 5) a significantly smaller cross section than ²A′ (blue, red): The coupling matrix element to continuum is at bent geometries smaller than for ²A′, the ²A′′ PES has a higher H + CO(d³Δ) dissociation limit, and the intersection with the ionic PES is located somewhat less favorable (see Fig. 2a and b). Moreover the spin–orbit matrix element mediating the coupling to the (dissociative) ⁴A′ state changes sign near the doublet–quartet intersection seam, thus this quartet contribution to the cross section is even smaller than for ²A′ → ⁴A′′. And the low limit C + OH, accessible also for ²A′′ symmetry, requires breaking the CO bond.

The electronic energies of the dissociation limits H + CO(d³Δ) and O + CH(X) correspond quite precisely to the initial level of HCO⁺ including its ZPE, but are a little higher. Asymptotically at least the vibrational ground state of the diatom has to be accessible, hence the contributions of H + CO(d³Δ) and O + CH(X) to the DR cross section vanish for ε → 0. As the vibrational frequency ω_e = 1172 cm⁻¹ for CO(d³Δ) is considerably lower than ω_e = 2859 cm⁻¹ for CH(X), the consequent decrease of σ(ε) in Fig. 5 is shifted to higher energies ε in the latter case (compare the low energy end of the cyan and orange curves in Fig. 5b). The weak structure in the O + CH curves may be interpreted as follows: Around ε ≈ 0.2 eV only ν_CH = 0 is within reach, whereas at ε ≈ 0.55 eV ν_CH = 1 begins to contribute.

According to the computations of Larson et al. (LStO)¹⁷ direct DR of HCO⁺ proceeds (for ε < 2.5 eV) mainly via the same (lowest) ²Π resonant state as is considered in our present study. But in contradiction to our results they predict the O + CH channel to be absolutely dominant for ε < 4 eV, and their total direct cross section is clearly lower than ours up to ε ≈ 0.6 eV (see Fig. 7 in Section 5). For example they determine the direct DR cross section to be about 1.4 × 10⁻¹⁶ cm² at ε = 0.15 eV, whereas our value at the same energy is 5.01 × 10⁻¹⁶ cm² (or even 5.45 × 10⁻¹⁶ cm² with the alternating propagation method). This discrepancy may be explained for one thing by the less advantageous intersection of their doublet (resonant) surfaces with the ionic PES, but more importantly by their weaker coupling V^cd (eqn (3)) to continuum.¹⁵ In ref. 17 the A′ and A′′ share turn out to be approximately the same size, which may be explained by the neglected dependence of Γ upon bending (see Section 2.3). In our calculations the A′ contribution is distinctly augmented by the increase of width Γ as a function of γ. According to ref. 7 (LStO) compute the direct to indirect effect ratio to be 17% at ε = 0.13 eV. And they find a minimum of the cross section resulting from the lowest ²Π resonant state at ε ≈ 0.5 eV, where our σ value is about twice as large and – in contrast to (LStO) – our total direct σ(ε) is monotonically decreasing for all examined electron energies and arrives at 10⁻¹⁷ cm² near 1.0 eV.

Interestingly (LStO) state the analogous dominance of the DC + O channel for direct DR of DCO⁺ as in the case of HCO⁺. But with a simple Landau-Zener⁵² model allowing coupling and flux redistribution to Rydberg states they find the low energy branching ratio to completely change to approximately 85% D + CO (in accordance with the experiments at the Heidelberg TSR¹⁶). Disregarding any coupling to Rydberg states we obtain 72% H + CO(a³Π) at ε = 0.1 eV for the contribution of the direct mechanism (solid blue curve in Fig. 5).

Higher resonant states might turn out to be important for direct DR with increasing ε. In preliminary calculations we have found the lowest ²Φ to be the most promising additional state. We obtain about 6000 cm⁻¹ as the vertical energy difference between ²Π and ²Φ at the equilibrium geometry of HCO⁺ and ²Φ crosses the ionic surface about 1500 cm⁻¹ above the ZPE of the ion.

4.3 Vibrational analysis of CO(a)

For ²A′ → H + CO(a³Π) a final state vibrational analysis was carried out by projecting the r dependence of the (asymptotic) wave packet onto a set of vibrational CO wave functions. As may be seen in Fig. 6 the population of the six lowest vibrational states is quite stable as a function of ε (though the results are rather uncertain for levels higher than ν = 5 due to the limited dimension of the r-basis). In the initial transition region the core-excited ²A′ PES features a pronounced slope towards larger r (see Fig. 2a), which might in combination with the excess energy explain why even for low ε some excited vibrational states are notably populated. Except for the population of the level ν = 0 our results are in qualitative agreement with the Rosati experiment²⁴ insofar as for the levels 1 to 4 approximately equal (not thermal) populations are predicted. For the population of the ground level however Rosati finds 45%, more than twice our value. But then our calculations do not exactly correspond to the setup of this experiment, e.g. the high vibrational ground state population might also be influenced by collisional relaxation in the flowing plasm. And of course we only cover the directly coupled contribution of the DR processes. Our restriction to total angular momentum J = 0 implies l = j (see eqn (5)). The limited diatom rotational quantum number (j ≤ 19) of the wave packet expansion might also have some influence on the calculated final state vibrational distribution.

Fig. 6 Calculated final state vibrational CO(a³Π) population as a function of ε for direct dissociation via the core-excited ²A′ state: ν = 0 black crosses, ν = 1 blue stars, ν = 2 green triangles, ν = 3 red disks, ν = 4 brown diamonds, ν = 5 gray squares, ν = 6 purple circles.

5. Comparison with experiment and with calculations including the indirect mechanism

There exist two series of experiments studying DR of HCO⁺: Merged beam studies by Padellec et al. from 1997 (Pa)¹² and extensive CRYRING storage ring experiments by Hamberg et al. (Ha),¹³ published in 2014. In the energy region up to 1.5 eV the latter date from 2006 and have been reproduced and cited as ‘unpublished’ by several authors. However, above ε ≈ 0.5 eV (where part of the signal is interpreted as background) the published DR data of (Ha) do not exactly correspond to their earlier version (in particular above 0.4 eV the ‘unpublished’ data points are below the error limits postulated in ref. 13). In both studies a power function decrease of the DR cross section σ is reported, between 0.01 eV and 0.3 eV as 1.0 × ε^−1.7 cm² (Pa) or 1.3 × ε^−1.3 cm² (Ha). At energies higher than ε ≈ 0.7 eV the results of both experiments become irregular, probably because of the limitation of their sensitivity for σ < 10⁻¹⁶ cm². In both papers it is explicitly stated that no DR effect could be detected above 1.0 eV (but (Ha) find resonance structures between 3 and 15 eV). At energy below 0.3 eV the cross sections measured by (Ha) are higher than those from (Pa) by a factor between 2 and 4; (Pa) report a broad peak between 0.3 and 0.8 eV, which however is absent in the results of (Ha). Based on quantum chemical calculations (Pa) believe that after electron capture the molecules dissociate to ground state H + CO, but they raise the question “how it is that CO(a³Π) molecules were seen as recombination products in the experiment of Adams and Babcock”.⁵¹ In addition a study by Nordhorn et al.¹⁶ from 2011 contributes branching ratios for the DR of DCO⁺ between 0 and 25 eV collisional energy (but no cross sections). (Ha) found that almost 90% of the DR events of H¹³ CO produce H + CO (not distinguishing ground and electronically excited CO fragments), plus rather more C + OH than O + CH.

Most theoretical studies of the indirect DR process predict cross sections roughly proportional to ε⁻¹ (caused by the k⁻² prefactor in the expressions for σ) up to ε ≈ 0.1 eV and higher^4–7,10 with slightly differing amplitudes. These methods only use the (spectroscopic or calculated) properties of the free HCO⁺ ion, and so produce smooth curves for σ(ε). The main difference lies in the value of ε (between 0.1 and 0.5 eV depending on the choice of internal degrees of freedom taken into account), where the predicted indirect cross sections fall to low values. Such drops are expected if the collision energy opens a new vibrational channel of the ion,^8,11 but they could not be observed unequivocally for HCO⁺, where they would appear near 0.1 eV (ω₂), 0.3 eV (ω₃) or 0.4 eV (ω₁). Therefore Schneider et al.⁸ ask “whether a direct process might also contribute” (and hence conceal the drop). It seems that an attempt to describe the charge distribution of the ion with higher electric multipoles⁷ has been given up, perhaps because a dipole field has been used instead of a quadrupole placed at the center of (positive) charge. None of the cited studies of the indirect process makes predictions about branching ratios (except assuming that dissociation leads to ground state H and CO fragments). Previous theoretical studies of the direct process culminate in the paper by Larson et al. (LStO),¹⁷ introduced in the foregoing sections.

Let us focus on the two most recent theoretical papers from 2014⁹ and 2015,¹⁰ written by the same authors and using similar methods (the low-energy scattering matrix calculated through the variational complex Kohn method^42,43), but arriving at quite different results. In the 2014 paper only the total DR cross section is documented, but it is stated that “the contribution of the direct process is only significant at energy above ∼0.1 eV”, from where on the direct (LStO) results are used (and thus the indirect part can be reconstructed). Between 0.001 eV and ∼0.07 eV the 2014 cross sections for the indirect process are about 2.5 times higher than the 2015 values. At higher ε this factor quickly increases (because the 2015 results drop at ω₂, those of 2014 at ω₁) up to about 0.14 eV, where the 2015 data – restricted to the indirect process – end to be documented. In both papers it is assumed that the experiments have been carried out with a mixture of 90% HCO⁺ and 10% HOC⁺, but the difference to calculations assuming pure HCO⁺ vanishes above 0.05 eV. Between ε = 0.005 eV and 0.1 eV the indirect DR contribution of 2015, obtained for the gas mixture, is quite close to the (Pa) experiments.

The total cross sections of the 2014 paper are the only ones that come fairly close to the CRYRING results. Up to ε ≈ 0.4 eV this is mainly due to the contribution of the indirect process, roughly parallel to the experiments and attaining about 75% of their value (below 0.3 eV proportional to ε^−1.4). Near 0.4 eV the drop of the indirect DR process sets in, due to its interaction with the highest vibrational mode of HCO⁺, and the cross section becomes dominated by direct DR. Combination with the (LStO) results for the direct process leads to a step downwards and then to continuation with some structure at a relatively constant level just below σ ≈ 10⁻¹⁶ cm² up to 1.5 eV, in contrast to the experiments decreasing to 10⁻¹⁷ cm² within the same interval.

If however the calculations for the indirect process are combined with our present results for direct DR the picture changes. Fig. 7 compares our calculated cross section results for the direct process with theoretical and experimental curves reconstructed graphically from Fig. 6 of ref. 9. At ε ≈ 0.1 eV the sum of the 2014 calculations for the indirect DR mechanism and our cross section for the direct process matches well with the (Ha) experiment. Between 0.2 and 0.4 eV this sum is even higher than the (Ha) data (perhaps because the drop of the indirect process near ω₃ has been neglected), but it never exceeds the upper error limit given in ref. 13. Up to ε ≈ 0.6 eV our direct mechanism cross section is distinctly larger than according to (LStO), as discussed in Section 4.2. Above 0.6 eV, in contrast to both the (LStO) calculations and the ‘unpublished’ experiments, only our regularly decreasing results for the direct (and thus total) DR match with the statement by (Ha) that “the number of DR events around 1–3 eV is close to zero”. Maybe in consequence of neglected higher resonance states our calculated curve falls somewhat too fast for ε > 0.6 eV (see the end of Section 4.2). Rather than a drop of the indirect contribution, above 0.6 eV the shape of our total DR reflects the decline of the direct process.

Fig. 7 Measured and calculated (total, direct, indirect) cross sections for the DR of HCO⁺ as functions of electron energy ε (adapted from Fig. 6 of ref. 9). Blue: CRYRING experiments;¹³ Red solid line: Fonseca et al. (ref. 9, total cross section); brown dot-dashed (and red solid line for ε > 0.6 eV): Larson et al. (ref. 17, direct process); green dashed: Fonseca et al. (ref. 9, indirect process, difference of the red and brown curves). Black: dashed line (and black solid line for ε > 0.6 eV): our present results (direct process); black solid line: sum of green and black dashed curves (total DR cross section).

6. Population and depopulation on quartet pes

Dissociation on the quartet surfaces hardly contributes to the DR cross section (see Section 2.2). Apart from dissociation spin–orbit coupling allows for population of the vibrational states in the potential wells on the lowest quartet PES. As the lower ones of these levels – in particular those of the ⁴A′′ state – are situated below the lowest quartet dissociation limit, further deactivation or dissociation is possible only after a second spin changing transition. Therefore these levels – if occupied – are expected to be long lived. In the present section we shall have a closer look at these processes.

The 1⁴A′′ PES shows two minima at bent geometry (see Fig. 4), whereas there is but one pronounced potential well on the 1⁴A′ surface. In Table 3 our calculated vibrational energies are given for the mentioned wells (b̃, c̃, c̃′, see Section 2.1) on the quartet surfaces and for the two linear minima of the ion ground state PES. The same type of basis set expansion (see eqn (5)) as for wave packet dynamics was applied to determine vibrational frequencies. For the three quartet minima we observe a rather strong mixing of the angular and radial vibrational modes ν₂ and ν₃. Our vibrational energies for the b̃ and c̃ minima are quite consistent with the theoretical (harmonic) frequencies found by Manaa and Yarkony.⁴⁷ The same applies to the frequencies of the ions HCO⁺ and HOC⁺, if compared with the values given in the Jacox tables.³⁰ To our knowledge for HOC⁺ only [small nu, Greek, tilde]₁ has been measured, but our fundamental frequencies are also in reasonable agreement with the theoretical results of Mladenović et al.⁵³ (predicting an unusually flat bending potential). The corresponding vibrational energies for the doublet PES and further details of the vibrational calculations will be given in a coming publication.¹⁹

A possible mechanism to release the excess energy of the quartet wave packet (introduced in Section 3.2) and populate a low-lying vibrational state Φ_vib is a dipole transition within the respective quartet state. In order to evaluate the transition dipole operator D⃑ between a wave packet and vibrational functions, accounting for nuclear as well as electronic contributions, the electronic dipole transition moment function was calculated at the CI level on a limited grid in the (geometrical) vicinity of the considered potential well. With Einstein's factor A of spontaneous emission a time-dependent rate is computed, which was integrated up to 36 fs, a characteristic timeframe after which the dissociation process on our coupled surfaces is essentially terminated. The wave packet energy normalization factor was omitted here, thus only relative populations of the individual vibrational levels after a (single) recombination event result. They are in the order of 10⁻¹⁵, depending only weakly on energy or level. Populations of the lowest five vibrational states for the two minima on the 1⁴A′′ PES are shown in part (A) of Fig. 8 as functions of the captured electrons energy ε. It seems to be rather plausible that well c̃ at γ ≈ 37° is more likely reached than well c̃′ on the HOC side. Only a tiny fraction of the initial wave packet on the doublet surface reaches the quartet state, where dissociation is the absolutely dominating process. This explains the very low probability for population of a bound vibrational state. On the other hand the lifetime of these vibrational levels is expected to be very long, provided that solely radiative processes of depopulation are to be considered.

Fig. 8 (A) Relative population of the lowest five vibrational levels (after a single recombination event of HCO⁺) for minima c̃ and c̃′ on the 1⁴A′′ PES. (B) Depopulation rates to H(1s) + CO(X¹Σ⁺) for potential well c̃ (1⁴A′′). For the vibrational ground state (ν₁ν₂ν₃) = (000) the specific contributions of individual final state vibrational levels ν_CO are shown by dashed lines.

Our calculated depopulation rates per unit energy interval for the five lowest vibrational levels of system c̃, w(E_k, ν₁ν₂ν₃ → CO(X)), are depicted in part (B) of Fig. 8 as functions of the fragments translational energy E_k. The spin–orbit operator H_SO leads for both the initial electronic state as well as the final state to first order expansion terms of interchanged multiplicity and spatial symmetry:³⁸

where the coefficients a_q, a_X, normalizing the electronic wave functions obtained with perturbation theory, are very close to unity. In both sums only open shell configurations with three (or more) singly occupied orbitals contribute. Thus the resulting (electronic) first order transition operator is composed of two types of terms containing dipole matrix elements either between doublet or between quartet states respectively.

In order to analyze the translational and vibrational energy distribution of the resulting CO(X¹Σ⁺) fragments, energy normalized Gaussian wave packets Ω_ν(E_k,t) were started on the surface backwards (and upwards) from the CO levels ν, at asymptotic distance (R ≈ 9a_o) with translational energy E_k. The vibrational dipole matrix elements were evaluated as functions of time t, whereby the maximum of |d⃑(t)|² is usually reached not too far from the turning point of the wave packet. Summation over 0 ≤ ν_CO ≤ 4 has already been carried out for the curves in Fig. 8(B). The dashed lines show some individual ν_CO contributions for the initial (vibrational) ground state. As the CO(X) bond length is significantly shorter than r at the minimum of well c̃, it is not surprising, that the final state turns out to be vibrationally hot. The rates per unit energy feature a rather broad structure with a maximum for E_k between 0.7 and 1.2 eV. Tendentially this peak is shifted towards higher E_k values with increasing initial (vibrational) excitation. The vibrational ground state of the c̃ potential well lies about 4.2 eV above the asymptotic limit H(1s) + CO(X, ν_CO = 0), whereas the vertical electronic energy difference down to the ground state PES is 3.1 eV. Thus radiation is predicted in the visible and near UV, and the (radiative) lifetimes resulting after integration over E_k are in the range 8 to 36 s.

7. Summary

The aim of this study was to reconsider the process of direct dissociative recombination of HCO⁺. To this end we have calculated new electronic PES for 2²A′ (B̃ state), 2²A′′ (C̃) as well as surfaces for 1⁴A′ (b̃), 1⁴A′′ (c̃) of HCO and for the ion HCO⁺, using a modified version of the MC-CEPA program ¹⁸ and preparing the required reference functions with iterated CI. The investigated DR mechanism involves the lowest core-excited ²Π state. The corresponding ²A′ and ²A′′ resonant state PES are merged diabatically into the B̃ and C̃ surfaces close by their respective (bent) minima upon bending the molecule. No additional diabatization was applied to connect the region of electron capture with the PES leading to dissociation. This situation is illustrated in Fig. 2 and 3. All examined core-excited doublet and quartet PES were found to intersect the ground state surface of HCO⁺ within the range of the vibrational zero point energy.

The electronic coupling of the initial state (free electron with energy ε plus ion) to the diabatic core-excited surfaces could be treated numerically in a restricted area with aid of an auxiliary wall potential (AWP)³⁶ allowing adjustment of ε (between 0.02 eV and 1.0 eV). We found that the coupling matrix elements increase (decrease) in value upon bending the molecule, for A′ (A′′) symmetry. At the equilibrium geometry of the ion our calculated autoionization width Γ(ε → 0) = 564 cm⁻¹ is almost twice the value obtained by Larson et al. (LStO) in 2012.¹⁷

Wave packet propagation on the diabatic surfaces has been carried out using extended wave function basis sets and transforming the Hamiltonian of nuclear motion to tridiagonal form. We start wavepackets at the HCO⁺ geometry only, not from the (by ca. 2 eV higher) HOC⁺ minimum. We estimate that at ε = 0.15 eV dissociation to H + CO(a³Π) is by a factor 2.5 to 5 more probable than to C + OH(X²Π). The dominance of the channel H + CO decreases with increasing ε. DR cross sections are shown in Fig. 5 and 7.

Transitions from doublet to quartet surfaces mediated by spin–orbit coupling, followed by dissociation, are a possible effect. We obtain coupling matrix elements in the order of about 50 cm⁻¹, and supported by a propagation method using coupled surfaces of the two multiplicities we predict that dissociation from 2²A′′ via 1⁴A′′ to H + CO(a³Π) is the most important pathway on a quartet surface, but its (relative) contribution to the DR cross section is only of the order of magnitude 10⁻⁴ compared to pure doublet dissociation.

Compared to the results of (LStO) our theoretical direct DR cross sections are significantly larger for ε < 0.4 eV, mainly due to the included angular dependence of autoionization width Γ (eqn (3)) and to the more favourable location of the intersection seam of the ²A′ PES with the ionic surface. At the lower end of our considered energy range (0.02 eV ≤ ε ≤ 0.2 eV) the larger contribution of the direct mechanism improves the agreement of theory and experiment. Thus the often stated dominance of the indirect process for HCO⁺ DR at ε ≤ 0.4 eV must be relativized. Our results support the supposition of Schneider et al.⁸ that the drop of σ(ε) predicted in ref. 9 is concealed by the direct DR mechanism. In other words, from the point where ε reaches the highest vibrational mode of HCO⁺ (3089 cm⁻¹) the direct process becomes dominant.

In contrast to the results of (LStO) we find the (lowest accessible) H + CO(a³Π) channel to contribute most to the cross section σ for ε < 0.5 eV. This seems reasonable from a classical mechanics point of view: for low ε the wave packet tends to move in the direction of the largest excess energy. For this DR channel our theoretical final state analysis leads to a similar (non-thermal) population of vibrationally excited states as the experiments of Rosati,²⁴ which however show a distinctly higher population of the ν = 0 level.

The low vibrational levels in the c̃ and c̃′ wells (see Table 3, c̃′ is a second minimum of the 1⁴A′′ PES on the HOC side) are of special interest. They may be populated by radiative transition from the small wave packet fraction transferred to the 1⁴A′′ state by spin–orbit interaction. Our calculations showed the population probabilities to be very small (compared to dissociation). But as the lowest quartet dissociation limits are higher than these vibrational levels, depopulation is expected to be very slow. Computation of radiative transition rates from the lowest vibrational c̃ levels to 1²A′ indicate that the 1⁴A′′ PES wells make up the best candidates for long living excited states of gaseous HCO, which might be relevant under certain astrophysical conditions.

Improved dynamics calculations including simultaneous propagation towards all three atom-diatom limits (and possibly coupling to the ²Φ and Rydberg states) might be the topic of a forthcoming paper,¹⁹ and we shall try to extend the range of collision energy ε. Time dependent calculations modelling the indirect DR process and its contribution to the branching ratios of the fragments are also desirable.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank the Ernst Miescher Foundation (Basel, Switzerland) for technical support. Moreover we wish to acknowledge the program development by M. Luo (spin–orbit coupling) and R. Guérout (auxiliary wall potential).

Notes and references

1. T. W. Hartquist, J. M. Pittard and S. A. E. G. Falle, Diffuse Matter from Star Forming Regions to Active Galaxies, 2007, ISBN 978-1-4020-5424-2.
2. D. R. Bates, Phys. Rev., 1950, 78, 492–493.
3. J. N. Bardsley, J. Phys. B, 1968, 1, 365–380.
4. I. A. Mikhailov, V. Kokoouline, Å. Larson, S. Tonzani and C. H. Greene, Phys. Rev. A, 2006, 74, 032707(9).
5. N. Douguet, V. Kokoouline and C. H. Greene, Phys. Rev. A, 2008, 77, 064703(4).
6. C. Jungen and S. T. Pratt, J. Chem. Phys., 2008, 129, 164311(9).
7. N. Douguet, V. Kokoouline and C. H. Greene, Phys. Rev. A, 2009, 80, 062712(6).
8. I. F. Schneider, N. Pop and C. Jungen, Phys. Rev. A, 2012, 86, 062706(11).
9. S. Fonseca dos Santos, N. Douguet, V. Kokoouline and A. E. Orel, J. Chem. Phys., 2014, 140, 164308(12).
10. N. Douguet, S. Fonseca dos Santos, V. Kokoouline and A. E. Orel, EPJ Web Conf., 2015, 84, 07003(10).
11. S. T. Pratt and C. Jungen, J. Chem. Phys., 2012, 137, 174306(6).
12. A. Le Padellec, C. Sheenan, D. Talbi and J. B. A. Mitchell, J. Phys. B, 1997, 30, 319–327.
13. M. Hamberg, et al. , J. Phys. Chem., 2014, 118, 6034–6049.
14. Å. Larson, S. Tonzani, R. Santra and C. H. Greene, J. Phys.: Conf. Ser., 2005, 4, 148–154.
15. Å. Larson and A. E. Orel, Phys. Rev. A, 2009, 80, 062504(9).
16. C. Nordhorn, et al. , J. Phys.: Conf. Ser., 2011, 300, 012004(8).
17. Å. Larson, M. Stenrup and A. E. Orel, Phys. Rev. A, 2012, 85, 042702(8).
18. R. Fink and V. Staemmler, Theor. Chim. Acta, 1993, 87, 129.
19. M. Jungen and M. Lehner, to be published.
20. Q. Chen, S. Han, X. Hu and D. Xie, Chin. J. Chem. Phys., 2022, 35, 303–310.
21. S. Han, et al. , Sci. Adv., 2019, 5, eaau0582.
22. G. Sun, et al. , Faraday Discuss., 2022, 238, 236–248.
23. G. Winterhoff, S. C. Galleguillos Kempf, P. Jensen and P. R. Bunker, J. Mol. Spec., 2018, 354, 71–82.
24. R. E. Rosati, M. P. Skrzypkowski, R. Johnsen and M. F. Golde, J. Chem. Phys., 2007, 126, 154302(8).
25. V. Staemmler, Theor. Chim. Acta, 1983, 64, 205–215.
26. J. M. Dyke, J. Chem. Soc., Faraday Trans., 1987, 83, 69.
27. H. J. Werner, C. Bauer, P. Rosmus, H. M. Keller, M. Stumpf and R. Schinke, J. Chem. Phys., 1995, 102, 3593–3611.
28. M. C. Chuang, M. F. Foltz and C. B. Moore, J. Chem. Phys., 1987, 87, 3855.
29. K. P. Huber and G. Herzberg, Constants of Diatomic Molecules, Van Nostrand, New York, 1979.
30. M. E. Jacox, J. Phys. Chem. Ref. Data, 2003, 32, 1–441.
31. T. A. Cool and X. M. Song, J. Chem. Phys., 1992, 96, 8675–8683.
32. P. J. Bruna, R. J. Buenker and S. D. Peyerimhoff, J. Mol. Struct., 1976, 32, 217–233.
33. K. Tanaka and E. R. Davidson, J. Chem. Phys., 1979, 70, 2904–2913.
34. P. O. Widmark, P. A. Malmqvist and B. O. Roos, Theoret. Chim. Acta, 1990, 77, 1007.
35. K. Kaufmann, W. Baumeister and M. Jungen, J. Phys. B, 1989, 22, 2223–2240.
36. R. Guérout, M. Jungen and C. Jungen, J. Phys. B, 2004, 37, 3043–3055.
37. W. Pauli, Z. Phys., 1927, 43, 601.
38. C. M. Marian, Rev. Comp. Chem., 2001, 17, 99–204.
39. A. I. Florescu-Mitchell and J. B. A. Mitchell, Phys. Rep., 2006, 430, 277–374.
40. A. E. Orel, Phys. Rev. A, 2000, 62, 020701.
41. M. J. Seaton, Comp. Phys. Commun., 1982, 25, 87.
42. C. W. McCurdy and T. N. Rescigno, Phys. Rev. A, 1989, 39, 4487.
43. A. E. Orel, T. Rescigno and B. Lengsfield, Phys. Rev. A, 1991, 44, 4328.
44. U. Galster, Phys. Rev. A, 2010, 81, 032517.
45. M. Lehner and M. Jungen, J. Phys. B, 2015, 48, 035101(9).
46. C. Lanczos, J. Res. Nat. Bur. Standards B, 1950, 45, 225.
47. M. R. Manaa and D. R. Yarkony, J. Chem. Phys., 1994, 100, 473–480.
48. C. A. Mead and D. G. Truhlar, J. Chem. Phys., 1982, 77, 6090–6098.
49. K. C. Kulander and E. J. Heller, J. Chem. Phys., 1978, 69, 2439.
50. C. W. McCurdy and J. L. Turner, J. Chem. Phys., 1983, 78, 6773.
51. N. G. Adams and L. M. Babcock, J. Phys. Chem., 1994, 98, 4564.
52. C. Zener, Proc. R. Soc. London, 1932, A137, 696.
53. M. Mladenović and S. Schmatz, J. Chem. Phys., 1998, 109, 4456–4470.

---