Collisional excitation of HCN by CO to refine the modeling of cometary comae

TL;DR

This paper provides the first dataset of collisional excitation and de-excitation rate coefficients for HCN induced by CO in cometary atmospheres, covering 5-50 K, to improve modeling accuracy of molecular distributions.

Contribution

It introduces a novel dataset of rate coefficients for HCN-CO collisions, derived using the SACM method validated against quantum calculations, enhancing cometary atmosphere models.

Findings

First dataset of HCN-CO collisional rates for 5-50 K.

Validated SACM approach against quantum calculations.

Improves modeling of HCN distribution in cometary atmospheres.

Abstract

We present the first dataset of collisional (de)-excitation rate coefficients of HCN induced by CO, one of the main perturbing gases in cometary atmospheres. The dataset spans the temperature range of 5-50 K. It includes both state-to-state rate coefficients involving the lowest ten and nine rotational levels of HCN and CO, respectively, and the so-called "thermalized" rate coefficients over the rotational population of CO at each kinetic temperature. The derivation of these coefficients exploited the good performance of the statistical adiabatic channel model (SACM) on top of an accurate interaction potential computed at the CCSD(T)-F12b/CBS level of theory. The reliability of the SACM approach was validated by comparison with full quantum calculations restricted at the lowest total angular momentum of the system. These results provide essential input to accurately model the distribution among the rotational energy levels and the abundance of HCN in cometary atmospheres, accounting for deviations from local thermodynamic equilibrium that typically occurs in such environments.

Figures

Tables

Table 1. Table 1: Geometric parameters and potential energy for stable structures in the PES. Energies are given relative to the asymptote. Units are Å , degrees, and cm -1 .

	GM 1	LM 1
$R$	$4.826$	$4.532$
${\theta }_1$	$180.0$	$0.0$
${\theta }_2$	$0.0$	$0.0$
$\phi$	—	—
$V$	$-604.13$	$-363.32$

Equations

$$k_{j_1{j}_2\to j_1^{\prime }{j}_2^{\prime }}(T)={\left(\frac{8}{\pi \mu k_B^3{T}^3}\right)}^{1/2}$$

$$\times {\int }_0^{\infty }{\sigma }_{j_1{j}_2\to j_1^{\prime }{j}_2^{\prime }}\left(E_c\right)E_c\exp \left(-E_c/k_{B}T\right)\mathrm{d}{E}_c,$$

$${\bar{k}}_{j_2\to j_2^{\prime }}(T)=\sum _{j_1}{n}_{j_1}(T)\sum _{j_1^{\prime }}{k}_{j_1{j}_2\to j_1^{\prime }{j}_2^{\prime }}(T),$$

$$n_{j_1}(T)=\frac{\left(2j_1+1\right)\exp \left(-\frac{E_{j_1}}{k_{\mathrm{B}}T}\right)}{{\sum }_{j_1^{\prime }}\left(2j_1^{\prime }+1\right)\exp \left(-\frac{E_{j_1^{\prime }}}{k_{\mathrm{B}}T}\right)},$$

Full text

1 Introduction

Cometary ices are among the most pristine materials of the Solar System, as they preserve the chemical conditions that existed in its earliest stages 1. Thanks to recent advances in observational techniques, we can now observe with very high spatial and spectral resolutions the molecular content of cometary comae, $i.e.$, the ensemble of gases which sublimate from cometary ices by solar radiation, forming an atmosphere around the cometary nucleus 2, 3, 4, 5. However, the accurate interpretation of these spectra requires sophisticated radiative transfer models. A key challenge in this regard is describing how a molecule’s energy is distributed among its rotational levels. At the typical low densities of cometary comae ($1<n<10^{-10}$ cm*-3*), this distribution can strongly deviate from the local thermodynamic equilibrium (LTE). In these cases, accurate datasets of collisional coefficients with the most abundant perturbing species are crucial (as an example, the reader is referred to Loreau et al. 6), and usually rely on state-of-the-art theoretical calculations.

However, when dealing with heavy projectiles such as \chH2O, \chCO and \chCO2 (the most abundant gases in cometary comae), collisional calculations notably complicate: the high densities of states of the projectiles dramatically increase the number of channels to be considered, which makes full quantum scattering calculations unaffordable in terms of computational time and memory requirements. To address this, the statistical adiabatic channel model (SACM) was recently tested as an alternative. Developed by Quack and Troe 7 and subsequently refined by Loreau et al. 8, this approach can generate reliable sets of collisional coefficients in favorable cases at significantly reduced computational cost. In terms of accuracy, the SACM is particularly effective for systems that form a stable intermediate complex, often showing discrepancies of less than a factor 2 with respect to the exact close coupling (CC) calculations 8, 9, 10, 11. For cometary applications, the SACM currently offers the most practical solution to overcome the prohibitive costs associated with heavy projectiles and to deal with the low temperatures ($\sim$ 5-50 K) that rule out the application of (quasi)classical methods ***The mixed quantum/classical theory (MQCT) offers a promising compromise in this direction, but it currently remains insufficient to accurately characterize collisional processes at temperatures below 50 K 12, 13.. For the \chCS–\chCO collisional system, it was recently shown that it can produce state-to-state outcomes differing by less than an order of magnitude from full quantum calculations 14.

In order to refine the modeling of cometary comae, we address in this work the collisional excitation of \chHCN by \ceCO. \chHCN is one of the most abundant nitrogen-bearing molecules in cometary comae 15. Detected in both millimeter and IR wavelengths, \chHCN primarily originates from direct sublimation of ices on the comet nucleus. Hence, it serves as a valuable tracer of the thermal and chemical history of pre-cometary material 16, 17. Understanding the abundance of \chHCN in comets compared to that in the interstellar medium provides insights into the early conditions of the Solar System and the possible delivery of prebiotic molecules to Earth 18, 19. Moreover, the abundance of \chHCN relative to other nitrogen-bearing species ($e.g.$, \chNH3, \chN2) can indicate the temperature at which cometary ices formed and the extent and nature of the chemical process they underwent 20, 21, 22. However, the only collisional investigations of cometary interest regarding \chHCN concern its excitation processes induced by \chH2O 23, 24. These data are not sufficient to model cometary comae at intermediate and large heliocentric distances, where the solar radiation is not intense enough to sublimate the water from the nucleus. Distant cometary comae are instead generally dominated by \chCO 25, 26, 27. These comets are particularly relevant for study of the composition of the early Solar System because they typically contain materials not previously altered by solar radiation. Moreover, the low densities that characterize cometary comae lead to strong deviations from LTE, making accurate collisional rate coefficients essential for reliable radiative transfer modeling.

This work aims to bridge this gap by providing the first dataset of collisional coefficients for the rotational (de)-excitation of \chHCN by \chCO. The methodology and results are presented as follows: the first Section describes the employed computational strategy. There, in the context of the Born-Oppenheimer approximation, the characterization of the interaction potential and the scattering calculations used to solve the nuclear Schrödinger equation are presented in two different subsections. The subsequent Section presents and discusses the results. Finally, in the last Section the main conclusions are drawn.

2 Computational Details

2.1 Potential Energy Surface

The Jacobi coordinates used to represent the four-dimensional (4D) HCN–CO interactions are depicted in Figure 1. Note that for the representation of the PES, CO is monomer 1, and HCN is monomer 2. $\vec{R}$ is the vector between the centers of mass of the two fragments, and $\vec{r}_{1}$ and $\vec{r}_{2}$ are vectors aligned with each molecule. Coordinate $R$ is the length of vector $\vec{R}$, while coordinates $\theta_{1}$ and $\theta_{2}$ represent (respectively) the angles between $\vec{R}$ and the vectors $\vec{r}_{1}$ and $\vec{r}_{2}$. The fourth coordinate is the dihedral (out of plane) torsional angle, labeled $\varphi$, which is the angle between the vectors $\vec{R}\times\vec{r}_{1}$ and $\vec{R}\times\vec{r}_{2}$.

In line with previous works on other van der Waals (vdW) dimers composed of two linear molecules,28, 29, 30, 31, 32, 33, 34, 35, 36 the analytical representation of the PES was constructed using an automated interpolating moving least squares (IMLS) methodology, freely available as a software package under the name autosurf.37 The fitting basis and other aspects of the IMLS procedure are similar to those employed for other previously treated systems that have been described in detail elsewhere.38, 39, 37, 40 The shortest intermonomer center-of-mass distance considered is $R=2.4$ Å. The short-range part of the PES was restricted by excluding regions with repulsive energies above a maximum of $6$ kcal/mol ($\sim 2\,000$ cm*-1*) relative to the separated monomers asymptote. The ab initio data coverage in the fitted PES extended to $R=25$ Å, while the zero of energy was set at infinite center-of-mass separation between the monomers. To guide the placement of high-level ab initio data, a lower-level guide surface was constructed using $3\,000$ points at the explicitly correlated CCSD(T)-F12a/VDZ-F12 level,41 distributed using a Sobol sequence 42 biased to sample the short range region more densely. Since the two fragments are not identical, and each has only C*∞v* symmetry, there is no additional symmetry to respect/exploit in the PES beyond that of the angle $\varphi$, for which energies were only computed in the reduced angular range: $0<\varphi<180^{\circ}$, as is the case in this treatment for all 4D systems.

The final high-level PES was computed using the CCSD(T)-F12b method, extrapolated to the complete basis set (CBS) limit. The CBS extrapolation was performed using the VTZ-F12 and VQZ-F12 bases 43 and the $l^{-3}$ formula. 44 All ab initio calculations were performed using the Molpro electronic structure code package.45

For construction of the PES, both monomers were held rigid. The bond distance for CO was fixed at $r_{\text{C}O}=1.128206$ Å, the vibrationally-averaged bond distance for the ground ro-vibrational state of CO, consistent with its rotational constant of $B(\text{CO})=1.9225$ cm*-1*. 46 For HCN, the equilibrium geometry is linear and the structural parameters employed in this study ($r_{\text{H}C}=1.0655$ Å and $r_{\text{C}N}=2.2187$ Å) are consistent with the experimental rotational constant, $B(\text{HCN})=1.4782$ cm*-1*. 47 Masses of $15.9949146221$, $12$, $1.007825032$, and $14.0030740052$ $u$ were used for 16O, 12C, 1H and 14N, respectively.

To represent the long range (out to arbitrary distances), the PES switches smoothly to an analytic expression representing electrostatic, induction, and dispersion interactions between the fragments (truncated at 8th order). For the high-level PES, $3\,841$ points were used, resulting in an estimated global root-mean-squared fitting error of $0.8$ cm*-1* (excluding the long range where the fitting error is extremely small). The estimated error is $0.1$ cm*-1* for geometries corresponding to energies in the wells (below the asymptotic zero). The analytical representation of the PES is available from the authors upon request.

Figure 3 shows radial cuts through the PES for nine different angular orientations. The points represent ab initio data (not included in the fit), while the lines plot the fitted PES. This helps one appreciate the anisotropy of the PES as reflected in the different values of $R$ at which the onset of a steep repulsive wall begins, as well as the accuracy of the fit.

The PES is characterized by two minima, the global minimum (GM), a collinear arrangement where the C-atom of CO approaches the H-atom of HCN, and a local minimum (LM), also collinear, where the O-atom of CO approaches the H-atom of HCN. These minima are seen across the middle of Figure 2, with GM appearing twice in the plot due to the use of extended angles. The energies and structural parameters are given in Table 1 where it is seen that GM ($E=-604.1$ cm*-1*) is much more stable than LM ($E=-363.3$ cm*-1*), although it is LM with a significantly shorter separation between monomers, due to the location of the CO monomer’s center-of-mass.

2.2 Scattering Calculations

The constructed state-of-the-art interaction potential was subsequently implemented in the scattering equations to solve the nuclear part of the Schrödinger equation and characterize the dynamics of the collisional system.

The most accurate method to solve scattering equations and to obtain the so-called scattering matrix, which contains all the collisional observables of the system, is the full quantum CC approach. However, for systems characterized by high densities of rotational states, combined with (a) deep well(s) of the intermediate complex (such as \chHCN interacting with \chCO), the number of coupled equations required for convergence becomes computationally prohibitive.

To address this, as introduced above, we decided to adopt the SACM approach. This approximation requires to compute only the adiabatic curves that are obtained by excluding the nuclear kinetic term and diagonalizing the interaction Hamiltonian in a basis of rotational functions for the two colliders. This noticeably simplifies the scattering calculations, as each adiabatic curve is independent of the energy and hence needs to be computed only once for each value of $J_{\text{tot}}$. The calculations were performed by using the HIBRIDON code 48, which contains an implementation that permits one to retrieve only the adiabatic states. The scattering matrices were then rebuilt by including all the open adiabatic channels at each total energy, assuming that they occur with the same probability. This yields to state-to-state collisional cross sections, which can be subsequently averaged over a Maxwell–Boltzmann distribution of energies to obtain the corresponding rate coefficients.

As mentioned in the Introduction, the SACM method assumes that the intermediate complex formed during each collision lives long enough for energies to be statistically redistributed. This is a good approximation for the \chHCN and \chCO system, where the hydrogen bonding interaction with \chHCN ($2.97$ D) stabilizes the complex significantly. In fact, the well depth of the \chHCN–\chCO complex is nearly three times greater than that of \chCS–\chCO, for which the SACM has already shown good accuracy.14 However, \chHCN has a lower density of rotational states compared to \chCS, and is therefore less efficiently stabilized in the potential well during collisions. On the other hand, this characteristic reduces the computational demand and enables the exploration of a broader energy range.

The cross sections were calculated for transitions between rotational levels up to $j_{2}=9$ for \chHCN and $j_{1}=8$ for \chCO, which are the most populated states up to 50 K. At this temperature, the adiabatic channels associated with energy levels up to $\sim 500$ cm*-1* would typically be included in the calculations. In this energy range, the accessible rotational levels of \chHCN and \chCO are $j^{\text{max}}_{2}=18$ and $j^{\text{max}}_{1}=15$, respectively. However, the inclusion of the full rotational basis led to a substantial memory demand, reaching over 15000 open channels even at only $J_{\text{tot}}=2$. In order to ease the calculations from a computational point of view, the rotational basis was further truncated to $j_{2}^{\text{max}}=12$ and $j_{1}^{\text{max}}=10$. The accuracy of this truncation was ensured by confirming that the convergence of the rate coefficients at $J_{\text{tot}}=2$ is within 20% with respect to the full rotational basis over the 5–50 K temperature range. Since the main contributions from higher values of $J_{\text{tot}}$ are expected to increase the centrifugal barrier for each collisional process, the resulting deviations can reasonably be considered within the uncertainty range of the fully converged calculations.

Using the reduced rotational basis, a partial wave expansion from $J_{\text{tot}}=0$ up to $J_{\text{tot}}=142$ was performed, ensuring convergence of the cross sections better than 5%.

Starting from the cross-sections for each $j_{1}j_{2}\rightarrow j_{1}^{\prime}j_{2}^{\prime}$ transition ($\sigma_{j_{1}j_{2}\rightarrow j_{1}^{\prime}j_{2}^{\prime}}$), the (de)-excitation rate coefficients ($k_{j_{1}j_{2}\rightarrow j_{1}^{\prime}j_{2}^{\prime}}$) from 5 to 50 K were derived by thermal averaging over the collisional energy ($E_{c}$):

[TABLE]

where $k_{B}$ is the Boltzmann constant.

In addition to the dataset of state-to-state collisional coefficients, thermalized rate coefficients for the (de)-excitation of \chHCN were also retrieved. These coefficients were obtained by summing the $k_{j_{1}j_{2}\rightarrow j_{1}^{\prime}j_{2}^{\prime}}(T)$ coefficients over all possible final states of \chCO and averaging over the initial rotational states of \chCO:

[TABLE]

where $n_{j_{1}}(T)$ is the relative population of the projectile:

[TABLE]

with $\left(2j_{1}+1\right)$ being the degeneracy of the rotational levels of \chCO.

Thermalized rate coefficients are commonly used in radiative transfer models 49 to address collisional excitation of the target molecule, under the assumption that the rotational temperature of the projectile is equal to the kinetic temperature of the surrounding gas. Therefore, this dataset can be applied to model the abundance of \chHCN in cometary comae thermalized with respect to \chCO.

To assess the accuracy of the SACM method, we compared a subset of rate coefficients with the results obtained with full quantum CC calculations. Due to the deep potential well of the system, a direct comparison of the inelastic cross sections could indeed be biased by the presence of resonances extending over a wide range of kinetic energies. To minimize this effect, instead of the inelastic cross sections, we directly compared the state-to-state rate coefficients, computed in an energy interval between 3 and 200 cm*-1*. The 3–20 cm*-1* range was sampled with very fine steps (0.2 cm*-1*), which gradually increased up to 2 cm*-1* above 100 cm*-1*. To optimize the rotational basis for both \chHCN and \chCO, we evaluated how the inclusion of different rotational levels affected the convergence of the inelastic cross sections and we selected the basis sets that ensured convergence within 2%. Specifically, from $2$ and $50$ cm*-1*, we included the first 22 rotational states of \chHCN and the first 19 states of \chCO. These were progressively increased up to 28 and 26, respectively, at $200$ cm*-1*. The angular part of the nuclear Schrödinger equation was restricted to a single partial wave $J_{\text{tot}}=0$ to keep the computational cost affordable in the cosidered energy range. The radial part was subsequently included by numerical propagation over the $R=2.9-42.0$ Å  range. We computed CC rate coefficients in the $5-50$ K temperature range and accounting for (de)-excitations among the rotational levels of \chHCN-\chCO below 80 cm*-1*. These were compared with their counterparts computed within the SACM method, using the same rotational basis used for data production but accounting only for the adiabatic channels up to $200$ cm*-1* and for $J_{\text{tot}}=0$. The comparison for $T=10,30$ and 50 K is reported in Fig. 4.

At these temperatures, most of the rate coefficients are reproduced within a factor of 2, and deviations do not exceed one order of magnitude, confirming the method’s reliability. This was somewhat expected, as the agreement is consistent with the findings of Godard Palluet et al. 14 for the \chCS and \chCO system.

The agreement between the CC and SACM approaches significantly improves when comparing thermalized rate coefficients, showing deviations of less than a factor of 2, as illustrated in Fig. 5. Thermal averaging over the rotational states of \chCO effectively smooths out the discrepancies between the two computational methods, making the SACM approach particularly accurate in capturing the overall collisional behavior of the system.

3 Results and Discussion

We derived the first dataset of (de)-excitation rate coefficients between the lowest ten rotational levels of \chHCN in collision with the first nine rotational levels of \chCO. While the full dataset will be made available through the EMAA†††https://emaa.osug.fr//, LAMDA50 and BASECOL51 databases, we summarize here the main outcomes and trends observed.

Fig. 6 shows the thermalized rate coefficients of de-excitation from the first ten rotational levels of \chHCN. For each initial $j_{2}$ level (indicated by a different color), the rate coefficients are plotted as a function of the $j_{2}^{\prime}$ final state. The weak dependence of these trends with respect to the kinetic temperature stands out: although absolute values of the rate coefficients are slightly higher at low temperature, the relative tendency remains largely consistent across the different thermal contributions. The same behavior (although not shown here) is observed for the state-to-state rate coefficients resolved for the rotational levels of the \chCO projectile. In general, the thermalized rate coefficients tend to decrease with increasing $\Delta j_{2}$, reflecting the growing energy gap that they address. However, exceptions to this trend emerge when the $j_{2}^{\prime}$ final state is equal to $1$ or $4$ . In these cases, regardless of the initial state, a slight increase of the rate coefficient value is observed, suggesting a favor for these particular transitions. At present, the authors have not identified a clear physical rationale for this behavior.

We now address the influence of the rotational state of the \chCO projectile, which is illustrated in Fig. 7. Here, the values of the state-to-state rate coefficients for de-excitations from the $j_{2}=9$ rotational state of \chHCN and for the transitions involving the lowest nine $j_{1}$ rotational states of \chCO are reported. The values are mapped with respect to $\Delta j_{2}$ (x-axis) and $-\Delta j_{1}$ (y-axis). Hence, the diagonal elements correspond to transitions where a de-excitation of \chHCN is matched by an excitation of \chCO of the same magnitude. Likewise, the elements of the first superdiagonal and subdiagonal represent transitions where the energetic jump between the levels of the two molecules differs by $\pm 1$, and so on. The largest values of the state-to-state rate coefficients are predominantly concentrated in this region, particularly for transition where $\Delta j_{2}=-\Delta j_{1},\,\,\Delta j_{2}=-\Delta j_{1}\pm 1,$ and $\Delta j_{2}=-\Delta j_{1}\pm 2$ . This behavior is a signature of near resonant energy transfer regimes 52, 53, 54, where the system tends to conserve its total internal energy during the collision. This phenomenon is observed when the colliding partners, such as \chHCN and \chCO, have similar rotational constants, allowing their energy levels to align closely in resonance. Noteworthy, the near resonant energy transfer is efficiently captured by the statistical SACM method, highlighting its robustness in describing even purely quantum effects.

We expect these results to be relevant to model \chHCN in cometary comae, where non-LTE radiative transfer is key for determining molecular abundances. As demonstrated in the previous section, the SACM approach provides sufficiently accurate rate coefficients, proving that the deep potential well and dense rotational structure of the system allow for a statistically valid redistribution of its energy after the formation of the intermediate complex.

4 Conclusions

The present study aims to support the interpretation of current and future observations of \chHCN in cometary comae. In these environments, \chHCN is largely abundant and widely spread, making it an important tracer of the thermal and chemical history of precometary material. Furthermore, \chHCN is considered a fundamental building block in prebiotic chemistry, and its presence in comets offers valuable insights into the processes that may have contributed to the origin of life in the primordial Earth. To contribute to this understanding, we have computed the first dataset of inelastic rate coefficients involving the ten lowest rotational levels of \chHCN in collision with the lowest nine states of \chCO, over a temperature range of 5–50 K.

We characterized the interaction potential using coupled-cluster electronic structure theory and an accurate fitting method. Subsequently, we performed scattering calculations by means of the SACM statistical approach. In order to facilitate direct application in non-LTE radiative transfer models, the computed state-to-state rate coefficients were thermalized over the rotational states of \chCO, assuming equilibrium between the rotational population of \chCO and the kinetic temperature.

The accuracy of the SACM method was assessed by comparison with a subset of full quantum CC calculations at $J_{\text{tot}}=0$. Given the complexity of the system and the lack of existing data, the SACM showed an encouraging agreement, reproducing state-to-state rate coefficients within a factor of 10 and thermalized rate coefficients within a factor of 2.

With the contribution of this work, accurate collisional datasets for \chHCN interacting with both \chH2O23, 24 and \chCO are now available. Given the ubiquity of \chHCN in cometary comae, these datasets offer a valuable test case to investigate the sensitivity of non-LTE effects under diverse physical conditions, such as different heliocentric distances. Worth to note, also \chCO2 is quite abundant in many cometary comae, especially at intermediate heliocentric distances. However, its high density of rotational states (resulting from a very low rotational constant) makes its collisional characterization significantly more demanding from a computational point of view. Although it is not expected to be the dominant collisional partner, computing rate coefficients with \chCO2 could nonetheless provide additional insights. All of these data will play a central role in driving the non-LTE distribution of population of \chHCN in cometary comae, bringing new knowledge to probe the chemical and physical conditions of the early Solar System.

{acknowledgement}

This work has been supported by Region Bretagne. The authors acknowledge the “Programme National de Planétologie” (PNP) under the responsibility of INSU, CNRS (France). FL acknowledges the Institut Universitaire de France. RD, EQS, and ABP are supported by the U.S. Department of Energy under Award DE-SC0025420.

5 TOC Graphic