Otoniel Denis-Alpizar

Thierry Stoecklin

Philippe Halvick

Marie-Lise Dubernet

Sarantos Marinakis

Potential energy surface and rovibrational energy levels of the H 2 -CS van der Waals complex

Keywords: Quantum molecular collision studies for processes of astrophysical interest potential energy surface, inelastic quantum dynamic, bound states, bending-rotation coupling, CS-H 2, HCN-H 2, HCN-He, C 3 -He

TITRE: Calculs de dynamique inélastique pour des collisions moléculaires d'intérêt astrochimique.

R ÉSUM É

L'analyse des conditions physico-chimiques régnant dans le milieu interstellaire (ISM) nécessite de connaître les constantes de vitesse de collision inélastique qui ont lieu plus fréquemment dans l'ISM. Nous avons à cette fin calculées les surfaces d'énergie potentielles ainsi que les états liés des complexes CS-H 2 , HCN-H 2 , HCN-He et C 3 -He. Nous avons déterminé pour la collision CS-H 2 les sections efficaces et les taux d'excitation collisionnels pour les premiers niveaux rotationnels. Des observations récentes suggérent que l'excitation des modes de pliage des molécules triatomiques doit être prise en compte dans les modèles astrochimiques. Nous présentons donc deux nouvelles approches théoriques permettant d'effectuer un traitement Close Coupling des collisions inélastiques d'un atome avec une molécule triatomique. Le couplage entre les mouvements de rotation et de pliage de la molécule est traité soit exactement dans le cadre de l'approximation du rigid bender (RBCC) ou de façon approximée en moyennant le potentiel d'interaction atome-molécule sur le mode de pliage de la molécule (RBAA). La méthode RBCC est appliquée à l'étude des collisions HCN-He et C 3 -He pour lesquelles les sections efficaces de transition entre niveaux rotationnels appartenant à des modes de pliage différents sont obtenues. Les résultats sont comparés avec ceux fournis par l'approximation du rotateur rigide linéaire. Dans le cas de la collision HCN-He ils sont aussi comparés avec ceux obtenus en utilisant l'approche RBAA. Nous montrons que les sections efficaces de transitions entre des niveaux rotationels appartenant à des niveaux de bending différents doivent être calculées au niveau RBCC. MOTS-CL ÉS: surface d'énergie potentielle, dynamique quantique inélastique, calcul d'états liés, couplage pliage-rotation, CS-H 2 , HCN-H 2 , HCN-He, C 3 -He.

Résumé du travail de thèse

La thèse est présentée sous la forme de 7 chapitres dont un d'introduction et un de conclusion. Hormis ces deux chapitres chacun est constitué d'une introduction puis de la ou des publications correspondantes, soit un total de six publications dans des revues à comité de lecture. Le chapitre 2 résume les connaissances théoriques fondamentales nécessaires pour aborder le traitement quantique d'une collision inélastique entre deux molécules. Il s'agit tout d'abord des méthodes ab initio de calcul de structure électronique qui permettent d'obtenir l'énergie d'interaction de deux molécules pour un ensemble de points dans l'espace des coordonnées du système bimoléculaire, dans le cadre de l'approximation de Born-Oppenheimer.

Ces énergies sont ensuite utilisées pour obtenir un modèle analytique de la surface d'énergie potentielle. La deuxième partie du chapitre traite de la théorie de la diffusion tout d'abord pour les collisions élastiques. Ce formalisme est ensuite étendu au cas des collisions entre deux molécules linéaires ou plus précisément de deux rotateurs couplés, puis au cas d'une collision entre un atome et une toupie symétrique. Le concept de résonance est introduit au passage ainsi qu'une méthode de calcul des états liés rovibrationnels du système basée sur une utilisation particulière des équations de diffusion.

Le chapitre 3 traite de la collision entre H 2 et CS. Il s'agit d'obtenir les sections efficaces d'excitation rotationnelle de CS avec les deux espèces para et ortho de H 2 . A cette fin, la première surface d'énergie potentielle complète est calculée pour ce système et un premier test de la qualité de cette surface est effectué par le calcul des états liés du complexe. Ce travail est présenté dans une première publication. La comparaison de ces résultats avec les données expérimentales s'avère excellente, prouvant ainsi la qualité de la surface obtenue. Puis les calculs des sections efficaces de collision sont effectués et comparés aux résultats précédents qui concernent la collision He-CS, ce qui fait l'objet d'une deuxième publication. Les astronomes ont en effet l'habitude de déduire les sections efficaces de collision avec para-H 2 de celles obtenues pour l'hélium en faisant une simple multiplication des sections efficaces de collisions par la racine carrée du rapport des masses relatives. Cette approximation, qui ne peut fournir les taux avec ortho-H 2 , s'avère correcte pour ce système. Pour les collisions de para-H 2 avec CS, une préférence pour les transitions rotationnelles avec j impair est observée pour les faibles énergies de collision tandis que cette tendance s'inverse pour les énergies de collision plus importantes. La même observation est faite pour les collisions impliquant ortho-H 2 , l'inversion se produisant toutefois à plus haute énergie de collision pour cette forme de H 2 .

Le chapitre 4 présente une étude similaire, mais effectuée pour les collisions entre H 2 et HCN. La surface d'énergie potentielle est calculée ainsi que les états liés du système, ce qui constitue la première publication. La dynamique fait l'objet d'une publication séparée qui n'était pas encore soumise au moment de la soutenance de thèse. Les états liés obtenus sont en excellent accord avec les mesures expérimentales démontrant la qualité de la surface calculée. Le calcul de l'énergie de dissociation du complexe H 2 -CS donne 37.79 cm -1 pour l'espèce para et 60.26 cm -1 pour l'espèce ortho de H 2 . Le nombre d'états liés supporté par le puits de potentiel est de 101 pour le complexe de CS avec para-H 2 et de 330 pour celui avec ortho-H 2 . Un résultat particulièrement intéressant de ce travail est la comparaison de la surface moyennée sur l'état rotationnel fondamental de para-H 2 avec la surface originale. En effet une expérience pour ce système avait montré que le complexe He-HCN devait être linéaire en contradiction avec les seuls résultats théoriques disponibles à ce moment là. La surface moyennée donne un complexe linéaire et confirme donc les résultats expérimentaux.

Avec le chapitre 5, nous entrons dans la partie la plus novatrice de cette monographie. Dans le cas des collisions non-réactives, le calcul des sections efficaces de transition rotationnelle est quasiment toujours réalisé dans l'approximation des monomères rigides. En effet, puisque les forces générées par les interactions intermoléculaires sont de bien plus faible intensité que les forces intramoléculaires, on peut considérer que le couplage entre les mouvements intramoléculaires et intermoléculaires est négligeable. En diminuant le nombre de degré de liberté, cette approximation réduit fortement la complexité et le coût des calculs de dynamique quantique collisionnelle. Cette approximation permet d'obtenir de très bons résultats lorsque les fréquences vibrationnelles intramoléculaire sont bien plus grandes que les fréquences associées aux modes de vibration intermoléculaires. On peut citer le cas du système H 2 + CO pour lequel un excellent accord a été obtenu entre les données expérimentales et le calcul quantique du spectre de transition infrarouge, ainsi que pour des sections efficaces de collision inélastique à basse température.

Tant que l'on considère des collisions entre molécules diatomiques, on reste dans le domaine d'application de l'approximation des monomères rigides. Mais si on s'intéresse aux molécules triatomiques (ou plus grandes) alors il est plus difficile de justifier cette approximation puisque certains modes de mouvement peuvent être associés à des fréquences de vibration faibles. C'est le cas des mouvements de pliage pour certaines molécules triatomiques, et de torsion pour des molécules de taille supérieure. Dans ce chapitre, nous nous intéressons au cyanure d'hydrogène (HCN) en collision avec He. La fréquence de pliage de HCN est ν 2 = 712 cm -1 . Dans le chapitre suivant, nous traiterons le cas du tricarbone, une molécule bien plus souple dont la fréquence de pliage est très basse, ν 2 = 63 cm -1 . HCN, et son isomère HNC, figurent parmi les molécules organiques les plus abondantes dans le milieu interstellaire. Les émissions de leurs transitions rotationnelles constituent un marqueur important du gaz moléculaire dense dans les galaxies qui sont particulièrement brillantes dans l'infrarouge. Pour remonter à l'abondance de HCN à partir des observations astronomiques, lorsque la densité du gaz est trop faible pour que s'établisse un équilibre thermodynamique, il est nécessaire de connaitre les constantes de vitesse associées à l'excitation (ou désexcitation) rotationnelle par collision avec les espèces les plus abondantes, H 2 et He. Pour étudier la collision de HCN avec He, nous avons défini un modèle mathématique de l'interaction entre HCN et He, où HCN est représenté par un modèle semi-rigide: l'angle de pliage est une variable dynamique alors que les distances internucléaires sont des constantes.

Les paramètres de ce modèle ont été ajustés par la méthode des moindres carrés sur un ensemble de 43015 points dans l'espace des 4 coordonnées du système He-HCN. Les énergies ab initio ont été calculées par la méthode CCSD(T) avec une base aug-cc-pVQZ à laquelle on a ajouté un jeu d'orbitales atomiques diffuses sur le centre de l'axe reliant HCN à He, ce qui permet d'améliorer la description de la liaison faible entre les deux monomères. Nous avons obtenu ainsi une surface d'énergie potentielle dont le minimum global a une énergie de dissociation de 30.35 cm -1 et une configuration géométrique linéaire He-H-C-N avec une distance R=7. 94 Bohr entre He et le centre de masse de HCN. Le calcul des états liés de ce système a été réalisé par deux méthodes différentes: la première utilise l'approximation des monomères rigides (RMA), c'est-à-dire que l'on considère HCN comme une molécule linéaire rigide, et la deuxième méthode prend en compte le mouvement de pliage de HCN par le calcul d'un potentiel d'interaction moyen sous l'effet de la densité de probabilité de la fonction d'onde de pliage (RBAA). Les deux méthodes donnent des énergies de transition proches des valeurs expérimentales, la deuxième méthode étant légèrement plus précise. Pour le calcul des sections efficaces de relaxation rovibrationnelle, nous avons développé une troisième méthode entièrement nouvelle dans laquelle le mouvement de pliage de HCN et ses couplages avec la rotation et le mouvement de He sont totalement pris en compte par un formalisme de collision quantique (RBCC). Le traitement exact du couplage rovibrationnel pour une molécule linéaire telle que HCN permet de d'accéder au calcul du l-doubling, c'est-à-dire la levée de dégénérescence des états de pliage vibrationnel qui sont doublement dégénérés en l'absence de rotation. Dans le cas des transitions purement rotationnelles de HCN, les sections efficaces que nous avons calculées ne montrent pas de différences majeures entre les trois méthodes RMA, RBAA et RBCC. Ceci montre que l'approximation des monomères rigides reste parfaitement valable pour une molécule triatomique dont la fréquence vibrationnelle de pliage n'est pas trop basse. Si on considère maintenant les transitions rovibrationnelles de HCN, alors on observe des différences significatives entre les sections efficaces calculées avec la méthode RBAA et celles avec RBCC. L'ensemble de ces résultats a fait l'objet de deux publications.

Le chapitre 6 est dédié au système He-C 3 . Comme dans le chapitre précédent, nous avons défini une surface de potentiel de haute qualité et calculé les états liés avec les méthodes RMA et RBCC. La différence avec HCN est que la molécule C 3 est très souple, puisque sa fréquence de vibration de pliage est assez proche des fréquences vibrationnelles associées aux modes de mouvements intermoléculaires. Les calculs ab initio de 26526 points dans l'espace des coordonnées ont été réalisés avec la méthode CCSD(T) avec une base aug-cc-pVQZ à laquelle on a ajouté un jeu d'orbitales atomiques diffuses centré sur la liaison intermoléculaire. La géométrie d'équilibre du système He-C 3 ainsi calculée correspond à une structure en T légèrement dissymétrique, où le monomère C 3 est légèrement non-linéaire et l'énergie électronique de dissociation est 26.9 cm -1 . En raison de cette structure en T, la rotation de He autour de C 3 correspond à un potentiel d'interaction très anisotrope. La différence entre les énergies des états liés calculées par les méthodes RMA et RBCC est assez faible, ce qui montre que malgré la faible fréquence vibrationnelle de pliage de C 3 , le couplage entre ce mouvement vibrationnel intramoléculaire et les mouvements intermoléculaires reste faible. Ce travail sur le système He-C 3 a fait l'objet d'une publication.

Dans le dernier chapitre, nous présentons la conclusion générale de l'ensemble des travaux présentés dans cette monographie. Quatre systèmes de van der Waals ont été étudiés. Pour chacun, nous avons modélisé l'énergie d'interaction entre les deux monomères et calculé par des méthodes quantiques les énergies des états liés et les sections efficaces de collision inélastique. Les deux premiers systèmes, H 2 -CS et H 2 -HCN ont été étudié dans le cadre de l'approximation des monomères rigides. Les deux suivants, He-HCN et He-C 3 ont fait l'objet de travaux théoriques plus approfondis. Le mouvement intramoléculaire de pliage est pris en compte dans la dynamique, aussi bien pour les états liés que pour les sections efficaces de collision inélastique. Nous avons ainsi montré que l'approximation des monomères rigides donne de bons résultats pour les énergies des états liés et pour les sections efficaces tant que l'on s'intéresse seulement aux transitions rotationnelles. A partir d'une énergie de collision suffisante pour ouvrir le premier niveau vibrationnel excité, alors il est nécessaire d'utiliser notre nouvelle méthode RBCC pour calculer correctement les sections efficaces de transitions rovibrationnelles.

List of Figures

CHAPTER 1

General Introduction

The knowledge of the interstellar medium (ISM) has changed dramatically in the last century. From the idea of vacuum among the star, nowadays the development of the astrochemistry has given us a clear understanding of the main chemical and physical processes in different interstellar environments.

In 1926, Eddington proposed the hypothesis of the existence of matter in the interstellar space [1]. However, the first observation of molecular lines from the interstellar medium was done in 1934 [2], and in 1935 Rusell [3] suggested that these lines had a molecular origin. Up to the 1950s, the molecules detected in the ISM were limited to a few ones. In the 1960s, the technical progresses of radio-astronomy allowed the detection of first interstellar OH [4] and since then of many other molecules.

New recent observational facilities have opened the doors to the golden age of the studies on the interstellar medium. The instruments of the Hubble Space Telescope are able to examine the far-ultraviolet light and the near-ultraviolet light. In the same domain is the FUSE (Far Ultraviolet Spectroscopic Explorer) which also provides a large amount of data. The SOFIA (Stratospheric Observatory for Infrared Astronomy) and the Herschel Space Observatory provide complementary access to the infrared and submillimeter wavelengths domain. Finally, we want to mention the capabilities offered by ALMA (Atacama Large Millimeter/sub-millimeter Array) interferometer, which exhibits an impressive sensitivity and a high spatial and spectral resolution. ALMA will be fully operational in 2015 and will provide every day a huge amount of data of molecular lines observations. About eighty years after the first observation of a molecule line in the ISM, the Cologne Database for Molecular Spectroscopy (CDMS) [5-7] listed around 180 molecules detected in the ISM or in circumstellar shells (March 2013).

The densest regions in the ISM are the molecular clouds. In these regions, the formation of molecules is possible while the ionized gas is predominant in other regions of the ISM. The molecular clouds are furthermore the first stage in the formation of stars and for this reason they are of great interest for the astronomers. The detected emission and absorption spectra of molecules are the main source of the information on the ISM. Modelling the intensity and observed profiles allows the determination of the main physicochemical characteristics, density, temperature, and chemical abundance of the interstellar molecules. The two more important processes leading to transitions among internal levels of molecules in a typical interstellar cloud are the radiative and the collisional one with the most abundant neutral species [8]. The analysis of the radiative processes requires the knowledge of the Einstein coefficients, which are usually well known while for the collisional processes, the collisional rate coefficients of the detected interstellar molecules with the most common colliders in the ISM (e.x., H, He and H 2 ) are required and are usually unknown. If we consider collisions to be the dominant process, it is possible to make the hypothesis of the local thermal equilibrium (LTE). In this case, the rotational population is described by a Boltzman distribution at a given temperature. As the ISM has a very low density, the interactions are rare and the LTE is not a good approximation. When the LTE does not apply, the analysis of observed spectra requires the knowledge of the collisional rate between specific quantum levels. The routine analysis of the observation made in the ISM by the astronomers relies on tracers molecules which can be in small quantity, but are easy to detect. Among the most used molecules as dense mass tracers are the CO, HCN, HCO + and CS [9] molecules. For this reason, the study of the collision dynamics of these molecules is of paramount importance.

The very specific physical conditions typical of the interstellar clouds make the experimental approach rather complicated while the very low temperature and density facilitates the theoretical studies by reducing the number of open channels [8]. These two observations motivated the develop-ment of several quantum mechanical approaches to generate the requested theoretical data. In the 1960s Arthurs and Dalgarno [10] developed the formalism of quantum scattering for an atom colliding a rigid rotor within a time independent approach. Shortly after, the generalization to the collision between two rigid rotors was done [11, 12] and the corresponding computational packages started to be developed. Today, while a few groups like us develop their own package [13, 14] many studies are based on the use of the MOLSCAT [15] or Hybridon [16] packages. The amount of theoretical data available for the (de-)excitation of simple molecule by collisions with He, H 2 and H has consequently grown over the years and the necessity was made to develop data bases like BASECOL [17] where these informations can be found. A recent review of Roueff and Lique [18] summarize these studies (see references therein). An example of such studies will be presented in Chapter 3 for the H 2 -CS collision.

While it was thought that these studies could be limited to the rotational degrees of freedom, owing the very low temperature of the interstellar clouds, several recent observations suggest that the vibrational excitation of triatomic molecules in the ISM, ex. HCN [19, 20] and C 3 [21, 22], at least in the bending motion needs to be considered in the collision mechanisms. The investigation of these new topics is an important part of the present PhD work which is founded by the European Astronet network CATS. The vibrational (de-)excitation of triatomic molecules by collision with an atom or a diatomic molecule has motivated the development of several theoretical approaches [23-25]. However, all these pioneering studies neglect the rotational degree of freedom or are based on the use of drastic approximations such as the infinite order sudden (IOS) approximation [26-30].

The most recent of these approaches, dedicated to the collision of CH 2 with He, was published at the same time than our first paper dedicated to the He-HCN collision. The authors treated both the rotational and the bending motion [31, 32] at the ab initio level but averaged the interaction potential over the bending wave function of the CH 2 molecule before performing the dynamics of the collision between He and the rigid equilibrium structure of CH 2 . We developed a similar approach in our first paper dedicated to the He-HCN collisions as it will be presented below, the main difference between these two studies being the equilibrium geometry of the molecule which is linear for HCN leading to l-doubling while it is planar for CH 2 . Another important difference between these two systems is the magnitude of the bending frequency which is 750 cm -1 for HCN [33] while it is 963 cm -1 for CH 2 [34]. The coupling between bending and rotation which furthermore leads to l-doubling is then expected to be important for HCN and required the development of new theoretical approaches going beyond the rigid bender averaged approximation. This exact treatment of the bending-rotation interaction was made even more necessary by the recent detection of the HCN l-type transitions in the ISM [35] which cannot be accounted by the averaged potential over the bending wave function. This new method will be presented and applied to HCN in Chapter 5. Among the small carbon chains expected to play an important role as building block in the formation of complex organic molecules in the ISM, the linear C 3 molecule has been detected through several of its rovibrational transition. Its very small bending frequencies, around ten times smaller than HCN, makes furthermore of the He-C 3 collision the ideal system to test our method as will be shown in Chapter 6.

In summary, the present work gathers the studies of the collisions of some interstellar molecules (CS, HCN and C 3 ) with H 2 and/or He. For the two triatomic molecules, the bending-rotation interaction was included for their collisions with He. The thesis is divided into seven chapters.

In Chapter 2, we introduce the general background of the electronic calculations and we present the scattering formalism of two rigid rotors in the space fixed frame. We also give a brief account of the close coupling bound state method which we use to calculate the bound states energies of the Van der Waals complexes.

In Chapter 3, we study the collision between CS and H 2 . Both molecules are treated as rigid rotors. A new potential energy surface (PES) for the CS-H 2 calculated at coupled cluster level is first presented. We then report the bound levels for the CS-para-H 2 and CS-ortho-H 2 complexes and the results of the close coupling calculations which are compared to previous works.

In Chapter 4, we present a new four dimensional analytical model of the PES for the HCN-H 2 system when treating both HCN and H 2 as linear rigid rotors. A first use of this PES to compute the bound energy levels of HCN-para-H 2 and HCN-ortho-H 2 complexes is reported. The accuracy of the PES is discussed by comparing the transition frequencies among the calculated levels with the experimental data.

In Chapter 5 which is the central one of the present manuscript, a new 4D PES for the He-HCN system is first reported which takes fully into account the coupling between the rotation and bending of HCN. The details of the two different methods which are developed to treat the dynamics of the He-HCN collisions using this surface are then presented. In the first approach which we call the rigid bender averaged approximation (RBAA) the interaction potential is averaged over the bending wave function before performing the dynamics while in the second, called rigid bender Close coupling (RBCC), the coupling between rotation and bending is treated exactly within the rigid bender approximation. The calculations of the dynamics is then performed using both the RBAA and RBCC approaches and the results discussed. The computation of the bound states of the He-HCN complex are also reported first when treating the HCN molecule like a rigid linear molecule while in a second time the bound levels of the complex are calculated using the interaction potential averaged over the bending wavefunctions. Both calculation are compared with the available experimental and theoretical data.

In Chapter 6, we present a new PES for the C 3 -He complex which describes the coupling between bending and rotation and compute the bound states using the RBCC method. We compare these levels with those obtained treating C 3 like a linear rigid rotor. The results of the dynamics are presented and discussed for this system. Finally, the conclusions and perspectives of this work are presented in the Chapter 7. In this chapter we will present a short overview of the theoretical methods used to calculate inelastic cross sections of two simple molecules. The first traditional approximation which is used takes advantage of the difference of masses between the electron and the nuclei to assume that the electrons can respond almost instantaneously to the displacements of the nuclei. The use of this approximation allows dividing the theoretical work in two parts. First the electronic energy has to be calculated for a grid of fixed positions of the nuclei and the different methods available to perform this computation will be presented in Section 2.2. The resulting grid of energy points which is called a potential energy surface (PES) is then fitted using an analytical model. These models being system specific will be presented in each chapter. The second main part of such studies is dedicated to the dynamics of the nuclei and will be presented in Section 2.3.

CHAPTER 2

Born-Oppenheimer approximation

The time-independent Schrödinger equation for a molecular system formed by N atoms, without external field and relativist effects, is

HΨ(q, Q) = {T N (Q) + T e (q) + V (q, Q)}Ψ(q, Q) = EΨ(q, Q).
(2.1)

Q are the nuclear coordinates while q correspond to all electron coordinates. T N (Q) and T e (q) are the kinetic operator of the nucleus and the electrons, respectively. Ψ is the wavefunction depending of q and Q, and E is the total energy of the system. V (q, Q) is potential energy of the system and include the nuclei-nuclei, electron-nuclei and electron-electron interactions.

The wavefunction can be expressed as the product of a function that depend of the electronic coordinates (and parametrically of the nuclear coordinates) and a function which depend only of the nuclear coordinates:

Ψ(q, Q) = ∞ i=1 ψ i (q, Q)φ i (Q).
(2.2)

With ψ i (q, Q)|ψ j (q, Q) q = δ ij , where the sub-index q means that the integration is over the electronic coordinate only. The ψ i (q, Q) are the eigenfunctions of the electronic Hamiltonian,

[T e (q) + V (q, Q)] ψ i (q, Q) = E e ψ i (q, Q). (2.3)
If we substitute 2.2 in 2.1, use 2.3, multiply by ψ * i (q,Q) on the left and integrate over the electronic coordinates, we get

T N (Q) + E e φ i (Q) + W(q, Q) = Eφ i (Q).
(2.4)

The term W is

W(q, Q) = ∞ j 2 M ψ i |∇ Q |ψ j ∇ Q + 1 2 ψ i |∇ 2 Q |ψ j φ j (Q). (2.5)
The quantity W is non-zero because ψ i (q, Q) depend on the nuclear coordinates and thus ∂ψ i (q, Q) ∂Q is non-zero. The nuclear mass M in the denominator makes W small, and then it can be neglected.

The first step is to solve the electronic Schrödinger equation, which can be done fixing the nucleus in a given configuration. If the process is repeated for various molecular geometries, we obtain a set of electronic energies depending of the nuclear coordinates. It is possible to used this set of energies to get an energy surface V e (Q) which depend only of the nuclear coordinates.

Ab-initio Approaches

The ab initio approaches solve the electronic Schrödinger equation without using any empirical parameter. Indeed, the term ab initio comes from the Latin words from the beginning. These methods choose a model for the electronic wavefunction and solve the electronic Schrödinger equation using only the values of the fundamental constants and the atomic numbers of the nucleus.

The ab initio methods have been discussed extensively in several books about electronic structure (ex. [36-39]), for that reason we will outline only some aspects here. These methods are accurate at the expense of a high computational cost. Several commercial package such as Molpro [40] and Gaussian [41] implement most of these methods.

Hartree-Fock Approximation

The starting point for most of the ab-initio methods is the Hartree-Fock (HF) or Self-Consistent Field (SCF) method. The electron-electron interaction is the main difficulty of the electronic calculations. The HF approximation considers each electron moving in the electrostatic field of the nuclei and the average field of the other N -1 electrons. The wave function is a Slater determinant constructed from a set of functions called spin-orbitals. A spin-orbital is the product of a spatial wave function and a spin wave function. This determinant satisfies the symmetry requirements imposed by the Pauli principle.

The solutions are generated by a variational procedure. A trial set of spin-orbitals is used to solve the HF equations. A system with N electrons can be expressed by m spin-orbitals (with m ≥ N ). The calculated spinorbitals are used to construct a new Fock operator, which in turn is used to determine new spin-orbitals. The procedure continues until the energy and the spin-orbitals coefficients converge. At the end of the calculation, the spin-orbitals are arranged in order of increasing energy. The first N spinorbitals are called the occupied spin-orbitals and the remaining are called virtual spin-orbitals.

Electron Correlation

The HF approximation gives a quite good description of the electronic wavefunction. Indeed, the HF energy is generally more than 99% of the total electronic energy. But this accuracy is not sufficient for the study of chemical processes. As a matter of fact, let's consider some molecule with a total energy of 100 E h (for example some diatomic molecule with C or N or O atoms). Then the missing 1% is 1 E h and this error is dependent of the nuclear coordinates. But the accuracy wanted for the study of chemical processes is about 1 cm -1 , i.e. 5 × 10 -6 E h . In other words, the HF method gives the total energy with 2 or 3 correct digits while 8 or 9 are needed. The difference between the HF energy and the exact electronic energy is called the correlation energy. This energy arises from the pairwise electron-electron repulsions which is not taken into account in the HF method because the quantum mechanical effect on electron distributions due to the other N -1 electrons is treated only in an averaged way. Many methods have been developed to calculate the contribution of electron correlation to the total energy. We briefly present the main features of some of them in the following paragraphs.

Configuration interaction

In order to calculate the correlation energy, the electronic wavefunction is expressed as a linear combination of the HF determinant and many other determinants built by substitution of one or more occupied spin-orbitals by virtual spin-orbitals. The energy is then calculated by the variational method. If only one substitution (or excitation) is made in all determinants, the method is called the configuration interaction with single-excitation (CIS). In the case of the inclusion of single and double excitations, then it is the CISD method. The full-CI includes all possible excitations and gives the best possible treatment for nonrelativistic molecular system within the limit of the basis of atomic orbitals used. However, it is very expensive in computation time and memory computing.

Møller-Plesset many-body perturbation theory

The correlation energy being significantly smaller than the total energy, it appears possible to apply the perturbation theory by considering the correlation as a perturbation. In the Møller-Plesset perturbation theory (MPPT), the zero-order Hamiltonian is the sum of one-electron Fock operators. The ground state wavefunction is an eigenfunction of this Hamiltonian with an eigenvalue given by the sum of the energies of all occupied spin-orbitals. The perturbation is the difference between the averaged potential and the exact electronic repulsion. The HF energy is the sum of the zero-order energy and first order correction energy. The correction to the energy arises from the second-order perturbation theory (MP2). The use of third-order (MP3) and fourth-order (MP4) is also usual. These calculations are faster than the variational CI calculations, but of course do not obey to the variational principle. This means in practice that we cannot assume that the lowest energy is the closest energy to the exact value.

The coupled-cluster method

The CISD method allows to compute a significant part of the correlation energy, but has a major drawback, the size-consistency error. It has been shown that the CISD energy of N non-interacting identical molecules is not proportional to N . The attempts to suppress this error have resulted in a new method, namely the coupled-cluster (CC) method. This is the method used in this work, because it has been recognized more accurate than CISD and perturbation methods. The CC methods introduce the cluster C operator. This operator links the exact electronic wavefunction and the HF wavefunction as:

Ψ CC = e C Ψ 0 , (2.6)
where the exponential operator e C is defined by the Taylor series expansion

e C = 1 + C + 1 2! C 2 + 1 3! C 3 + • • • (2.7)
and where C is the sum of the one-electron excitation operator C 1 , twoelectron excitation operator C 2 until the N -electron excitation operator:

C = C 1 + C 2 + • • • + C N . (2.8)
The operator C i act over the HF wave function Ψ 0 and generates a sum of Slater determinants where i occupied spin-orbitals have been replaced by virtual spin-orbitals. In the case of single and double excitations we can write

C 1 Ψ 0 = a,p t p a Ψ p a (2.9) C 2 Ψ 0 = a,b,p,q t pq ab Ψ pq ab , (2.10)
where the t p a and t pq ab are numerical coefficients called single-excitation amplitudes and double-excitation amplitudes, respectively. Substituting the expressions 2.6-2.8 in the Schrödinger equation and taking only a few terms of the operator C, gives a set of not-linear equations which are then solved for the excitation amplitudes. With these amplitudes, one can determinate the CC wave function and the energy.

The CC method has different expansion orders, when C is approximated only with C 2 is called CCD. If it is included C 1 and C 2 , the method is called CCSD. The CCSDT furthermore includes the term C 3 , and the CCSD(T) means that the triple excitations are included perturbatively rather than exactly. This last method computes accurate energies with a tractable computational time.

As this method uses the spin-orbitals computed within the single-reference HF wave function, it is dependent of the quality of these spin-orbitals. Consequently, it is only when the electronic wave function is well represented by the HF determinant that the CC method can reach its maximum accuracy. Two criterion have been proposed to quantify the quality of the HF wave function [42, 43]. The first practical indication of this quality, based on single-excitation amplitudes, is the T 1 diagnostic. The study of Lee and Taylor [42] for several systems concluded that a T 1 diagnostic less than 0.02 indicates that the CCSD(T) energy is as accurate as possible. However the study of molecular systems with a large number of electrons has shown that the T 1 diagnostic is less reliable than for the small systems. This happen when the major part of a molecule is well described by the HF wave function and that only a minor part is poorly described. To address this shortfall, a more sensible criterion was proposed, the D 1 diagnostic [43].

Finally, let us mention the explicitly correlated method. When applied to the CCSD(T) method, this yield the CCSD(T)-F12 approximation [44]. In general, the explicitly correlated methods are based on the fact that the Slater determinants fail to model the exact wave functions at short interelectronic distances [45]. To have a better representation of the correlation, a term that depends explicitly of the inter-electronic distances is included in the wave function.

Multiconfiguration and multireference methods

The electron correlation methods discussed so far are mono-reference methods. The starting point is the HF determinant. The HF orbitals remain fixed in the subsequent calculations. There are many chemical cases where the HF wave function is a poor representation of the electronic wave function. Let's mention, in a non-exhaustive list, the cases of excited electronic states, degenerated or quasi-degenerated states, biradical systems, bond breaking or bond formation. In order to represent properly the electronic wave functions for such cases, multiconfigurational methods are necessary.

In the multiconfiguration self-consistent field method (MCSCF) the wave function is a linear combination of Slater determinants. The optimization process is performed on both the orbitals and configurations coefficients simultaneously. The simultaneous optimization makes this method computationally expensive. One of the most efficient approach for solving the MCSCF equations is the complete active-space self-consistent field method (CASSCF). The orbitals are divided into three classes: inactive orbitals, active orbitals and virtual orbitals. The inactive orbitals are occupied by two electrons, a full-CI expansion is performed within the active space which is span by the active orbitals, and the virtual orbitals are unoccupied.

Multireference (MR) methods include the instantaneous interactions among the electrons. Multiconfigurational wave functions are the star-ting point of the MR methods. The multiconfigurational wave function can be computed with a MCSCF, ex. a CASSCF method. Multireference configuration interaction (MRCI) and multireference perturbation theory (CASPT2) are commonly used while multireference coupled cluster methods are still under development.

Basis sets

In all ab-initio methods, the selection of the basis sets has an important role. Indeed, we need to expand the molecular orbital in a finite basis set. That creates an error called basis-set truncation error and make necessary to check the convergence of the basis set used.

Another important point to take into account is the form of functions. The most common are the Slater-type orbitals (STOs) and the Gaussiantype orbitals (GTOs). Owing to its mathematical form, a comparable representation of the wave function requires a greater number of GTOs than the STOs. In spite of this, the GTOs are more commonly used. The relative ease of calculation of the integrals that involve these functions makes their use more computationally efficient.

The basis set are classified according the number of functions used to represent the atomic valence shells. In the minimal basis set, only one function is used to represent each valence orbital. In the double-zeta (DZ) basis set, two basis functions are used for each valence orbitals, and this yield a significant improvement. Further basis sets, which can be formed in the same way, are the called triple-zeta (TZ) basis set, quadruple-zeta (QZ) basis set and quintuple-zeta (5Z) basis set.

The bond formed in molecules distort the atomic orbital due to the adjacent atoms. Considering this displacement of the center of charge, it is necessary to include polarization functions. These functions can be orbitals with high values of l. Adding such polarization functions to a DZ basis creates the called double-zeta plus polarization (DZP) basis set.

In the correlation-consistent (cc) basis set, the functions are designed for a better description of the correlation energy. The basis set can be improved also with the addition of diffuse functions (aug-). These are functions with a small exponent that allow a better description of the wavefunction far from the nucleus. For example, this can be necessary for anions, rydberg states, and long range interaction between two systems. Among all available basis set with diffuse and polarization functions, let's mention the widely used augmented correlation consistent polarized triple, quadruple, and quintuple zeta basis sets (aug-cc-pVnZ, with n = T, Q, 5).

In addition of diffuse functions, one can use bond functions [46] in the calculation of weakly bound species such as van der Waals systems. In the latter, the bond length is significantly large and using bond functions, i.e. functions positioned at the middle of the bond length, can improve the description of the interaction energy.

It is usual to compute the intermolecular interaction energy as the difference between the energy of the system and the energy of the isolated molecules (supermolecular approximation). One disadvantage of this method is the impossibility to distinguish between the various intermolecular interactions (electrostatic, induction and dispersion). Moreover, the calculated bond energy is overestimated. This error is called basis set superposition error (BSSE) and arises from the fact that each monomer is better described in the interacting dimer than when it is isolated. A method commonly used to minimize the BSSE is the counterpoise method of Boys and Bernardi [47]. This consists in the calculation of the energy change of each monomer when the basis set change from the monomer basis set up to the dimer basis set.

The ab-initio calculations give in general very reliable results for the energy, especially for small molecules. We want to close this section with a comparison of the relative accuracy of the results for different singleconfiguration methods presented by Young in [39]:

HF MP2 < CISD ∼ = MP4 ∼ = CCSD < CCSD(T) < CCSDT < Full CI

Scattering formalism

The time independent approach of quantum scattering theory was developed in the second half of the last century. The formulation proposed by Arthurs and Dalgarno [10] or Curtiss and Adler [48] for atom-diatom collisions are the most commonly used [49]. Subsequent generalizations of their methodologies have been done for larger systems [50] and for various open shell cases [51, 52]. The main differences between the different formulations come from the use of different frames: the body fixed (BF) or the space fixed (SF) coordinates.

Pack [53] and Alexander and DePristo [54] compared the scattering formulations using the SF and BF and concluded that both formulations have their own advantages. The BF is the natural framework to describe the potential energy surface and gives a simpler mathematical form of the rotational sudden approximations. The SF formulation expresses the scattering asymptotic boundary conditions more conveniently and does not require any change of coordinates to obtain the S matrix. In this section, we describe the scattering between two systems in the SF frame as it is the one we use and which is implemented in our scattering code.

Single-channel scattering

The main concepts of the scattering formalism can be introduced by analysing first the elastic scattering of two structureless particles. Our starting point is then the time-independent Schrödinger equation for single-channel scattering

- 2 2µ ∇ 2 + V (r) Ψ = EΨ, (2.11)
where µ is the reduced mass of the system, V (r) is the interaction potential depending of the relative position of the particles which we assume to vanish faster than r -1 for large value of r, and E the total energy. Before the collision, i.e. when the particles are infinitely separated, the interaction potential is zero and this equation reduces to the one of a plane wave directed towards the target. The Z axis is put along the incident direction and the incoming asymptotic solution is written:

φ(r) r→∞ = exp (ikz), (2.12)
where k is the wavevector associated with the relative motion such as E = k 2 2 2µ . After the collision, the cylindrical symmetry is lost as all the scattering directions are now possible and the outcoming asymptotic solution at infinite separation becomes now a spherical wave:

ϕ(r) r→∞ = exp (ikr) r . (2.13)
Taking into account that the amplitude of the outgoing wavefunction varies with the scattering direction, a peculiar solution called the stationary scattering wave function (SWF) of the Schrödinger equation satisfies these asymptotic conditions:

Ψ(r) r→∞ = A exp ikz + f k (θ, ϕ) exp (ikr) r .
(2.14)

In the expression 2.14, the scattering amplitude f k (θ, ϕ) contains all the information about the collision process while the constant A is arbitrary and independent of r and the SWF is not square integrable .

The cross section is the measurable quantity characterising a collision processes. The differential cross section is defined as the outgoing flux of particles scattered through the spherical surface r 2 dΩ divided by the incident flux [38]. The flux density is defined as usual as

J = 1 2µ (Ψ * ∇Ψ -Ψ∇Ψ * ). (2.15)
The flux density for the incident and for the scattered wave functions are easily shown to be

J in = k µ , (2.16
)

J sc = k |f k (θ, ϕ)| 2 µr 2 .
(2.17)

The differential cross section is then dσ dΩ = |f k (θ, ϕ)| 2 , and the total cross section is obtained by integration over the scattering angles,

σ tot = π 0 2π 0 |f k (θ, ϕ)| 2 sin θdθdϕ.
(2.18)

It is worth noting that these expressions of the cross sections are independent of the choice of the arbitrary coefficient A.

Partial wave analysis

For a central potential, i.e. a potential depending only on r, the Hamiltonian commutes with both the angular momentum operators L 2 and L z . The scattering wave function can then be expended in spherical harmonics.

Ψ = 1 r ∞ l=0 m u lm (k, r)Y m l (θ, ϕ). (2.19)
The Schrödinger equation 2.11 is satisfied by the u lm functions and is independent of m.

- d 2 dr 2 -k 2 + l(l + 1) r 2 + V (r) u lm (k, r) = 0. (2.20)
Therefore u lm is usually replaced by u l . We can now find the boundary conditions respected by this wave function by looking again at the problem of the free particle (i.e. when V (r) = 0)

The previous equation reduces to

- d 2 dr 2 -k 2 + l(l + 1) r 2 u 0 l (k, r) = 0. (2.21)
There are two well known sets of linearly independent solutions of 2.21: the (zj l (z) ,zn l (z)) and (zh

(1) l (z) ,zh (2)
l (z)) sets, with z = kr. j l (z) and n l (z) are the spherical Bessel and Neumann functions while zh The free particle wavefunction can be expressed as a linear combination of any of these two sets as for example u 0 l = A 0 l zj l (z) + B 0 l zn l (z). If we impose the condition that u l (r = 0) = 0, we obtain B 0 l = 0 and thus:

u 0 l (r) ≈ A 0 l (k)(kr)j l (kr) = A 0 l (k) sin kr - lπ 2 . (2.22)
If V (r) tends towards zero faster than 1/r 2 when r → ∞, then the equation 2.20 is asymptotically identical to 2.21. The asymptotic solutions of the equation 2.20 can be written

u l (r) r→∞ ≈ kr [A l (k)j l (kr) + B l (k)n l (kr)] .
(2.23)

Which can be written as,

u l (r) r→∞ ≈ A l (k) sin kr - lπ 2 + B l (k) cos kr - lπ 2 .
(2.24)

or

u l ≈ C l (k) sin kr - lπ 2 + δ l (k) , (2.25)
where the phase shift δ l is

tan δ(k) = B l (k) A l (k) . (2.26)
The term in n l (kr) appears due to the action of the potential, and results in a phase shift with contains the information of the collision. When the potential is zero for all value of r, the phase shift is also zero. If the potential is attractive, the phase shift is positive and the radial wave function is pulled in with respect to the free radial wave function. In the case of a repulsive potential the phase shift is negative and the radial wave function is pushed out with respect to the free radial wave function.

The phase shift δ l can be related with the cross section by rewriting the first term in equation 2.14 as exp (ikz) = exp (ikr cos θ) = (4π)

1 2 l i l (2l + 1) 1 2 j l (kr)Y l0 ( r) = l i l (2l + 1)j l (kr)P l (cos θ), (2.27)
which asymptotically is exp (ikz) = l i l (2l + 1) sin kr -lπ 2 kr P l (cos θ).

(2.28)

Because we have chosen k to lie along the Z axis, the scattering amplitude f k (θ, ϕ) depends only on θ. Then we can expand the scattering amplitude in Legendre polynomials,

f k (θ) = ∞ l=0 f l (k)P l (cos θ). ( 2
i l (2l + 1) sin kr -lπ 2 kr + f l (k) exp (ikr) r = C l (k) r sin kr - lπ 2 + δ l (k) . (2.30)
After some algebraic work on the later expression, and grouping the term with the same sign in the exponents one obtains the expressions:

C l (k) = i l (2l + 1) k exp(iδ l ), (2.31)
f l (k) = 2l + 1 2ik (exp (2iδ l ) -1) . (2.32)
The factor S l = exp(2iδ l (k)) is called the scattering matrix element. The transmission matrix is defined by T l (k) = 1 -S l (k). Both S l (k) and T l (k) are complex quantities, but a real matrix can also be defined which is called the reactance matrix K l = tan(δ l (k)). The T and K matrices are related by

T l = -2iK l (k)[1 -iK l (k)] -1 .
The asymptotic radial wave function can be written as a function of K l as

u l (r) ≈ sin kr - lπ 2 + K l (k) cos kr - lπ 2 .
(2.33)

The scattering amplitude 2.18 is readily obtained to be

f k (θ) = 1 2ik ∞ l=0 (2l + 1) [exp (2iδ l ) -1] P l (cosθ) (2.34)
which give the following expression for the cross section as a function of the T matrix:

σ tot (k) = l σ l = π k 2 l (2l + 1)|T l (k)| 2 , (2.35)
or as a function of the phase shift:

σ tot (k) = l σ l = 4π k 2 l (2l + 1) sin 2 (δ l (k)).
(2.36)

Resonances

The variation of the cross section as a function of collision energy is usually quite slow. However, this variation can become rapid in a certain energy interval. These peaks of the cross section are called resonances. Some of them can be very simply related with the shape of the effective potential and for this reason they are called shape resonances. Indeed, equation 2.20 can be rewritten as a function of an effective potential made of the sum of the potential plus the centrifugal term. The centrifugal term is repulsive for all values of r while the potential can be repulsive in some regions and attractive in some others. For this reason at intermediate distances, a local minimum can be found which is followed by a barrier for larger distances.

The particles can then be trapped for some time between the repulsive barrier at short range and the centrifugal barrier before tunnelling out. These states are also known as quasi-bound levels.

When the scattering energy crosses the energy of a quasi-bound state, the resonant phase shift δ l increases of π, and δ l = π(n + 1/2) where n is natural number. In this domain of energy, we can expand cot (δ l (k)) in the vicinity of the resonance energy E r ,

cot (δ l (k)) = cot (δ l )| E=Er +(E-E r ) d(cot (δ l )) dE E=Er +O[(E-E r ) 2 ]. (2.37)
The first term is zero and we can define Γ as:

d(cot (δ l )) dE E=Er = - 2 Γ . (2.38)
If Γ is very small, the variation of δ l (k) is very fast. Rewriting equation 2.36 in function of cot (δ l (k)) and using 2.37 and 2.38, we get for the resonant partial wave contribution to the cross section near the resonance:

σ l = 4π k 2 (2l + 1)(Γ/2) 2 (E -E r ) 2 + (Γ/2) 2 .
(2.39)

This expression is known as the Breit-Wigner form. As the others phase shifts vary slowly in this domain of energy, the peak of σ l causes a rapid variation of the total cross section. Another frequently met type of reso-nance is known as a Feshbach resonance, and will be briefly described in the next subsection.

Collision of two diatomic molecules

We can now generalize the single channel scattering to the multichannel case. The Schrödinger equation for the collision between two linear rigid rotors, in the SF coordinates located at the center of mass of the system, can be written as [11]

H( R 1 , R 2 , R) -E Ψ( R 1 , R 2 , R) = 0, (2.40)
where

H = H 1 ( R 1 ) + H 2 ( R 2 ) + V ( R 1 , R 2 , R) + T(R). (2.41)
R are the coordinates of the vector linking the centers of mass of the two molecules, R i are the angular coordinates of the molecule i, V is the interaction potential, E is the total energy of the system, T is the intermolecular kinetic operator and H i is the rigid rotor Hamiltonian of the diatomic molecule i. The Schrödinger equation for each isolated diatomic molecule is

H i Y j i m i ( R i ) = E i Y j i m i ( R i ), (2.42)
with E i = B i j i (j i + 1) and B i is the rotational constant. The rotational states of |j 1 m 1 and |j 2 m 2 of the diatoms can be changed by the collision process while the quantum numbers associated with the total angular momentum and its projection along the Z space fixed axis are conserved. The total angular momentum J is obtained by first coupling the two rotors angular momenta: j 12 = j 1 + j 2 , and by coupling the resulting angular momentum with the one associated with the relative angular momentum of the two rotors: J = j 12 + l. A similar procedure is followed to define an angular basis set which is a set of eigenfunctions of J and J z :

Φ JM j 1 j 2 j 12 l ( R 1 , R 2 , R) = m j 1 m j 2 m j 12 m l j 1 m j 1 j 2 m j 2 |j 12 m j 12 j 12 m j 12 lm l |JM Y j 1 m 1 ( R 1 )Y j 2 m 2 ( R 2 )Y lm ( R). (2.43)
We use the latter basis set to expand the wave function of the system:

Ψ( R 1 , R 2 , R) = 1 R JM j 1 j 2 j 12 l G JM j 1 j 2 j 12 l (R)Φ JM j 1 j 2 j 12 l ( R 1 , R 2 , R). (2.44)
Substituting this expansion in equation 2.40, multiplying by 2.43 on the left, integrating over all angular coordinates and using γ ≡ j 1 j 2 j 12 l, we get the radial coupled equations:

d 2 dR 2 + k 2 γ - l(l + 1) R 2 G Jγ (R) = 2µ 2 γ W J γγ (R)G Jγ (R), (2.45)
where

W J γγ (R) = d R 1 d R 2 d RΦ J γ ( R 1 , R 2 , R)V ( R 1 , R 2 , R) × Φ J γ ( R 1 , R 2 , R).
(2.46)

and where

k 2 γ = 2µ 2 (E -E 1 -E 2 ).
The quantum number M has been dropped as the equation 2.45 is independent of it.

For a given total energy, these equations are solved for each value of J giving a non-zero contribution to the cross section and the size of the rotational basis set is varied up to convergence. If no other approximation is used one speaks of close coupling calculations. The potential is expanded in an angular basis set similar to the one proposed by Green [11]

V ( R 1 , R 2 , R) = l 1 l 2 l A l 1 l 2 l (R)I l 1 l 2 l (( R 1 , R 2 , R) (2.47)
I l 1 l 2 l ( R 1 , R 2 , R) = m 1 m 2 m l 1 m 1 l 2 m 2 |lm Y l 1 m 1 ( R 1 )Y l 2 m 2 ( R 2 )Y * lm ( R). (2.48)
The angular dependence of the interaction potential can equivalently be expressed as a function of body fixed angles i.e in the body fixed frame which is defined by setting the body fixed Z-axis along the collision coordinates R = (Θ, Φ) or in other words by using the angles that the axis of the target and of the projectile make with the intermolecular axis, θ 1 and θ 2 such as cos(θ i ) = R i • R and the associated dihedral angle ϕ = ϕ 1ϕ 2 .

The two frames being simply related by the Euler rotation R(Φ, Θ, 0), it is straightforward to obtain the following expression:

I l 1 l 2 l ( R 1 , R 2 , R) = 2l + 1 4π 1 2 l 1 l 2 00|l0 P l 1 0 (θ 1 )P l 2 0 (θ 2 ) + m (-1) m 2 l 1 l 2 m -m|l0 P l 1 m (θ 1 )P l 2 m (θ 2 ) × cos(mϕ) . (2.49)
This body fixed basis set is the most commonly used in the analytical models of potential energy surfaces developed for collision between two linear rigid rotors. The coupling matrix elements W J γγ (R) can be obtained by integrating over all angular coordinates:

W J γγ (R) = l 1 l 2 l A l 1 l 2 l (R)(-1) J+j 1 +j 2 +j 12 (4π) -3 2 ×([l] 2 [l 2 ][l 2 ][j 1 ][j 2 ][j 12 ][L][j 1 ][j 2 ][j 12 ][L ]) 1 2 l L L 0 0 0 l 1 j 1 j 1 0 0 0 × l 2 j 2 j 2 0 0 0 L L l j 12 j 12 J      j 12 j 2 j 1 j 12 j 2 j 1 l l 2 l 1      . (2.50)
The matrix elements W J γγ (R) are independent of M and vanish for different initial and final values of the spectroscopic parity, P = (-1) J+j 1 +j 2 +L . This makes also P to be a good quantum number.

The radial part of the scattering wave function is propagated from the classically forbidden region, where it is initialised to zero, up to an intermolecular distance large enough for the interaction potential to be negligible and where it is requested to satisfy the boundary conditions of the scattering process. These conditions are presented below for the open channels (for k γ > 0) by generalizing the equation 2.33 and writing it in a matrix form proposed by Launay [12]

G J (R) = J J (R) + K J N J (R), (2.51)
where J J γj 12 lγ j 12 l (R) = δ γγ δ j 12 j 12 δ ll k 1/2 γ Rj l (k γ R) (2.52)

N J γj 12 lγ j 12 l (R) = δ γγ δ j 12 j 12 δ ll k 1/2 γ Rn l (k γ R). (2.53)
Recalling that T = -2iK(1 -iK) -1 , we can write the partial cross section averaged over the initial states as

σ J γ→γ = π(2J + 1) k 2 γ (2j 1 + 1)(2j 2 + 1) j 12 l j 12 l |T J γj 12 lγ j 12 l | 2 , (2.54)
and the total cross section is

σ γ→γ = ∞ J=0 σ J γ→γ .
(2.55)

Before closing this section, we can mention the Feshbach resonances. These resonances, which are associated with the coupling between the open and the closed channels, happen when the collision energy of a given open channel allows the near opening of a closed quasi-bound state.

A symmetric top molecule in collision with an atom

The collision of a symmetric top molecule with an atom was presented long ago by Green [55]. It is shortly overviewed here as the same formalism will be used in the Chapter 5 dedicated to the collision between an atom and a triatomic rigid bender. The total Hamiltonian of the system is written again in SF coordinates:

H = H rot ( Ω) + T(R) + V ( Ω, R), (2.56)
where R is the collision coordinate, T(R) is the relative kinetic energy and the rotor orientation is defined by the Euler angles Ω = (αβγ) that rotate the SF axis to the principal axis of inertia of the symmetric top. The BF frame is again defined by the rotation aligning the SF z-axis along the intermolecular coordinate.

The term H rot in 2.56 is the symmetric top Hamiltonian. The eigenfunction of this Hamiltonian can be labelled by j, m and k which are the quantum numbers associated with the total rotor angular momentum and its projections over the SF z-axis and the BF z -axis. The latter eigenfunction is proportional to the Wigner matrix associated with the Euler rotation R(αβγ):

|jkm = 2j + 2 4π 2 D j * km (αβγ).
(2.57)

The eigenvalues of H rot for a symmetric top, in which the z -axis is chosen as the symmetric axis and I 1 and I 3 are the moment of inertia about the principal axis x and z , are

E jk = 2 1 2I 1 j(j + 1) + 1 2I 3 - 1 2I 1 k 2 .
(2.58)

The degeneracy in k and -k is used to define a symmetrized linear combination of the symmetric top eigenfunctions:

|j kmp = 1 √ 1 + δ k0 |j km + (-1) p |j -km , (2.59)
where k is the absolute value of k and p = 0 or 1 except for k = 0 where only p = 0 is possible. The parity of these symmetrized functions is (-1) (j+p) . The formalism presented here can be straightforwardly extended to the case of an atom colliding with an asymmetric top because the asymmetric top wave functions are linear combinations of the symmetric top ones.

The interaction potential is conveniently expressed in terms of spherical harmonics in the SF frame,

V (R , Θ , Φ ) = λµ υ λµ (R )Y λµ (Θ , Φ ).
(2.60)

Using the rotation which links the SF and BF coordinates, we can write the same expansion as a function of body fixed angles

V ( Ω, R) = λµν υ λµ (R)D λ µν ( Ω) * Y λµ ( R).
(2.61)

As spherical harmonics are complex, it is more convenient to use instead the following real symmetrized combinations:

V (R , Θ , Φ ) = λ≥µ≥0 υ λµ (R ) (1 + δ µ0 ) Y λµ ( R ) + (-1) µ Y λ,-µ ( R ) . (2.62)
The rotational angular momentum is coupled with the orbital angular momentum l in order to obtain eigenfunctions of J 2 and J Z .

|JM jkl = m jmlm l |JM |jkm |lm l , (2.63)
where |lm l = Y lm l (Θ, Φ) and all the projections are taken along the SF z-axis. The scattering wave function is expanded in the latter angular basis set

Ψ JM jkl = 1 R j k l G JM jkl j k l (R)|JM j k l .
(2.64)

The resulting coupled equations can be written

d 2 dR 2 + k 2 j k l - l(l + 1) R 2 G JM jkl j k l (R) = 2µ × j k l JM j k l |V |JM j k l G JM jkl j k l (R), (2.65) with k 2 j k l = 2µ 2 (E + E jk -E j k )
, where E is the kinetic energy and E jk is given by 2.58. The coupling matrix elements are

JM jkl| λµ υ λµ (R )Y λµ ( R )|JM j k l = λµ υ λµ (R)(-1) j+j +k-J × [j][j ][l][l ][λ] 4π l l λ 0 0 0 j j λ k -k µ (2.66)
where the notation [i] = 2i + 1 is used. As can be seen, the matrix elements are independents of M and there is no coupling between different value of J.

If we use the symmetrized symmetric top eigenfunctions 2.59 these matrix elements become

υ λµ Jjksl|[Y λµ + (-1) µ Y λ-µ ]/(1 + δ µ0 )|Jj k s l = υ λµ (1 + (-1) s+s +j+j +λ+µ ) 2 (1 + δ k0 )(1 + δ k 0 ) (ω Jjkl|Y λ±µ |Jj k l +(-1) s Jj -kl|[Y λµ |Jj k l ) , (2.67)
the first term contributes only if µ = ±(kk), and

ω = 1, k -k ≥ 0 (-1) µ , k -k < 0.
(2.68)

The second term in 2.67 contributes only if µ = k + k.

Bound states calculations

The bound states we are interested in are those which are associated with the trapping of the particles by the interaction potential. Such states exist only if the potential is attractive enough. Efficient variational methods (based on DVR/FBR methods) can be used to perform bound states calculations. However, we will present here another approach based on the use of the scattering theory. It is based on the fact that the coupled equations needed for scattering calculations are identical to those for bound states, the only difference being the applied boundary conditions. For large values of the collision coordinate R, the bound state wave function must vanish whereas the scattering wave function keeps oscillating. This method is, for example, coded in the BOUND program developed by Hutson [56] for the log-derivative propagator and is described in great detail in the corresponding paper. A first propagation is done outward starting from the classically forbidden region using the usual nullity initial condition. It is stopped in the region of the well at a value R = R Bound and the corresponding R-matrix R F orward (R Bound ) or log-derivative matrix Y F orward (R Bound ) is stored. A downward propagation is then performed starting in the asymptotic region and using WKB initial conditions down to the same value of R where the R-matrix R Backward (R Bound ) or log-derivative matrix Y Backward (R Bound ) is obtained. The continuity of the two kinds of propagated wave functions and their derivatives is then imposed by diagonalising the matrix equal to the difference between the two R-matrices or log-derivative matrices. The total energy is varied from the bottom of the well up to zero and the bound state energies are located by looking at the change in sign of the eigenvalues.

The convergence of the propagation method is a function of the step size h. The error for the usual R-matrix and log-derivative propagators are respectively proportional to h 2 and h 4 . In the case of the log-derivative propagator, the energy using a step h is

E(h) ≈ E(h = 0) + E (4) h 4 .
(2.69)

The calculations are then done for two values of the step size, h 1 and h 2 giving the energies E 1 and E 2 . The converged energy can be estimated from a Richardson extrapolation:

E(h = 0) = E 1 - h 4 (E 2 -E 1 ) h 4 2 -h 4

Introduction to the study

Among the molecules containing sulfur detected in the ISM, we will focus our attention in this study on one of the simplest: carbon monosulfide [57]. This molecule has been extensively observed not only in dense molecular clouds [58-60] but also in diffuse interstellar clouds [61]. CS is one of the molecules used commonly as a tracer to estimate the density in several interstellar regions [62].

The first theoretical study of the collisions of CS with H 2 was done by Green and Chapman [8] using a PES adapted from the Gordon and Kim electron gas model. Later, Turner et al [63] extended their work by including more rotational states in the calculation, using the same PES. Recently the collision of CS with He was studied by Lique et al [64, 65] using a PES computed at coupled cluster level. The comparison of their results with those of Green and Chapman [8] and Turner et al [63], based on the relative mass ratio, shows remarkable differences. Astrophysicists often use scaled He collisional rate coefficients to approximate excitation by H 2 . The use of this approximation has previously been found to be inaccurate for other systems like SiS [66, 67] or SO [68]. A new study of the collisions of CS with H 2 was then necessary.

Close coupling calculations were performed using the rigid rotor approximation for both H 2 and CS following the lines presented in the previous chapter. The ab-initio energies for the CS-H 2 complex were determined at CCSD(T) level, using aug-cc-pVQZ basis set with bond functions. These calculations were done for a grid of points defined in the space of the four coordinates presented in Fig. 3.1. All the calculations were made using the Molpro package [40].

The computed ab-initio energies were fitted to the analytical function previously used for the H 2 -HF system [14]. This function has the general form,

V (R, θ 1 , θ 2 , ϕ) = 6 l 1 =0 12 l 2 =0 min(l 1 ,l 2 ,2) m=0 f l 1 l 2 m (R) × P m l 1 (cos θ 1 ) P m l 2 (cos θ 2 ) cos(mϕ), (3.1)
where l 1 is restricted to even values and P m l (cos θ) are associated Legendre polynomials. In a first step, ab-initio data was fitted to the angular function for each of the R values belonging to the grid, using a least squares procedure. The radial coefficients were then interpolated using the cubic spline method. The extrapolation was performed using two different radial functions for the short range and the long range parts which take the following form,

f SR l 1 l 2 m (R) = τ m l 1 l 2 exp(-γ m l 1 l 2 R), f LR l 1 l 2 m (R) = β m l 1 l 2 R -α m l 1 l 2 . (3.2)
The parameters of the long range part were adjusted to fit the ab-initio grid in the asymptotic region while the β m l 1 l 2 and α l 1 l 2 parameters were determined by requiring the continuity of the radial functions and their first derivatives at the smallest value of R computed.

The equilibrium structure of the CS-H 2 complex is found to be linear, with the carbon atom pointing toward H 2 . The bound levels of the CSpara-H 2 and CS-ortho-H 2 complexes have been computed. This is the first prediction of these bound levels. We then computed the cross section involving 16 rotational levels of CS. These cross sections were also used to determine the rate coefficients. In the joined publication we compare the present results with those previously published for the collisions of CS with para-H 2 . We report also the first cross sections and rate coefficients for the collisions with ortho-H 2 . The approximation using the square root of the relative mass of the colliders to obtain the rate coefficients of the collision between a molecule and H 2 from the scaling of the data available with He is found to be a good qualitative approximation for the CS-para-H 2 system. In contrast, the rate coefficients based on the crude electron gas model show strong discrepancy with the new results and have to be scaled by a factor of 2 to get a good qualitative agreement.

Publications

I. INTRODUCTION

Carbon monosulfide (CS) is the first molecule containing sulfur that has been detected in the interstellar space. 1 Owing to its strong emission, it is used as a tracer of molecular gas in various regions of the interstellar medium in our galaxy and in external galaxies. The scientific applications may range from determining sulfur isotopic ratios in the galactic disk 2 to study the anatomy of a filament 3 or tracing physical conditions in protostellar envelopes or in hot cores. 4 Analysis of observed spectra requires the knowledge of accurate spectroscopic data such as frequencies and Einstein coefficients when the molecules are in local thermal equilibrium (LTE), i.e., when the local gas density is above a critical density. When LTE does not apply the analysis of spectra must consider collisional processes that lead to population or depopulation of the energy levels of the observed molecule. The physical conditions where CS is most commonly used as a tracer correspond to regions where H 2 is the most abundant collider and therefore the collisional rate coefficients of CS by H 2 are included in the radiative transfer analysis of spectra.

The only available rate coefficients for the excitation of CS with H 2 are those of Turner et al. 5 Their dynamical calculations were based on the crude potential energy surface a) Electronic mail: t.stoecklin@ism.u-bordeaux1.fr.

(PES) of Green and Chapman 6 adapted from a Gordon and Kim 7 electron gas model for CS-He, treating the CS molecule as a rigid rotor. For an intermolecular separation R > 8 a o , the lowest 4 terms of the PES expansion were modified to join smoothly with asymptotic electrostatic interaction 8 of CS-H 2 . Turner et al. extended the dynamical calculations of Green and Chapman and provided rate coefficients calculations for transitions between the first 21 rotational levels of CS and for temperature between 20 K and 300 K. Those rate coefficients have been used very recently to analyse hot cores. 4 Astrophysicists often use scaled He collisional rate coefficients to mimic excitation by para-H 2 (pH 2 ) and ortho-H 2 (oH 2 ) as it has been done recently in the study of the molecular abundances of the inner layers of the carbon star IRC + 10216 (Ref. 9) where the recent rotational and rovibrational rate coefficients of Lique and Spielfiedel 10,11 were used for temperature up to 1500 K. While He and pH 2 (j = 0) results are rather similar for heavy diatomic molecules, several recent results for neutral diatomic and polyatomic molecules [12][13][14][15] show that scaled He calculations cannot mimic collisions with oH 2 (j = 1) as the latter have much larger rate coefficients. Recent results on H 2 O-H 2 (Ref. 16) show that rate coefficients for collisions with pH 2 (j = 0, 2, 4) are found to be of the same order of magnitude as rate coefficients for collisions with oH 2 (j = 1) at temperature where those levels are populated and that the influence of the rotational level of H 2 is much less pronounced for temperatures above 800 K. It is likely that we can generalize this behaviour to diatomic molecules and therefore the use of either scaled-He or pH 2 (j = 0) calculations might be not adequate in regions corresponding to temperatures in the range of 200-800 K, which is quite a large part of the hot core regions and AGB stars envelops.

Considering the importance of CS for astrophysical applications, we have decided to provide a new state-of-the-art PES allowing the computation of collisional rate coefficients for the excitation of CS including the influence of the rotational level of H 2 (j = 0, 1, 2). As a first use of this new PES, we computed the rovibrational bound states energies of the H 2 -CS complex as well as the interaction second virial coefficient for the H 2 -CS pair. The experimental infrared spectrum of H 2 -CS is not known. This first theoretical work will not only help assigning the infrared spectrum when it will be measured but aims to provide a good test of the PES which is presented in this paper and which will be used to calculate rotational excitation cross sections.

II. THEORY

A. Potential energy surface

The interaction potential of CS with H 2 has been calculated within the supermolecular approach using the coupledcluster method with single and double excitations and a perturbative treatment of triple excitations (CCSD(T)) and using the augmented correlation consistent polarized valence quadruple zeta (aug-cc-pVQZ) basis set, 17 supplemented with a set of midbond functions. 18 The interaction energy, corrected for the basis set superposition error, was obtained by subtracting the sum of the energies of the isolated monomers H 2 and CS from the energy of the H 2 -CS complex, each monomer's energy being calculated with the full basis set used for the complex. This corresponds to the counterpoise procedure of Boys and Bernardi 19 applied to the rigid rotor case. All calculations were performed with the MOLPRO package. 20 The bond length of the H 2 monomer was set equal to the vibrationally averaged value in the rovibrational ground state, r o = 1.448736 a o . 21 For the CS monomer, the vibrationally averaged bond length is expected to be very close to the equilibrium bond length. Indeed, the vibrationally averaged bond length of CO (Ref. 21) is close to the equilibrium bond length: the difference is 7 × 10 -3 a o . For CS, this difference is expected to be even smaller because both the harmonic frequency and the anharmonicity constant are smaller. Therefore we simply use the equilibrium bond length r e = 2.900619 a o . 22 We used body fixed coordinates in which R, θ 1 , θ 2 , and ϕ are the coordinates which describe the relative orientations of the two molecules as shown in Fig. 1. The radial grid includes 73 points ranging from 6 to 30 a o along the R coordinate. By taking into account the exchange symmetry of the H atoms and the global inversion symmetry, we have reduced the intervals of variation of the angles θ 1 , θ 2 , and ϕ to [0, 90 • ], [0,180 • ], and [0, 90 • ], respectively. For all angles, steps of 15 • were used. Moreover, a supplementary set of points for R ranging from 35 up to 50 a o by step of 5 a o was calculated. These points, which have not been included in the fitting process, were used to check the accuracy of the long-range part of the fitted PES.

The functional form of the PES and the fitting procedure have already been described in a previous work. 23 Let us recall here only some of the essential features. The angular part of the functional form is represented by a product of normalised associated Legendre functions and a cosine function,

ȳm l1,l2 (θ 1 , θ 2 , ϕ) = P m l1 (cos θ 1 ) P m l2 (cos θ 2 ) cos(mϕ), (1)
with 0 ≤ l 1 ≤ 6 and restricted to even values, 0 ≤ l 2 ≤ 12, and 0 ≤ m ≤ 2. For each point of the radial grid, a least-square procedure has been performed to compute the coefficients of the development on the angular functions. Then the radial functions were obtained by cubic spline interpolation. Outside the range of the radial grid, the radial functions were extrapolated for shorter distances by exponential functions and for longer distances by reciprocal power functions.

B. Bound states calculations

The method we use is widely employed to calculate the rovibrational levels of atom-diatom or diatom-diatom 24 van der Waals complexes. It is based on the fact that the coupled equations needed for scattering calculations differ from those for bound states only in the applied boundary conditions. This method, which is particularly well suited to obtain the higher rovibrational bound levels which are difficult to converge using a variational method, is for example coded in the BOUND program developed by Hutson with the log-derivative propagator and is described in great detail in the corresponding paper. 25 We implemented the bound state calculations in our diatom-diatom close coupling scattering code according to the recommendation of Danby 26 for an R-matrix propagator as our program uses the Magnus propagator and tested it by reproducing the spectra of the H 2 -CO complex calculated by Jankowski and Szalewicz. 21 We previously used the same approach to calculate the bound states of several atomdiatom complexes [27][28][29] and the details of the implementation of the bound state calculations in our code can be found in Ref. 27. Briefly, our code solves the rovibrationally inelastic close coupling equations in the space fixed frame and uses the Jacobi coordinates. The vibrational levels of the H 2 -CS van der Waals complex are calculated by performing calculations for all the values of the total angular momentum J and parity ε leading to bound states. The rigid rotor calculations were performed separately for ortho-and para-H 2 and the two lowest rotational states of ortho/paraH 2 were included in each case. Fifteen rotational states were included in the basis set describing the CS molecule. We used the rotational constant B = 60.853 cm 1 for H 2 and B = 0.8200462 cm -1 for CS. The maximum propagation distance was 50 a o . The calculations were performed for two values of the propagator step size (0.1 and 0.05 a o ) and the values of the bound state energies were obtained from a Richardson extrapolation as suggested by Hutson. 25

C. Second virial coefficient

In order to provide data for a future test of the quality of our PES we have calculated the cross second virial coefficient B 12 (T) within the classical approximation. The value of B 12 for a given temperature T was calculated with the formula 30

B 12 (T ) = N A 4 2π 0 dϕ π 0 dθ 1 sin θ 1 π 0 dθ 2 sin θ 2 ∞ 0 R 2 dR × 1 -exp -V (R, θ 1 , θ 2 , ϕ) k b T (2) in which V (R, θ 1 , θ 2 , ϕ)
is the potential energy and k b and N A are, respectively, the Boltzman and Avogadro constants. This integral was evaluated numerically, using a Gauss-Chebyshev quadrature of 20 points for the integration over ϕ and a Gauss-Legendre quadrature with 20 points over θ 1 and θ 2 . The radial integration was performed by the Romberg method, with the integrand approximated by 1 for R < 3 √ 3 Å and by 0 for R > 100 Å.

III. RESULTS AND DISCUSSION

A. Potential energy surface

We found that the equilibrium structure of the H 2 -CS van der Waals system is linear, with the carbon atom pointing toward H 2 . This configuration is similar to the one of the isovalent system H 2 -CO. 21 The energy and the geometry associated with the minimum are, respectively, D e = 173 cm -1 , internuclear axis is rotated while the H 2 molecule is aligned along the intermolecular axis. In Figs. 3 and4, the H 2 internuclear axis is rotated while the CS molecule is aligned along the intermolecular axis with the S atom (Fig. 3) or the C atom (Fig. 4) pointing toward H 2 . In Fig. 5, both H 2 and CS internuclear axis are rotated, while the intermolecular distance R is relaxed and the dihedral angle ϕ is fixed to 0 • . We can see two equivalent principal minima in Figs. 4 and 5 as a consequence of the symmetry of permutation of hydrogen atoms.

θ 1 = 0 • (or 180 • ), θ 2 =
The principal minimum is also seen in Fig. 2. No secondary minimum has been found. The minimum point which is seen in Fig. 3 is actually a saddle point as shown by Fig. 5. In this last figure, we can observe that a high potential ridge separate the two equivalent minima. Consequently, the most favorable pathway to convert from one minimum to the other goes through a combined rotation where H 2 do a half turn while CS completes a full turn. Let us note that displacing ϕ from 0 • should not change strongly this topology because the potential energy is independent of ϕ in all the configurations where at least one diatomic axis is aligned along the intermolecular axis.

The quality of the fitted PES has been checked by evaluating the root-mean-square (rms) of the differences between the ab initio energies and the fitted energies. We have divided the range of R spanned by the ab initio data in three regions: , respectively. This can be compared to the rms of the potential energies in the same regions which are 704.2, 25.5, and 1.28 cm -1 , respectively. The FOR-TRAN code of the PES is available from the authors upon request.

B. Rovibrational bound states

The results of the bound states calculations for the pH 2 -CS and oH 2 -CS complexes with J ≤ 2 are presented in Tables I andII, respectively. All other energy levels with J ≥ 3 are available as supplementary material. 31 The parameters of the scattering calculations (i.e. the basis set, the propagator step size) ensure that an accuracy of about 10 -2 cm -1 is reached for the energy of bound states. However, because the energy gap between some states can be as small as 10 -4 cm -1 , the search of the sign change in the difference between inward and outward R-matrices 28 was performed with a 10 -6 cm -1 threshold. The maximum value of the total angular momentum J leading to bound states was found to be 10 for pH 2 -CS and 12 for oH 2 -CS. The rovibrational energies are given relative to the ground state energy of infinitely separated CS and pH 2 . Hereinafter we denote by j 1 and j 2 the rotational angular momentum of H 2 and CS respectively. The spacing between the levels j 1 = 0 and j 1 = 2 in the para form and between the j 1 = 1 and j 1 = 3 levels in the ortho form are significantly larger than the well depth of the potential. Therefore all the ortho/para-H 2 -CS bound states wavefunctions are predominantly expanded on the ortho/para-H 2 rotational ground state. In Table I, we report the energy, the quantum numbers J and ε, and the approximate quantum numbers L and j 2 TABLE I. Calculated rovibrational bound states of pH 2 -CS for J ≤ 2. For each state, we report the energy in cm -1 , the total rotational quantum number J, the parity ε, the CS rotational quantum number j 2 , the orbital quantum number L, and the percentage weight (w) of the leading basis set function. For some states, several basis functions need to be given in order to distinguish them from lower states with same J and ε.

State

State State associated with each pH 2 -CS bound state. The value of j 1 is assumed to be 0 for all the states. We report also the relative weight, given in percent, of the dominant basis set functions in the rovibrational wavefunction. The relative weights have been obtained by multiplying by 100 the square of the largest expansion coefficients. In most cases good quantum numbers could be assigned to a given para bound state. There are, however, a few exceptions for J = 0, 1, 2, 3, 4 where a single attribution was not possible. In Table II, we report only the energy and the quantum numbers J and ε for each oH 2 -CS bound state, while j 1 is assumed to be 1. It was not possible to assign more quantum numbers because we did not find any angular expansion coefficients of the wavefunction larger than 0.2. Thus the oH 2 -CS spectrum appears considerably more complex than the pH 2 -CS spectrum. The total number of bound states supported by the present PES is 186 for pH 2 -CS and 553 for oH 2 -CS. There are about two times more states for the ortho form than for the para form. This is related to the fact that with j 1 = 0, a single value of j 12 is possible for any given value of j 2 , while with j 1 = 1, three values of j 12 are obtained. The same ratio was observed previously for the H 2 -CO complex. 24 However, the number of bound states is considerably larger in the case of H 2 -CS, about five times more, thereby suggesting a much higher complexity of the infrared spectra. The dissociation energies of the ground states of the para and ortho species are found to be only 35.87 cm -1 and 49.88 cm -1 , respectively. Thus the energies of both ground states are a significative part of the potential well depth, simi-larly to the case of H 2 -CO. Let us note that the oH 2 -CS ground state is found for J = 2, while it is J = 1 for oH 2 -CO. 24 In the case of pH 2 -CS, for L = 0, the spacing between bound states associated with two successive values of j 2 are reported in Table III. As it can be seen, the spacings are very close to the rotational spacing of the free CS diatom. This suggests that the CS monomer behaves like a weakly hindered rotor. Indeed, Fig. 5 shows that the potential energy is below the energy of the ground state of pH 2 -CS in almost the whole coordinate space. Again, this is similar to the dynamics of the H 2 -CO system. 21 We also report in Table IV, for L = J, the energy spacings between the para levels associated with successive values of L. If we assume that the energy levels are those of a free IV shows that these values are distributed in a small range around an average value of 0.47 cm -1 . Furthermore, if we consider the system H 2 -CS at the equilibrium intermolecular separation R e = 8.601 a o , we can calculate the value of B rel with the equation B rel = ¯2/2μR 2 e . We obtain a value of 0.42 cm -1 which is in qualitative agreement with the average value deduced from the energy spacings. This indicates that the respective orbital motion of the two monomers is close to a free motion, and therefore, the orbital motion is weakly interacting with the rotation of CS. Again, this leads to conclude that the CS monomer is almost freely rotating.

J ε j 2 L w Energy J ε j 2 L w Energy J ε j 2 L w Energy 0 + 0 0 - 35

C. Second virial coefficient

The calculated cross second virial coefficient B 12 (T) is shown in Fig. 6. There are unfortunately no experimental data available for a mixture of H 2 and CS which could be compared with our results. Whenever available, they will provide a good test of the accuracy of our surface. As a matter of fact, we found recently 32 that the B 12 based scaling of the H 2 -CO potential energy surface suggested by Jankowski and Szalewicz 24 was crucial to get good agreement between theoretical and experimental inelastic cross sections. The quantum corrections of B 12 where not calculated here since these small corrections are expected to be significant only at low temperatures. 33

IV. CONCLUSION

A new four-dimensional analytical PES based on a large grid of high-level ab initio calculations was obtained for the H 2 -CS van der Waals complex. The single equilibrium structure of the H 2 -CS complex was found to be linear with the carbon pointing toward H 2 . The corresponding well depth is -173 cm -1 . Within the rigid-rotor approximation, all the rovibrational bound states of the H 2 -CS complex have been calculated. The dissociation energies of the ground states of the para and ortho species are found to be only 35.87 cm -1 and 49.88 cm -1 , respectively, indicating the energies of both ground states are a major part of the binding potential energy.

The PES supports a large number of bounds states, 136 for the para form and 553 for the ortho form. Rotational quasi-quantum numbers L and j 2 have been unambiguously assigned to most of the bound states of pH 2 -CS. For the states of oH 2 -CS, this assignment have not been possible, owing to the strong mixing between basis set functions. This higher complexity of the oH 2 -CS bound states results from the significantly larger number of bound levels, which is a consequence of the fact that j 1 = 1 allows to build two more states than j 1 = 0.

In the case of the pH 2 -CS system, a simple examination of some levels spacing indicate that the CS monomer behaves like a weakly hindered rotor. This is a consequence of the large zero point energy which allows the pH 2 -CS system the move freely in almost the whole angular coordinates space.

I. INTRODUCTION

The emission spectrum of the CS molecule has been extensively observed in many interstellar objects since its first detection in 1971. 1 The large number of rotational transitions detected is used to trace the different regions of circumstellar envelopes as detailed for example in a recent study dedicated to the carbon star IRC +10216. 2 The analysis of the observed spectra requires the knowledge of accurate spectroscopic data such as radiative and collisional rates with the most abundant atomic and molecular species of the interstellar medium: He and H 2 . The rate coefficients for the CS-He system are now well documented after the recent work by Lique et al. 3 for the first 31 rotational levels of CS and for temperatures ranging from 10 to 300 K using a potential energy surface calculated at the coupled-cluster level. On the contrary, the accuracy of the rate coefficients available for the collisions of CS with H 2 can be questioned. The old calculations of Green and Chapman, 4 performed at the coupled states level, were based on a semi empirical potential energy surface (PES) adapted from a Gordon and Kim 5 electron gas model for CS-He. More recently, using the same semi-empirical surface, Turner et al. 6 extended the dynamics calculations of Green and Chapman and provided rate coefficients calculations for transitions between the first 21 rotational levels of CS and for temperature between 20 K and 300 K still at the coupled states level. The comparison between the rate coefficients obtained by Lique et al. for CS-He and the results of Turner et al. for CS-para-H 2 (j = 0) show remarkable differences between the two sets of rate coefficients, which cannot be accounted only by the difference of mass of the colliders. A new study of the H 2 -CS collisions based on an accurate PES appeared to be necessary a) E-mail: t.stoecklin@ism.u-bordeaux1.fr to make a more realistic comparison between the two systems. A new four-dimensional intermolecular potential energy surface for the CS-H 2 collision was then developed and presented in a recent publication. 7 This surface based on a large grid of ab initio energies calculated at the coupled-cluster level was used to compute the bound state energies of the H 2 -CS complex. In the present work we present the cross sections and rate coefficients for the rotational transitions of CS in collisions with para-and ortho-H 2 for the 15 first rotational levels of CS. We compare them with those available for this system and with the recent calculations of Lique et al. for He-CS. 3 This paper is organized as follows. In Sec. II, a brief account of the parameters of the close coupling calculations is given, while the results are discussed in Sec. III.

II. METHOD

A. Potential energy surface

The surface used in the present work was recently described in details by Denis-Alpizar et al. 7 This is a 4D rigid rotor PES, computed for the H 2 distance being fixed at its vibrationally averaged value in the rovibrational ground state (r 0 = 1.448736 a 0 ), while the equilibrium bond length is used for CS (r e = 2.900619 a 0 ). 8 The equilibrium structure of the H 2 -CS complex is found to be linear with the carbon pointing toward H 2 at the intermolecular separation of 8.6 a 0 and the well depth is 173 cm -1 . The corresponding energy curve and those associated with secondary minima or saddle points are presented in Figure 1. The level of accuracy of this surface should be similar to the one obtained for similar systems like H 2 -CO. 9 Recent progresses in experimental techniques 10 conditions. The agreement obtained between theory and experiment by Chefdeville et al. 10 for the H 2 -CO inelastic collision cross sections in this domain of energy was quite good and we expect to reach the same level of accuracy for H 2 -CS.

B. Scattering calculations

The close coupling calculations were performed using our code DIDIMAT which was first applied to the H 2 -HF collisions 11 and also used recently to calculate the bound states of the H 2 -CS van der Waals complex. 7 The T matrix elements are obtained and used to calculate the cross sections for the transitions from an initial set of rotational levels of H 2 and CS labeled by the quantum numbers (j H2 , j CS ) to a final set labeled by the quantum numbers (j H2 , j CS ):

σ jH 2 ,jCS→j H 2 ,j CS = π (2j H2 + 1)(2j CS + 1)k 2 jH 2 ,jCS × ∞ J =0 (2J + 1) jH 2 +jCS j12=|jH 2 -jCS| j H 2 +j CS j 12 =|j H 2 -j CS | × |J +j12| l=|J -j12| |J +j 12 | l =|J -j 12 | T J jH 2 ,jCS,j12;j H 2 ,j CS ,j 12 2 ,
where J, j 12 , l, j 12 , and l are the quantum numbers associated with the total angular momentum J and with the operators j 12 = j H2 + j CS ; l = Jj 12 ; j 12 = j H2 + j CS ; l = Jj 12 . The magnitude of the wave vector is defined by

k 2 jH 2 , jCS = 2μ ¯2 [E -ε jH 2 , jCS
], ε jH 2 ,jCS being the eigenenergy associated with a given set of initial rotational state (j H2 , j CS ), E is the total energy, and μ the relative mass of the system. In Sec. III, we also report the rotational quenching cross sections of CS in collisions with H 2 which are defined as follows:

σ jH 2 ,jCS = j CS <jCS σ jH 2 ,jCS →jH 2 ,j CS . TABLE I.
Comparison between the inelastic cross section of CS colliding with para-H 2 and ortho-H 2 for different rotational transitions of CS. In each column the rotational basis set used for H 2 is reported. The rotational basis set used for CS included j CS values ranging from 0 to 18, while the kinetic energy was set to 10 cm -1 .

Transition

Para-H 2 Ortho-H 2 (j cs → j cs ) From the inelastic cross section, the corresponding thermal rate coefficients at temperature T are readily obtained from

j H 2 = 0 j H 2 = 0, 2 j H 2 = 1 j H 2 = 1, 3
k jH 2 ,jCS →j H 2 ,j CS (T ) = 8 πμ (k B T ) -3 2 ∞ 0 dE C σ jH 2 ,jCS →j H 2 ,j CS (E C )e - E C k B T ,
where E c is the collision energy of the channel labeled by j CS and j H2 . Full close coupling calculations were performed using the rigid rotor approximation from 10 -2 up to 1600 cm -1 with the rotational constants B CS = 0.8200462 cm -1 and B H2 = 60.853 cm -1 . Transitions among all levels up to j cs = 15 were computed for collisions involving ortho-and para-H 2 .

The highest rotational level used in any calculation was j CS = 30. For H 2 , the convergence of the rotational basis set was investigated. The comparison of some state selected crosssections computed for a selected total energy and a selected number of CS channels is presented in Table I. For the kinetic energy considered, adding a second rotational function of H 2 appears to change the cross section by 4% at the most in comparison with the calculation with only one rotational function. Thus, in order to reduce computer time, only one rotational level of H 2 was included in the basis set for total energies lower than the opening of the first excited para or ortho states of H 2 . The propagation was carried out to a maximum distance of 100 a o for the lowest energy, and convergence was checked as a function of the propagator step size.

III. RESULTS AND DISCUSSION

Figure 2 illustrates the energy dependence of the collisional de-excitation cross sections of CS(j CS = 5) in collision with para-H 2 . Apart from the usual many resonances associated with the H 2 -CS van der Waals well, a propensity to favor odd j CS over even j CS can be noticed at low collision energy while even j CS become favored over odd j CS at higher collision energy. This tendency was analyzed long ago by Mc-Curdy and Miller 12 in terms of an interference effect related to the odd or even anisotropy of the PES. The same tendency can be observed in Fig. 3 for collisions involving ortho-H 2 . However, the inversion in this case takes place at higher collision energy.

The dependence of the propensity rule on the collision energy and on the initial state of H 2 can also be understood by considering both the matrix elements of the potential (see, for example, Green 13 ) in a space fixed angular basis set and the expansion coefficients V m l 1 ,l 2 (R) of the potential in the Green angular basis set S m l 1 ,l 2 :

S m l 1 ,l 2 = P m l1 (cos θ 1 ) P m l2 (cos θ 2 ) cos(mϕ).

In this expansion, θ 1 and θ 2 are the polar angles of the H 2 and CS molecules, respectively, and ϕ is the relative azimuthal angle, 7 while P m l (cos θ ) are normalized associated Legendre polynomials. Our model PES 7 uses even values of l 1 ranging from 0 to 6 for the H 2 molecule, while the range is 0 ≤ l 2 ≤ 12 for CS and 0 ≤ m ≤ 2. At the total energies allowing only the H 2 state j H2 = 0 to be populated, only the l 1 = 0 components of the potential expansion give non-zero contributions to the potential matrix elements. At larger total energies, when the channel j H2 = 2 opens, the contributions due to the l 1 = 0, 2, 4 components give non-zero matrix elements and lead to the change of propensity rule appearing in Fig. 2. On the contrary, in the ortho case, the l 1 = 0 and 2 components already contribute when only the channel j H2 = 1 is open. At higher energy, l 1 = 0, 2, 4, 6 contribute when the channel j H2 = 3 opens. Another element which explains the differences of behavior observed between the ortho and para case is the fact that the j H2 = 2 channel is 365 cm -1 above the j H2 = 0 channel, while the j H2 = 3 channel is 600 cm -1 above the j H2 = 1 channel. This is illustrated in Fig. 4 which displays the cross sections associated with the opening of the j H2 = 2 channel in the para case and the j H2 = 3 channel in the ortho case. This difference of energy spacing is the main reason explaining why the cross sections of rotational excitation of H 2 in the para case are larger than those of the ortho case. The rotational quenching cross section of CS in collision with para-H 2 (j H2 = 0) and ortho-H 2 (j H2 = 1) for initial values of j CS = 1, 2, 3, 4, and 5 are presented, respectively, in Figs. 5 and6. While the very low collision energy interval [0.01; 1] cm -1 is not very useful to calculate rate coefficients at temperatures typical of the interstellar medium (ISM) it is however represented in these figures as it shows a resonance in this domain of energy for the rotational quenching cross 4 He (green lines) and with para-H 2 (j H 2 = 0) (red lines) as a function of the collision energy. The rotational transition j CS → j CS is reported in each panel. section of CS(j = 1). The different number of partial waves contributing to the intermolecular potential matrix elements in the para and ortho cases which we discussed above explains the larger cross sections obtained in the ortho case. This analysis may be supplemented by taking advantage of the H exchange symmetry of the H 2 molecule and averaging our PES over the ortho and para rotational wave functions of H 2 . We plotted the corresponding averaged potentials in Figure 7(a) for the para-H 2 case and Figures 7(b) and 7(c) for the ortho-H 2 case. The ortho-H 2 averaged potentials are deeper than the ones obtained for para-H 2 . This result explains the larger number of resonances obtained in the ortho case. In Fig. 8, we compare the cross sections obtained for the collision of CS with para-H 2 (j H2 = 0) with those calculated by Lique et al. for the 4 He-CS system for selected rotational transitions of CS. As seen in this figure, the largest differences between the two types of cross sections are observed at low energies, while at energies above 10 cm -1 the spacing between the two curves is more or less constant. The value of the ratio between the para-H 2 -CS and the He-CS cross sections calculated for the transitions reported in Fig. 8 and averaged over a set of collision energies (1, 50, 100, 500, 900, 1100 cm -1 ) is found to be equal to 1.34, which is very close to the square root of the ratio of the masses of the targets (i.e., 1.4). This approximation, which is frequently used by astronomers, is then valid for this collision for the energies above 1 cm -1 .

σ (a 0 2 ) E (cm -1 ) (1,5) → (1,0) (1,5) → (1,1) (1,5) → (1,2) (1,5) → (1,3) (1,5) → (1,
σ (a 0 2 ) E (cm -1 ) (0,5) → (2,4) (0,5) → (2,7) (0,4) → (2,3) (0,3) → (2,2) (0,2) → (2,4) (1,5) → (3,4) (1,5) → (3,7) (1,4) → (3,3) (1,3) → (3,2) (1,2) → (3,
The rate coefficients at 100 K for the de-excitation transitions of CS(j CS = 15) in collisions with para-H 2 and ortho-H 2 are presented in Fig. 9, along with the scaled He-CS results of Lique et al. available in the BASECOL2012 database. 14 The scaled He-CS rate coefficients are smaller than the H 2 -CS coefficients for small j CS , and larger for large j CS . The largest error factor is about 3. However, the agreement is good for the small odd j CS for para-H 2 -CS.

In the same figure we also present the rate coefficient calculated by Turner et al. taken from BASECOL2012. These rate coefficients are too small for all j CS by a factor ranging from 1.2 to 2.5. This discrepancy is a consequence of the approximations made by Turner et al. They used the crude PES of Green and Chapman adapted from a Gordon and Kim electron gas model for CS-He calculations and the coupled states approximation for the scattering calculations. From a comparison between the electron gas model and a realistic PES for the CO-He system, Turner et al. reported that rotational excitation rates are in error by a factor of 2 or 3, which is similar to the error factor we have obtained for the CS-H 2 system.

IV. CONCLUSION

A new set of rigid rotor close coupling cross sections for the inelastic collisions of CS with H 2 was computed using a new 4D potential energy surface. For the collision of CS with para-H 2 , within the energy range from 0.1 to 1000 cm -1 , we observe a propensity to favor odd j CS over even j CS at low collision energy while at energy close to 1000 cm -1 , the propensity rule is reversed. The same tendency is observed for collisions involving ortho-H 2 . However, the inversion takes place at higher collision energy. This effect is related to the opening of the second channel of H 2 (j H2 = 2 for para-H 2 or j H2 = 3 for ortho-H 2 ).

The cross sections for the rotational transition of CS in collision with para-H 2 were compared with the scaled crosssections available for the collision of CS with He. The usual square root of the relative mass of the colliders, which is used by astronomers to obtain the cross sections between a molecule and H 2 from the data available with 4 He, is found to be a good qualitative approximation for this system.

The rate coefficients calculated with the crude electron gas model show more discrepancy with our accurate results. Nevertheless, the variation of the rate as a function of j CS is correctly reproduced. Scaling up this rate by a factor of 2 gives a good qualitative agreement. The ab-initio points were computed by the group of François Lique (Université du Havre) at the CCSD(T)-F12 level, using the aug-cc-pVTZ basis set. Furthermore, they analysed the long range behaviour of the grid of points computed with the latter ab-initio method.

We performed the fit of the ab-initio points to an analytical function using the least-square procedure. The angular part of the function is represented by a product of the normalized associated Legendre functions and a cosine function:

ȳm l 1 ,l 2 = P m l 1 (cosθ ) P m l 2 (cosθ) cos(mϕ) (4.1)
where the angles θ and θ correspond to the relative orientation of the HCN and the H 2 molecules, and ϕ is the azimuthal angle. The first radial functions, called ν m l 1 ,l 2 are represented in Fig 4 .1. The ν 0 0,0 is the isotropic term. The equilibrium structure of the system is linear HCN-H 2 .

The PES of the HCN-H 2 system have then been used for the calculation of the bound states of the HCN-para-H 2 and HCN-ortho-H 2 complexes separately. In these calculations, both molecules were treated as linear rigid rotors. Also, we averaged the PES over the rotational wave function of the H 2 molecule to explain the remarkable difference between the dissociation energies of the HCN-para-H 2 and HCN-ortho-H 2 species.

Finally, we used the bound levels to compare the calculated frequencies transition with the experimental values reported in the literature. The agreement found was quite good, confirming the good quality of the PES.

Publication

THE JOURNAL OF CHEMICAL PHYSICS 139, 224301 (2013)

A new ab initio potential energy surface for the collisional excitation of HCN by para-and ortho-H 2

Otoniel Denis-Alpizar, 1,2,a) Yulia Kalugina, 3,4 Thierry Stoecklin, 1 Mario Hernández Vera, 3,5 and François Lique We present a new four-dimensional potential energy surface for the collisional excitation of HCN by H 2 . Ab initio calculations of the HCN-H 2 van der Waals complex, considering both molecules as rigid rotors, were carried out at the explicitly correlated coupled cluster with single, double, and perturbative triple excitations [CCSD(T)-F12a] level of theory using an augmented correlationconsistent triple zeta (aVTZ) basis set. The equilibrium structure is linear HCN-H 2 with the nitrogen pointing towards H 2 at an intermolecular separation of 7.20 a 0 . The corresponding well depth is -195.20 cm -1 . A secondary minimum of -183.59 cm -1 was found for a T-shape configuration with the H of HCN pointing to the center of mass of H 2 . We also determine the rovibrational energy levels of the HCN-para-H 2 and HCN-ortho-H 2 complexes. The calculated dissociation energies for the para and ortho complexes are 37.79 cm -1 and 60.26 cm -1 , respectively. The calculated ro-vibrational transitions in the HCN-H 2 complex are found to agree by more than 0.5% with the available experimental data, confirming the accuracy of the potential energy surface. © 2013 AIP Publishing LLC.

[http://dx.doi.org/10.1063/1.4833676]

I. INTRODUCTION

Hydrogen cyanide (HCN) is a central molecule in the physical-chemistry of many media. For example, HCN is an important intermediate in the combustion reactions of hydrocarbon flames containing a nitrogen source. 1 HCN is also an important constituent of earth and planetary atmospheres 2 and plays a major role in the physical chemistry of these atmospheres. For example, in Titan atmosphere, HCN rotational emission is the dominant cooling process and is therefore responsible for the thermal structure of the thermosphere. 3 From the astrophysical point of view, HCN is ubiquitous. Hydrogen cyanide is one of the most observed molecule in the interstellar medium (ISM). Owing to its large dipole moment, this molecule and its closely related molecule HNC are frequently used by astronomers to determine the physical and chemical conditions in many regions of the ISM. [4][5][6][7][8][9] They are also powerful probes of high density gas, as at the opposite of molecules like CO or CS, they do not seem to deplete on grain surfaces in the denser cold part of prestellar cores. 10 The analysis of the astronomical HCN emission spectra requires the knowledge of accurate spectroscopical data, as well as collisional rate coefficients with the most abundant interstellar species. Indeed, collisional processes contribute, in competition with the radiative processes, to populate the rovibrational levels of the molecules. Consequently, the calcua) otonieldenisalpizar@gmail.com b) francois.lique@univ-lehavre.fr lation of rate coefficients for the rotational excitation of HCN by H 2 (H 2 being the most abundant collisional partner is the ISM) has been a major goal for astrochemistry.

The first work dedicated to the rotational excitation of HCN by collisions with He atoms (as a model for H 2 ) was performed by Green and Thaddeus 11 in 1974. Then, many collisional studies of astrophysical interest have been devoted to the HCN-He collisional system and we just mention here the most recent ones. [12][13][14] However, in all these studies, He was used as a model for H 2 and recent results on rotational excitation of SO, 15 SiS, 16,17 or HNC 18 have pointed out that rate coefficients for collisions with H 2 are generally different from those for collisions with He, particularly for the case of collision with ortho-H 2 .

Recently, Ben Abdallah et al. 19 have studied the rotational excitation of HCN by collisions with H 2 and provided rate coefficients for the rotational and hyperfine (de-)excitation of the HCN by para-H 2 at low temperature. However, these authors computed and used a potential energy surface (PES) averaged over the H 2 rotation, resulting in simplified dynamical calculations only valid for the collisional excitation by H 2 (j = 0) at low temperatures. Then, the validity of the approximation used by Ben Abdallah et al. 19 needs to be checked by comparing with calculations involving a global PES. Moreover, as HCN is also observed in high temperatures astrophysical environments, it seems also important to compute a new global HCN-H 2 PES in order to provide new rate coefficients for the rotational excitation of HCN by H 2 (para-and ortho-H 2 ) at high temperatures. These new data will provide the astronomical community with the necessary tools to interpret HCN emission that will be observed with high spatial and spectral resolutions from the Atacama Large Millimeter/sub-millimeter Array (ALMA) interferometer.

From the experimental point of view, to the best of our knowledge, collisional excitation studies of HCN by H 2 are not available. However, several studies dedicated to the measurement of the HCN-H 2 van der Waals spectrum were performed during the last decade. [20][21][22][23] High resolution infrared 21,22 and millimeter-wave spectroscopy combined with a pulsed supersonic jet technique 20 have been used to study the van der Waals complex structure in great detail. Very recently, Fourier-transform microwave spectroscopy was applied to observe rotational transitions in the HCN-H 2 complex (both para-and ortho-H 2 species) in order to improve hyperfine molecular constants of the complex. 23 From these studies, Ishiguro et al. 20 and Moore et al. 21 were able to show that ortho-H 2 and para-H 2 bind at opposite ends of HCN. Moore et al. 22 used an average one-dimensional ab initio model of the PES to justify this analysis.

In this paper, we present a new four-dimensional (4D) potential energy surface for the HCN-H 2 system obtained from first principle calculations. As a first application and in order to check the accuracy of the new PES, we compute the bound energy levels of HCN-para-H 2 and HCN-ortho-H 2 complexes and compare the calculated HCN-H 2 frequencies with the experimental ones. 20,23 This paper is organized as follows. The ab initio calculations and the analytical representation of the PES are described in Sec. II. Features of the PES are discussed in Sec. III. The bound levels of the complexes and the comparison of calculated and experimental transitions frequencies are presented in Sec. IV.

II. Ab initio CALCULATIONS AND ANALYTICAL REPRESENTATION

A. Ab initio calculations

In the present work, we use the Jacobi coordinate system presented in Fig. 1. The origin of coordinates is placed at the center of mass of the HCN molecule and the vector R connecting HCN and the center of the H 2 bond is directed along the z-axis. The rotation of HCN and H 2 molecules is defined by θ and θ angles, respectively. ϕ is the dihedral angle.

A very recent study dedicated to the ro-vibrational excitation of HCN by He 14 showed that neglecting the bending of HCN is a very reasonable approximation for treating the FIG. 1. Jacobi coordinate system of the HCN-H 2 system. pure rotational rotational excitation at low and moderate temperatures. Indeed, some of us computed a 4D PES taking in account the bending angle of the HCN molecule and found that the results for the pure rotational excitation were in good agreement with those considering the collisional partner as rigid. Consequently and as we are mainly interested in the pure rotational excitation of the title system, the collision partners are then considered as rigid (i.e., we neglect effects of the normal mode vibrations of the HCN and H 2 ).

As shown before, 24 a better description of the intermolecular potential is obtained by fixing the molecular distance at its averaged value in the ground vibrational level rather than at the equilibrium distance. Accordingly, we used a H 2 bond distance r H-H = 1.448736 a 0 . 25 For HCN, we used the following interatomic distances r CN = 2.17923 a 0 , r CH = 2.01350 a 0 (corresponding to the linear equilibrium geometry of HCN 26 ).

Ab initio calculations of the PES of the HCN-H 2 van der Waals complex being in its ground electronic state were carried out at the explicitly correlated coupled cluster with single, double, and perturbative triple excitations [CCSD(T)-F12a] 27 level of theory using an augmented correlation-consistent triple zeta (aVTZ) basis set. 28 The calculations were performed using Molpro 2010 package. 29 Additionally to the aVTZ basis set, density fitting and resolution of the identity techniques (needed for F12 calculations) require auxiliary basis sets. In the present study, we used the default corresponding auxiliary basis sets 30 from Molpro basis set library.

In all calculations, the interaction potential V (R, θ, θ , ϕ) was corrected by the basis set superposition error (BSSE) with the Boys and Bernardi counterpoise scheme 31

V (R, θ, θ , ϕ) = E HCN-H 2 (R, θ, θ , ϕ) -E HCN (R, θ, θ , ϕ) -E H2 (R, θ, θ , ϕ) , (1)
where the energies of HCN and H 2 are computed in the full basis set of the complex.

In order to evaluate the accuracy of CCSD(T)-F12a method in comparison with the standard CCSD(T) method using the aug-cc-pVXZ (X = D, T, Q) basis sets (hereafter aVDZ, aVTZ, aVQZ) and including results extrapolated to the Complete Basis Set (CBS) limit, 32 we have carried out additional calculations for selected angular orientations of the HCN-H 2 complex. The results are reported in Fig. 2.

As can be seen, the CCSD(T)-F12a together with diffuse basis set of aVTZ quality fully reaches the CCSD(T)/CBS accuracy as previously found for other van der Waals systems such as HCl-He 33 or O 2 -He. 34 This comparison shows again that the CCSD(T)-F12a leads to an excellent description of the interaction energies using a relatively limited atomic basis set.

In order to accurately describe the anisotropy of the potential energy surface, the calculations were carried out for a large grid of angular orientations: we vary the θ angle from 0 • to 180 • with a step of 10 • , the θ angle from 0 • to 180 • with a step of 15 • , and the ϕ angle from 0 • to 90 • with a step of 15 • . R-distances vary from 3.5 to 30 bohrs in order to get 39 radial grid points for each angular orientation.

B. Analytical representation

The angular part of the analytical representation of the PES is represented by a product of normalized associated Legendre functions and a cosine function. Hence, the potential is expressed such as

V (R, θ , θ, ϕ) = 4 l1=0 17 l2=0 min(l1,l2,3) m=0 ν m l1,l2 (R) ×P m l1 (cosθ ) P m l2 (cosθ ) cos(mϕ), ( 2
)
where ν m l1,l2 (R) are the radial functions. Only even values of l 1 are included in this expansion because of the H exchange symmetry of H 2 . For each value of the radial grid, we fitted the ab initio points using a least squares method. We used a weight function for energies higher than 5000 cm -1 as

w(E(θ , θ, ϕ)) = V 0 E(θ , θ, ϕ) . ( 3
)
By trial and error, the values of the weight function parameter were found to be V 0 = 5000 cm -1 . In total, we get 135 angular coefficients for each radial point.

The analytical representation of the radial coefficients ν m l1,l2 (R) was divided in three parts, corresponding to the shortrange, the intermediate, and the long-range regions as detailed in the three following expressions:

F sr (R) = 10 n=0 e -αR R n C m,n l1,l2 R ≤ 6a 0 , ( 4
)
F ir (R) = 11 k=4 C m,k l1,l2 R k 6 < R < 11a 0 , ( 5
)
F lr (R) = 8 k=4 t k (βR) R k C m,k l1,l2 R ≥ 11a 0 , ( 6
)
where the t k (x) is the Tang-Toennies damping function

t k (x) = 1 -e -x k i=0
x i i! (7) and α = 1.8 a -1 0 , β = 4.0 a -1 0 . These analytical functions allowed us to get an accurate representation of the PES using a least-square procedure to compute the final coefficients. The quality of the fitted surface was checked by evaluating the root-mean-square (RMS) of the differences between the ab initio energies and the fitted energies. For negative energies, the RMS was on the order of 10 -4 cm -1 . For energies in the range 0 ≤ E ≤ 1000 cm -1 , the RMS was 0.008 cm -1 while for 1000 < E ≤ 10 000 cm -1 we found the RMS = 3.656 cm -1 . The analytical representation of the PES is available from the authors upon request.

III. FEATURES OF THE POTENTIAL ENERGY SURFACE

The global minimum (D e = -195.20 cm -1 ) of the 4D PES corresponds to a linear geometry of the HCN-H 2 complex with θ = 0 • , θ = 0 • or 180 • , and R = 7.20 a 0 . A secondary minimum of -183.59 cm -1 is found for the T-shape configuration associated with θ = 90 • , θ = 180 • , and R = 7.75 a 0 . This finding is in good agreement with the experimental one of Moore et al. 22 which concludes that linear HCN-H 2 is the most stable geometry of the complex.

The contour plots in Figs. 3 and4 show the anisotropy of the interaction with respect to the HCN and H 2 rotation. In Fig. 3, we present the anisotropy of the potential as a function of θ for fixed values of θ and ϕ. As one can see, there is relatively strong anisotropy of the PES with respect to the θ Jacobi angle. Figure 4 shows the interaction energies for ϕ = 0 • with R relaxing in the range [7.0, 8.0] a 0 . Here again, we found a relatively strong anisotropy of the PES with respect to these two Jacobi angles. Therefore, we can anticipate that rotational state of H 2 will significantly influence the magnitude of the HCN excitation cross sections.

In Fig. 4, we can see the secondary minimum associated with the T-shape configuration (with the H of HCN pointing to the center of mass of H 2 ) which is only 11.61 cm -1 above the global minimum. Then, we expect this configuration of the HCN-H 2 complex to be almost equally important for the dynamics of the system.

The analysis of the PES features may be supplemented by taking advantage of the H exchange symmetry of the H 2 molecule and averaging our PES over the para and ortho rotational wave function of H 2 such as

V av = Y * j1,k (θ , ϕ)|V (R, θ , θ, ϕ)|Y j1,k (θ , ϕ) , (8)
where j 1 denotes the rotational angular momentum of H 2 . The projection of j 1 on the intermolecular axis is k. This integral was evaluated numerically, using a Gauss-Chebyshev quadrature of 20 points for the integration over ϕ and a Gauss-Legendre quadrature with 10 points over θ . of the complex is obtained for the linear configuration of HCN-H 2 and for an intermolecular separation of 8.14 a 0 . This surface can be compared with the approximate PES of Ben Abdallah et al. 19 Our averaged surface, with a minimum of -79.23 cm -1 , is 8.46 cm -1 deeper than the one presented by Abdallah et al. 19 while the secondary minimums obtained in both cases for a value of θ close to 60 • , differ only by 1.1 cm -1 . But the most interesting result is that they found the most stable configuration for θ = 164 • while we obtained it for a linear configuration of HCN in agreement with the conclusions of the experiments of Ishiguro et al. 20 These differences result from the fact that they averaged over ϕ and θ using only three different angular configurations of the complex.

In the ortho-H 2 (j 1 = 1) case, the wave function can be (j 1 = 1, k = 0) or (j 1 = 1, k = ±1). We plotted the two corresponding averaged potentials in Figs. 5(b) and 5(c). These averaged potential are deeper than the one obtained for para-H 2 . In both cases, the minimum corresponds to a linear configuration of the complex, but with θ = 0 • or 180 • . These two configurations of the HCN-ortho-H 2 complex are corresponding to the H atom and to the N atom of HCN pointing towards H 2 , respectively, while experimental observations of Moore et al. 21,22 lead to a linear configuration of the complex with the H atom of HCN pointing towards H 2 .

Then, in order to understand the nature of interactions in HCN-H 2 system, we have performed analytical calculations of interaction energy based on multipolar expansion. In the framework of the long-range approximation, 35 the potential energy of two interacting systems can be written as

V = E elec + E ind + E disp , ( 9
)
where E elec , E ind , and E disp are the electrostatic, induction, and dispersion contributions to the total interaction energy of the complex. For interacting HCN and H 2 molecules, the electrostatic, induction, and dispersion terms from R -4 through the order of R -8 are a Properties were calculated with respect to the origin at the center of mass with N atom along the positive direction of z axis. b Theoretical values by Maroulis and Pouchan. 40 where α = (α xx + α yy + α zz )/3 is a mean polarizability. Superscripts A and B denote molecules HCN and H 2 , respectively. In the present work, we adopted the tensor notations of Buckingham. 35 Definition of multipole moments and polarizabilities is presented in Table I. Tensor T αβγ . . . ν = ∇ α ∇ β ∇ γ . . . ∇ ν R -1 is the symmetric tensor relative to the permutation of any pair of subscripts. The tensor T of rank n, T (n) , is proportional to R -(n+1) . There is a summation over repeated indexes. Expression (12) was obtained in "constant-ratio" approximation, 36 which allows to evaluate dispersion contribution through static properties of subsystems and isotropic C 0 6 dispersion coefficient. One should note, that all properties in Eqs. ( 10)- (12) are represented in the coordinate system of the complex.

E elec = - 1 3 T αβγ μ A α B βγ + 1 9 T αβγ δ A αβ B γ δ - 1 105 T αβγ δ μ A α B βγ δ - 1 45 T αβγ δ A αβγ B δ + 1 315 T αβγ δ ϕ A αβ B γ δ ϕ + B αβ A γ δ ϕ (10) - 1 945 T αβγ δ ϕν μ A α B βγ δ ϕν - 1 1575 T αβγ δ ϕν A αβγ B δ ϕν - 1 2835 T αβγ δ ϕν A αβγ δ B ϕν , E ind = - 1 2 T αβ T γ δ α B αβ μ A β μ A γ + 1 3 T αγ T βδφ α B αβ μ A γ A δφ - 1 18 T αγ δ T βφ α B αβ A γ δ A φ + α A αβ B γ δ B φ (11) - 1 15
T αγ T βδφ α B αβ μ A γ A δφ - 1 15 T αϕ T βγ δ E B α,βγ δ μ A ϕ μ A - 1 6 T αβϕ T γ δ C B αβ,γ δ μ A ϕ μ A , E disp = - C 0 6 6α A α B T αβ T γ δ α A αγ α B βδ - 2 3 T αβ T γ δφ α B αγ A A β,δφ + 1 3 T αβγ T δφ α A αδ C B βγ,φ (12)
+ α B αδ C A βγ,φ + 2 15 T αβ T γ δφ α A αγ E B β,δφ + α B αγ E A β,δφ ,
The values of multipole moments and polarizabilities used in analytical calculations are presented in Table I. These values were calculated in the coordinate system of monomers at the CCSD(T)/aV5Z level of theory using finitefield method of Cohen and Roothaan 37 since it was not possible to compute them at the CCSD(T)-F12 level. The isotropic dispersion coefficient C 0 6 = 34.148 E h a 6 0 was calculated at the CCSD PROPAGATOR method 38 using the Molpro routine. 39 In Fig. 6, we present the major contributions to the interaction energy in long-range approximation for the equilibrium configuration of HCN-H 2 complex. One can see that the electrostatic interactions are dominant for this complex. The leading electrostatic term is the H 2 quadrupole-HCN dipole interaction proportional to μR -4 ; the H 2 quadrupole-HCN quadrupole and H 2 quadrupole-HCN octupole interactions proportional to R -5 and R -6 , compensate each other. The other contributions are relatively small and do not contribute much to the behavior of the long range PES.

Figure 7 shows the curves of the interaction energy for different orientations of monomers at long-range separations. There is a good agreement between the CCSD(T)-F12a results and analytical calculations for R > 13 a 0 , which confirms FIG. 6. Different contributions to the interaction energy of HCN-H 2 system for equilibrium configuration with θ = 0 • and θ = 0 • . Energy is in cm -1 . Solid black line -CCSD(T)-F12a calculations; dashed black line -total interaction energy in long-range approximation; solid color lines -major contributions to interaction energy from Eqs. ( 10) to (12).

that both analytical and ab initio calculations have a correct physical behavior. This also means that the present PES can be used with confidence to describe cold molecular collisions.

IV. BOUND STATES CALCULATIONS

We used a close coupling approach to compute the rovibrational energy levels of the HCN-H 2 complex. As described in a previous study, 41,42 we implemented the bound state calculations in our diatom-diatom Close Coupling scattering code according to the approach proposed by Danby 43 and Hutson 44 for the R-matrix and log-derivative propagators, respectively. Recently, we used the same approach to calculate the bound states of several diatom-diatom complexes. 45,46 Briefly, our code solves the rovibrationally inelastic Close Coupling equations in the space fixed frame using the Jacobi coordinates. The vibrational levels of the van der Waals complex are calculated by performing calculations for all the FIG. 7. Comparison of ab initio and analytical calculations of the PES for fixed angular arrangements. Energy is in cm -1 . Solid lines -CCSD(T)-F12a calculations; dashed lines -analytical calculations. values of the total angular momentum J and parity leading to bound states.

In the present study, we performed the diatom-linear molecule calculations using for both HCN and H 2 , a rigid rotor description. The calculations were done separately for ortho-and para-H 2 . We used the vibrationally averaged rotational constants for both molecules (B = 59.322 cm -1 for H 2 47 and B = 1.47822 cm -1 for the HCN 48 ).

The convergence of the bound levels of HCN-H 2 complex with respect to the HCN and H 2 rotational basis was tested. Twelve rotational states were included in the basis set describing the HCN molecule while for para-H 2 we included two rotational states. In the case of HCN-ortho-H 2 , we also included two rotational states for ortho-H 2 and 13 for HCN. The maximum propagation distance was 50 a 0 .

The calculations were performed for two values of the propagator step size (0.1 and 0.05 a 0 ) and the values of the bound state energies were obtained from a Richardson extrapolation as suggested by Hutson. 44 The results of the lower bound states calculations are presented in Tables II andIII. The other energy levels are available in the supplementary material. 49 The rovibrational energies are given relative to the ground state energy of infinitely separated HCN and para-H 2 . The spacing between the levels j 1 = 0 and j 1 = 2 in the para form and between the j 1 = 1 and j 1 = 3 levels in the ortho form are larger than the well depth of the potential. Consequently, all the HCN-para-H 2 bound states are associated with j 1 = 0 while those of HCN-ortho-H 2 involve only j 1 = 1. For ortho-H 2 , the quantum number associated with the angular momentum j 12 = j 1 + j 2 , where j 2 designates the HCN rotational state and is also necessary to classify the levels as j 1 is non-zero.

In Table II, we report, for the HCN-para-H 2 complex, the energy, the total angular momentum J, and the parity ε (where ε = (-1) j1+j2+L ) of the lower bound levels (E < 20 cm -1 ). We also assign the approximate quantum number j 2 and the orbital quantum number L. In Table III, we report the same TABLE II. Lower bound levels (E < 20 cm -1 ) of the HCN-para-H 2 van der Waals complex. The approximate quantum numbers j 2 and L are also given. information for the lower bound levels (E < 77 cm -1 ) of the the HCN-ortho-H 2 complex. In addition, we also assign the approximate quantum number j 12 . This attribution allowed us to determine the transitions of the system in the ground state of HCN and to compare them with the experimental results.

Energy Energy State State J ε (j 2 , L) ( c m -1 ) J ε (j 2 , L) ( c m -1 ) 0 + (0,0) -37.79 2
The maximum value of the total angular momentum J leading to bound states is 9 for HCN-para-H 2 and 11 for HCN-ortho-H 2 . The total number of bound states supported by our PES is 101 for HCN-para-H 2 and 330 for HCN-ortho-H 2 . There are around three times more states for the ortho form than for the para form. This is related to the fact that with j 1 = 0, a single value of j 12 is possible for any given value of j 2 , while with j 1 = 1, three values of j 12 are obtained. A similar ratio between the number of ortho and para states was observed previously for the H 2 -O 2 complex. 45 The dissociation energies of the ground states of the para and ortho species are 37.79 cm -1 and 60.26 cm -1 , respectively. These values are at least more than 4 times the dissociation energy found for the HCN-He complex (8.986 cm -1 ) 14 confirming that van der Waals complexes are generally more bonded with H 2 than with He. A comparison with the result of the bound levels calculations for the isoelectronic CO-H 2 complex, 47 also shows a smaller dissociation energy for the para-H 2 complex than for the ortho-H 2 complex. However, the HCN-H 2 complexes have a larger dissociation energy than the complex CO-H 2 (23.7 cm -1 for CO-para-H 2 and 30.8 cm -1 for CO-ortho-H 2 ), as we expected given the difference between the wells depths of these two surfaces.

The bound state energies were used to determine the transitions frequencies in the HCN-H 2 van der Waals complex. In the frequencies calculations, we considered HCN in its ground rotational state. These transitions are compared in Table IV with the spectroscopic data available published by Ishiguro et al. 20,23 These authors reported several lines including the splitting into several hyperfine components due to the spin angular momentum of the nitrogen nucleus (I = 1). As our calculations do not include the hyperfine structure and because the spin splitting is very small in comparison with the spacings of the rotational lines, we compare our results with those of Ishiguro et al. 20,23 averaged over the hyperfine components. The agreement between our results and experimental ones is better than 0.5% in all cases. It confirms the accuracy of the new HCN-H 2 PES.

V. CONCLUSION

We have developed a new 4D analytical PES for the HCN-H 2 van der Waals complex based on a large grid of ab initio points obtained at CCSD(T)-F12a level and using an aVTZ basis set. The equilibrium structure of the HCN-H 2 complex was found to be linear with the nitrogen pointing towards H 2 . The corresponding well depth is 195.20 cm -1 . We found a secondary minimum only 11.61 cm -1 above the global minimum in which the H atom of HCN is pointing towards the center of mass of H 2 .

As a first application, the rovibrational bound states were computed within the rigid-rotor approximation. The total number of bound states supported by our PES is 101 for HCN-para-H 2 and 330 for HCN-ortho-H 2 . The dissociation energies of the ground states of the para and ortho species are 37.79 cm -1 and 60.26 cm -1 , respectively. The calculated transitions frequencies are found to be in very good agreement with the experimental available data. This level of agreement suggests that our PES is accurate enough for computing accurate inelastic cross sections.

Finally, the present PES was also compared with the previous existing models. The differences which were found and taking into account the astrophysical importance of this molecule give a supplementary interest to new Close Coupling calculations of the inelastic cross sections which will be performed using this new surface. In particular, inelastic cross sections and rate coefficients will be provided for the first time for the ortho-H 2 collisional partner. These results will be very valuable for the astrophysical modeling of "hot environment" where ortho-H 2 is abundant. in the ISM. For these reasons, we decided to study the collision of He with two molecules often observed in the ISM, namely HCN and C 3 , and to include the coupling of the lowest frequency vibrational mode with the intermolecular motion. When this work was started there was few examples of such study. Let's mention the pioneering works of Clary [24, 25], using the IOS approximation, and later, Lan et al [86]. Recently, the group of Dagdigian and Alexander [31, 32] computed the bending wave function of the non-linear triatomic system CH 2 , averaged the interaction potential over the bending motion and did the calculation of the collision with He. In our first paper dedicated to the He-HCN collision, we developed a similar approach which we called rigid bender average approximation (RBAA). This method is detailed below. In the rigid bender approximation, a triatomic molecule is a 1D vibrational system: the bond lengths are fixed while the bending motion is allowed. The main differences between the previous two studies is that the equilibrium geometry of CH 2 is a symmetric top and its bending frequency is 963 cm -1 [34] while HCN is linear at equilibrium and its bending frequency is 712 cm -1 [33]. The coupling between rotation and bending in a linear molecule leads to l-doubling. Indeed, l-type transitions are observed in the ISM for the HCN molecule [35]. It appeared then necessary to check the reliability of the RBAA method for describing these transitions between almost degenerate levels by comparing with the results of a more accurate method taking into account exactly the coupling between the collisional motion and the bending motion. We developed such a method using the rigid bender approximation within the close coupling framework (RB-CC). This new method is presented below.

Rigid bender HCN molecule

The first step in the theoretical treatment of the collisions of He with rigid bender HCN is to determine the bending energies and wave functions of the triatomic molecule. The use of the rigid bender approximation requires using internal coordinates which were very popular in the eighties. Different Hamiltonians in internal coordinates were proposed at that time for different kind of triatomic systems [87, 88]. All these Hamiltonians have no singularities for linear geometries and can then be used to treat linear molecules. In our first work dedicated to HCN we neglected the coupling between bending and rotation and calculated separately the bending ener- gies and wave functions of HCN for j=0. For non zero values of j we simply added the linear rotational energies and used the same bending functions. Within this approach, for j=0, we used a rigid bender restriction of the Hamiltonian developed by Carter and Handy [87] and written in terms of the two bond lengths and the bending angle. For the HCN molecule, the stretching frequencies are quite larger than the bending one. Thus the rigid bender approximation can be safely applied as long as we consider only vibrational transitions with an energy smaller than the first excited stretching level. For j=0, the total Hamiltonian of the triatomic molecule can be expressed as

H = H V + V (γ), (5.1)
where H V is the kinetic part of the Hamiltonian and V(γ) is the potential which is a function of the bending angle γ.

An analytical form of this potential is obtained from its expansion in Legendre polynomials, V (γ) =

3 n=0
C n P n (cos γ). This potential, shown in Fig. 5.1, was fitted from a grid of 22 ab-initio points computed at the CCSD(T) level using an aug-cc-pVQZ basis set.

The diagonalization of the rigid bender Hamiltonian for j = 0 in a Legendre polynomials basis set gives the energies ε n and the wave functions χ n (γ) which are then used to perform the rigid bender average calculations as detailed below.

In our second paper dedicated to HCN, the coupling between rotation and bending was treated exactly within the rigid bender approximation using the Hamiltonian developed by Sutcliffe and Tennyson [89, 90]. This Hamiltonian is written

H = H (1) V + H (2) V + H (1) V R + V (γ).
(5.

2)

The two first terms are two contributions to the kinetic energy which are independent of the rotation and H RV is the rovibrational term. The bending part of the wave function is expended in normalised associated Legendre polynomials χ ν = n C ν n P K n (cos γ) while the rotational part is expected to be well described by an asymmetric top wave function, K≥0)

A p k |j KM, p for any bent configuration, where j is the total angular momentum of HCN, and K and M are its projections along the BF z -axis and SF z-axis, respectively. K is the absolute value of K. In order to dispose of wave functions of a given parity, one uses the symmetrised basis set proposed by Sutcliffe and Tennyson [90] |jn KM, p = P K n (cos γ)|j KM + (-1) p P -K n (cos γ)|j -KM , (

where p takes the values 0 or 1 for K = 0, and can only be equal to zero for K = 0. As the action of the parity operator on a symmetric top wave function is Π|j KM = (-1) K |j -KM , one obtains, using the property P K n (cos γ) = (-1) K P -K n (cos γ), that the total parity of the state |jn KM, p is (-1) j+p . This basis set allows to derive the analytical matrix elements of the Hamiltonian. From the diagonalisation of the Hamiltonian, we obtain for each value of j and p, the rovibrational energies ε p νj and the corresponding wave functions which can be written as (in the second paper, there is a misprint in equation 17),

χ νp jM (γ) = K≥0 Γ νp j, K (γ) 2(1 + δ K0 ) |j KM + (-1) p+ K |j -KM , (5.4)
where Γ νp j,

K (γ) = n≥ K C νjp n K P K n (cos γ).
For K > 0, the coupling between rotation and bending results into the splitting of each level in two non degenerate components. This splitting is called l-doubling.

Rigid bender close coupling equations

The Schrödinger equation for the collision of HCN with He can be written

- 2 2µR ∂ 2 ∂R 2 + L 2 2µR 2 + H HCN + V int (R, γ, θ, ϕ) ψ = E tot ψ, (5.5)
where µ is the collisional reduced mass, L 2 is the relative angular momentum, H HCN is the Hamiltonian of the isolated HCN molecule, and V int is the interaction potential between the atom and the rigid bender molecule. This potential depends on the bending angle γ, the distance R from the center of mass of HCN to the He atom, the angle ϕ of rotation around the axis defined by the H atom and the center of mass of CN, and the angle θ between the latter axis and the axis defined by the He atom and the center of mass of HCN as can be seen in Fig. 5.2 We developed two different methods to solve this equation respectively presented in the two following joined publications. In the first approach called rigid bender averaged approximation (RBAA), the rotation is decoupled from the bending motion and the equation 5.5 is solved using the PES averaged over the bending motion:

V ave (R, θ) = χ n |V int (R, θ, ϕ = 0 • , γ)|χ m . (5.6)
The problem becomes formally equivalent to the one of an atom colliding with a rovibrating linear molecule, the vibration of the molecule being in fact its bending vibration. The averaging of the potential can be done either using the rigid bender wave functions of the j=0 Hamiltonian or those of the Hamiltonian 5.1. In the first case, the rovibrational energies used for the dynamics are the bending energies added of the rotational energy E νj = Bj(j + 1) + ε ν where B is the rotational constant of the rigid linear HCN. In the second case, we use the rovibrational energies ε p νj obtained from the diagonalisation of the Hamiltonian 5.2.

Within the second approach called the rigid bender close coupling approach (RB-CC), equation 5.5 is solved exactly within the rigid bender approximation. The angular basis set used to expand the wave functions of equation 5.5 in the space fixed coordinates can be written as

|νjlpJM = m j m l 2J+1 4π jm j lm l |JM Y m l l ( R) K≥0 Γ νp j, K (γ) √ 2(1+δ K0 ) × D j * m j , K (φ R HCN , θ R HCN , 0) + (-1) p+ K D j * m j ,- K (φ R HCN , θ R HCN , 0) (5.7)
In this expression, we coupled the HCN total angular quantum number j with the relative momentum l to get angular wave functions for given values of the total angular momentum J and of the parity = (-) j+l+p .

The resulting Close Coupling equations presented in the second joined paper as well as those of the RBAA approach were implemented in the Bordeaux scattering code Newmat. The use of these two methods to calculate bound states of the complex using the method presented in section 2.4 was also included.

In the first joined paper, we determined the bound levels of the complex in the rigid monomer approximation (RMA), in which HCN is treated as a rigid linear molecule and we compared the results with those obtained when using the RBAA method. The energies obtained at the RBAA level are found to be systematically above those obtained using the RMA. This can be understood by reminding that the minimum of the PES is associated with the linear configuration of HCN as can be seen in Figure 5.1 while in the RBAA we averaged over the bending angles and then included contributions from less attractive parts of the PES. Using these bound level energies, we determined the transitions frequencies and compared them with the experimental data available for this system. The RBAA method was found to give a slightly better agreement with experiment than the RMA.

In the same paper, the cross section for several rovibrational transitions were also computed using this RBAA and RMA approaches. The very similar curves representing the cross section obtained show that the RMA approach is quite satisfactory for the computation of pure rotational transitions of the HCN molecule. In the second paper, the cross sections for several rovibrational transitions were computed also using this RB-CC approach. A comparison of the RMA, RBAA and RB-CC results shows that the pure rotational transition in the ground vibrational level can be accurately studied at the RMA level while the transitions involving different vibrational levels should be studied using the RB-CC method as the RBAA fails to give an accurate estimate of the magnitude of these cross sections. In Figure 5.3, we present the rate coefficient for several transitions between 4 ≤ T ≤ 100 K. We included the l-doubling transition ν = 2 → ν = 1 for j = 1. The cross section used for determine these rate coefficient were calculated at the RB-CC level. We found that the rate coefficients associated with the l-doubling transitions are of the same order of magnitude than the pure rotational transitions in the ground vibrational state in this domain of temperature.

Before concluding this summary, it is worth mentioning that the z -axis which is used in the calculations of the HCN rovibrational wave functions differs by an angle α from the z-axis used to define the θ angle for the dynamics (see Fig 5 .2). This change of axis was called dynamic axis switching by Band and Freed in their study of photodissociation of linear molecule [91]. In the present work, we have taken α = 0 which is valid only for the linear configuration of HCN. When we average over the bending angle, as it is done in the RBAA method, this identity is verified on the average. However at any given instantaneous displacement from linearity the axis system differs slightly by the angle α.

α(γ) = arctan     sin γ 1 + CG CN CH 2 -2 CG CN CH cos γ     , (5.8)
For γ in the interval [120,180] • , this function is well approximated by the linear law α = 0.64 [180γ]. This means that the value of α which is zero for γ = 180 • reaches a value of 12 degrees for γ = 160. As bending vibration is generally small in amplitude in the lowest bending states the previous neglect seems to be a reasonable approximation if we consider the lower bending states of a triatomic molecule which equilibrium geometry is linear.

Publications

THE JOURNAL OF CHEMICAL PHYSICS 139, 034304 (2013)

The interaction of He with vibrating HCN: Potential energy surface, bound states, and rotationally inelastic cross sections

Otoniel Denis-Alpizar, 1,2 Thierry Stoecklin, 1,a) Philippe Halvick, 1,b) and Marie-Lise Dubernet 3,4 A four-dimensional potential energy surface representing the interaction between He and hydrogen cyanide (HCN) subjected to bending vibrational motion is presented. Ab initio calculations were carried out at the coupled-cluster level with single and double excitations and a perturbative treatment of triple excitations, using a quadruple-zeta basis set and mid-bond functions. The global minimum is found in the linear He-HCN configuration with the H atom pointing towards helium at the intermolecular separation of 7.94 a 0 . The corresponding well depth is 30.35 cm -1 . First, the quality of the new potential has been tested by performing two comparisons with previous theoretical and experimental works. (i) The rovibrational energy levels of the He-HCN complex for a rigid linear configuration of the HCN molecule have been calculated. The dissociation energy is 8.99 cm -1 , which is slightly smaller than the semi-empirical value of 9.42 cm -1 . The transitions frequencies are found to be in good agreement with the experimental data. (ii) We performed close coupling calculations of the rotational de-excitation of rigid linear HCN in collision with He and observed a close similarity with the theoretical data published in a recent study. Second, the effects of the vibrational bending of HCN have been investigated, both for the bound levels of the He-HCN system and for the rotationally inelastic cross sections. This was performed with an approximate method using the average of the interaction potential over the vibrational bending wavefunction. While this improves slightly the comparison of calculated transitions frequencies with experiment, the cross sections remain very close to those obtained with rigid linear HCN. © 2013 AIP Publishing LLC.

[http://dx.doi.org/10.1063/1.4813125] I. INTRODUCTION The rigid monomer approximation (RMA), where the monomer's geometry is assumed to be independent of the dimer configuration, is commonly used to simulate the dynamics of systems governed by weak intermolecular interaction and where no breaking or formation of chemical bonds takes place. The decoupling of intramonomer and intermonomer motions reduces the dimensionality and thus simplifies greatly the calculation of the dynamics. Intermolecular bound states or cross sections for low collision energies can be calculated within this approximation. The quality of the approximation can be improved by using the average of the intramonomer coordinates over the internal stretching motions. 1 When the ratio of intramonomer over intermonomer vibrational frequencies is large (about 100), the RMA is very reliable. This has been demonstrated 2, 3 by the excellent agreement between experiment and calculation of the infrared spectrum of the H 2 -CO complex. For the same system, a good agreement has also been obtained between calculations and the first low temperature experimental inelastic cross section. 4 a) Electronic mail: thierry.stoecklin@u-bordeaux1.fr b) Electronic mail: philippe.halvick@u-bordeaux1.fr However, in the case of a triatomic (or larger) monomer, the RMA can be questioned because the coupling between the internal bending motion and the intermonomer motion may not be negligible. Bending motion may have large amplitude and low frequency, inducing a significant change of the electronic cloud, and consequently, a significant change of the intermolecular forces. While the RMA can be questioned for very floppy monomer (e.g., C 3 ), it is not known if this method can be accurate for rigid or semi-rigid monomers with vibrational bending mode.

Hydrogen cyanide (HCN) and isocyanide (HNC) are among the most abundant organic molecules in the interstellar medium. Owing to a large dipole moment, both molecules decay fast in their rotational energy ladder. The rotational emission lines of HCN and HNC are considered to be a major tracer of dense molecular gas (star-forming molecular clouds) in luminous and ultraluminous infrared galaxies. [5][6][7][8] Rotationally excited HCN and HNC suggests an excitation mechanism fast enough to counter the decay, such as frequent collisions with He and H 2 in dense clouds. Consequently, the estimation of abundances of both isomers in the interstellar clouds has motivated theoretical studies of the rotational excitation in collisions with He 9-12 and H 2 . 13 In these studies, the HCN or HNC molecule was always considered as a linear rigid rotor.

However, vibrational excitation of HCN has been observed in the interstellar medium. The rotational transitions of vibrationally excited HCN have been used to probe 14 the dust formation region around the carbon-rich star IRC +10216. The high vibrational levels are populated by radiation and by collision, owing to the high temperature, high gas density and high radiation flux prevailing in the circumstellar envelope. Vibrationally excited HCN in the ν 2 = 1 state has been also observed 15 in the nucleus of the luminous infrared galaxy NGC 4418. Most likely, the molecule is pumped to the excited level by infrared radiation and returned to the vibrational ground state with rotational excitation. 16,17 These observations suggest that the vibrational excitation of HCN, at least in the bending motion, deserves to be considered in the collision mechanisms.

The first studies dedicated to the rotational excitation of rigid linear HCN (l-HCN) by collisions with He atoms were based on the potential energy surface (PES) of Green and Thaddeus. 9 This primitive PES was obtained using the uniform electron gas model. Several new intermolecular potentials were later published in the last 20 years for the l-HCN-He system. Drucker et al. 18 calculated one at the MP4 level and reported the first theoretical determination of the highresolution microwave and millimeter spectrum. Later Atkins and Hutson 19 obtained two empirical PESs based on two different functional forms using the experimental data available. Toczyłowski et al. 20 reported a theoretical PES calculated at the CCSD(T) level, hereafter denoted by S01, which was found to describe correctly the internal-rotational band measured by Drucker et al. 18 and with a global minimum of -29.90 cm -1 . The most recent studies of the rotational excitation of l-HCN by He done by Sarrasin et al. 11 and Dumouchel et al. 12 used this last surface. The latest PES published for the l-HCN-He system is a semi-empirical one by Harada et al. 21 , denoted by S02, which was obtained by modifying the S01 surface in order to reproduce the experimental transitions frequencies.

The present paper focuses on the development of a PES describing the collision between He and HCN considered as a rigid bender. The vibrational bending motion of HCN is treated quantally while the CH and CN bond lengths are set to constant values. As a first test of this new PES, we determined the rovibrational energy levels of the l-HCN-He system and compared it to the existing theoretical and experimental data. We also computed the l-HCN-He inelastic cross sections and compared it with the theoretical data of Sarrasin et al. In the second part of this work, the effects of the vibrational bending of HCN have been investigated by using an interaction potential averaged on the bending wavefunctions. 22,23 Again, we calculated the energies of the rovibrational bound states and the inelastic cross sections and we compared these last results with the previous ones.

II. AB INITIO CALCULATIONS AND POTENTIAL FUNCTIONAL FORM

The body-fixed coordinates used in this work are shown in Fig. 1. R, θ , and ϕ are the intermonomer coordinates which describe the relative positions of the HCN molecule and He atom, while γ is the intramonomer coordinate which describes the bending angle of HCN. R is the distance from the center of mass of the HCN to the He atom. ϕ is the angle of rotation around the axis defined by the H atom and the center of mass of CN. θ is the angle between the latter axis and the axis defined by the He atom and the center of mass of HCN. The C-H and C-N rigid bond lengths have been set to the experimental values 24 (r CH = 2.0135 a 0 , r CN = 2.1792 a 0 ) at which we have added the corrections for the averaging over the ground vibrational state. 25 This results to r CH = 2.0286 a 0 and r CN = 2.1874 a 0 .

The interaction potential of HCN with He has been calculated in the framework of the supermolecular approach with the coupled-cluster method with single and double excitations and a perturbative treatment of triple excitations (CCSD(T)). The interaction energy was corrected at all geometries for the basis set superposition error (BSSE) with the counterpoise procedure of Boys and Bernardi. 26 A comparison of the interaction energies calculated with basis sets 27 of triple, quadruple, and quintuple-zeta quality is shown in Table I, with or without an additional set of bond functions 28 centered at middistance between the He atom and the HCN center of mass. The interaction energy, calculated at a configuration close to the equilibrium geometry, is quite stable in respect of the size of the basis set and the use of bond functions. For the largest basis set, it is safe to assume that the convergence of the oneelectron basis is close to the complete basis set limit. Considering the computational cost associated with the various basis sets, we have chosen the quadruple zeta basis set with bond functions. TABLE I. CCSD(T) interaction energy of the l-HCN-He system at R = 7.97 a 0 and θ = 0 • . The use of bond functions is denoted by +bf.

Basis set

Energy (cm The interaction energy was computed over a dense fourdimensional grid of points defined by the product of four onedimensional grids associated to a single coordinate. The radial grid included 35 points ranging from 3.8 a 0 to 20.8 a 0 . The bending grid included 11 points between 180 • and 110 • . The angular grids were spaced uniformly in steps of 10 • for θ and 30 • for ϕ, both in the range [0 • ,180 • ]. The total number of points was 43015. All calculations were carried out with the MOLPRO package. 29 The ab initio energies were fitted to a parameterized functional form defined as a sum of a short-range and a long-range contribution,

V int (R, θ, ϕ, γ ) = S(R) 14 l=0 min (l,3) m=0 F SR lm (R, γ ) Plm (θ ) cos (mϕ) +(1 -S(R)) 5 l=0 min (l,3) m=0 F LR lm (R, γ ) Plm (θ ) cos (mϕ). (1)
Here, Plm are normalized associated Legendre polynomials. F SR lm , F LR lm , and S are the short-range radial functions, the long-range radial functions, and the switching function, respectively,

F SR lm (R, γ ) = e -αR 9 n=0 R n 3 j =0 C lmnj Pj (cos γ ), (2)
F LR lm (R, γ ) = 8 k=6 t k (βR) R k 3 j =0 D lmkj Pj (cos γ ), (3)
S(R) = 1 2 [1 -tanh (A 0 (R -R 0 ))], (4)
where Pj are normalized Legendre polynomials and t k is the Tang-Toennies damping function,

t k (x) = 1 -e -x k i=0 x i i! . ( 5
)
The nonlinear parameters α, β, A 0 , and R 0 were set to the values α = 1.91 a 0 -1 , β = 1.06 a 0 -1 , A 0 = 1.69 a 0 -1 , and R 0 = 10.58 a 0 . These values were determined by the trial and error method. The linear parameters C lmnj and D lmkj were calculated with the weighted linear least squares method. On each ab initio point, we applied a weight w depending both on the interaction energy E and the angle γ ,

w = γ 0 (τ -γ ) 2 min 1, V 0 |E| , ( 6
)
with V 0 = 1000 cm -1 , γ 0 = 100 • , and τ = 181 • . Let us note that the ab initio grid is restricted to γ ≥ 110 • . Indeed, the rigid bender approximation used for HCN is expected to be reliable only for the ground and the first excited bending states, and possibly for the second excited state. Moreover, the potential energy of the HCN molecule at γ = 120 • is 7130 cm -1 . This value is much larger than the energy at which the rigid bender approximation remain reliable, if we remind that ω 2 is slightly larger than 700 cm -1 . Therefore, because there is no need to represent the interaction energy for γ ≤ 120 • , this value is used as a cut-off limit. Below this limit, the interaction energy is set equal to its value at γ = 120 • .

The total PES is the sum of the interaction energy of the He-HCN complex plus the one-dimensional bending energy of the isolated HCN molecule. The latter was calculated with the same ab initio method and same basis set which were used for the former. The ab initio calculations have been performed for a grid of 22 values of the bending angle γ in the range of 50 • -180 • , with the same C-H and C-N rigid bond lengths as defined previously. These bending potential energies were fitted to a linear combination of four Legendre polynomials.

III. BOUND STATES AND SCATTERING CALCULATIONS

We used the close coupling method to calculate both the rovibrational energy levels and the inelastic cross section of the He-HCN system. The coupled equations needed for scattering calculations are identical to those for bound states, the only difference being the applied boundary conditions. In this study we compare two approaches. In the first one, the bending motion is completely neglected and we use only the linear configuration of HCN and perform usual atom linear molecule calculations using for HCN a rigid rotor description. In the second one, we fix the value of ϕ to 0 as the potential varies slowly as a function of this angle and we calculate for each value of the intermolecular coordinate R used in the dynamics calculations the following expansion of the interaction potential in a Legendre polynomial P l (cos θ ) basis set along a grid of the bending angle γ ,

V int (R, θ, ϕ = 0, γ ) = l D l (R, γ ) P l (cos θ ). (7)
We then calculate the rigid bender energies and wavefunctions of HCN in internal coordinates using the bending potential of HCN described in Sec. II and the Hamiltonian of Carter and Handy, 30

H J =0 RB = - ¯2 2 1 μ 1 R 2 1 + 1 μ 2 R 2 2 ∂ 2 ∂γ 2 + cot γ ∂ ∂γ - ¯2 2M C R 1 R 2 ∂ 2 ∂γ 2 + cot γ ∂ ∂γ cos γ + cos γ ∂ 2 ∂γ 2 + cot γ ∂ ∂γ + l C l P l (cos γ ) , (8)
where

1 μ1 = 1 MH + 1 MC , 1 μ2 = 1 MC + 1
MN , and R 1 and R 2 are respectively the CH and CN bond lengths.

It may be confusing to compare this Hamiltonian with the different Hamiltonians published at that time [30][31][32] as some other terms are present in some of the three references (sometime with different signs) and are not in others. This is probably because the matrix elements of these missing terms in a Legendre polynomial basis set do compensate each other. The matrix elements of this rigid bender Hamiltonian in a Legendre polynomial basis set are

P l |H J =0 RB |P k = ¯2 2 1 μ 1 R 2 1 + 1 μ 2 R 2 2 δ kl 2l (l + 1) (2l + 1) - ¯2 2M C R 1 R 2 2l 3 > δ k,l±1 (2l > + 1) (2l > -1) +2 n C n l n k 0 0 0 2 , ( 9
)
where l > = max(l,k). The diagonalisation of this matrix gives the rigid bender energies n and wave functions χ n (γ ) as a function of the bending angle for the HCN rotational angular momentum j = 0. We take the same bending wavefunctions for all the values of j since the variation of the bending wavefunctions as a function of j is expected to be weak, at least when j is not too large. The wavefunctions describing the HCN motion within this very simple approach are then the product of a bending wavefunction by a spherical harmonics describing the rotation. Consequently, the energies of HCN are

E nj = B H CN j (j + 1) + n . (10)
For each value of the intermolecular coordinate R considered in the course of the propagation, the coefficients calculated in (7) are then averaged over the bending wavefunctions,

χ n |V int (R, θ, ϕ = 0, γ )| χ m = l dγ {χ n (γ ) D l (R, γ ) χ n (γ )} P l (cos θ ) = l Dn,m l (R) P l (cos θ ). ( 11
)
The problem is now formally equivalent to an atom colliding with a fictitious vibrating diatomic molecule where the vibration of the diatomic molecule is in fact the bending vibration. Using this very simple approach denoted in the following Rigid Bender Averaged Approximation (RBAA), we can obtain state to state cross sections for the transition between two different bending and rotational levels of HCN as well as bending averaged energies for the He-HCN complex. We use our NEWMAT code both for the scattering and the bound states calculations. This is a close coupling code working in the spaced fixed frame which has been described in some of our recent works. 33,34 The rotational basis set for HCN included 20 functions and the rotational constant of HCN was set to its experimental value 35 B HCN = 1.478 22 cm -1 . The maximum propagation distance was 80 a 0 and two values of the propagator step size (0.05 a 0 and 0.01 a 0 ) were used for the bound state calculations. The final bound state energies of the He-HCN complex were obtained from a Richardson extrapolation as recommended by Hutson. 36

IV. RESULTS AND DISCUSSION

A. Potential energy surface

The functional form defined above allowed us to obtain an accurate representation of the PES. The root mean square (RMS) of the differences between the ab initio and the interpolated total potential energies E is 0.016 cm -1 for the energies E ≤ 0 cm -1 (E = 0 cm -1 corresponds to the infinite separation of monomers). For 0 < E ≤ 1000 cm -1 , the RMS of the relative errors is below 1%, and for 1000 < E ≤ 3000 cm -1 , it is about 2%.

Contour plots of the interaction PES are shown in Fig. 2 for several values of the bending angle γ and for ϕ fixed at 0 • . The selected values of γ are lying in the range assumed to be spanned by the first excited vibrational function. We observe that the bending of HCN has a visible effect in the bottom of the potential well and in the repulsive short-range interaction. The long-range part of the potential is hardly changed by the bending. For γ = 180 • , the potential is, by definition, isotropic versus ϕ. In the range 150 • ≤ γ ≤ 180 • , the potential remains nearly isotropic in respect of ϕ (not shown here), except in the short-range repulsive region. Contour plots for ϕ = 180 • are shown in Fig. 3. For the same value of γ , the comparison with the contour plots at ϕ = 0 • shows a significant change only for γ = 150 • and R ≤ 7 a 0 .

The global minimum of the total PES has a depth of 30.35 cm -1 and a linear He-HCN configuration: γ = 180 • , θ = 0 • , and R = 7.94 a 0 . It is 0.45 cm -1 deeper than for the S01 PES. The latter was calculated with a triple-zeta basis set, while we use here a quadruple-zeta one. The discrepancy observed in the well depth is mainly a consequence of the basis set quality, and this is confirmed by the data shown in Table I. The difference due to the different analytical representations is probably not significant. Moreover, the present PES has a well depth of only 0.15 cm -1 deeper than the one of the semi-empirical surface S02, which is a S01 PES modified in order to improve the agreement with the experimental millimeter-wave spectrum. A secondary minimum with a depth of 22.08 cm -1 and a bent configuration is found at γ = 180 • , θ = 110.4 • , and R = 6.78 a 0 . This secondary minimum is very similar to the global minimum of the He-CN PES. 37 The FORTRAN 90 code of the PES is available on request by electronic mail to the authors.

B. Bound states and spectrum

The bound levels calculated in the RMA and RBAA with the present PES are collected in Table II. The approximate rotational quantum number of HCN and orbital quantum number, j and l, respectively, are also reported in this table. The

RMA

RBAA S02 a State l J ε Energy (cm -1 ) Energy (cm -1 ) Energy (cm -1 ) energies calculated using the RBAA are systematically above those obtained using the RMA. This is not surprising as the most attractive bending angle is associated with the linear configuration of HCN. The maximum value of the total angular momentum J leading to bound states is 5 in both cases. The potential well supports 19 bound levels and the dissociation energy is 8.986 cm -1 . Harada et al. 21 obtained a larger dissociation energy of 9.420 cm -1 using the S02 PES, which was optimized in order to reproduce the experimental transitions frequencies. All the bound state energies calculated by Harada et al. 21 are lower than those calculated in the present work by about a half cm -1 and they obtained three more bound states.

ν s = 0, j = 0 0 0 + - 8
A comparison between the PES presented in this work and the S02 PES is shown in Fig. 4. While the depth of the global potential well (θ = 0 • ) is almost the same in both PES, the S02 PES is deeper by more than 1 cm -1 for the other angular configurations. This allows the bound levels of the S02 PES to be slightly lower than those of the present PES.

The calculated transitions frequencies using the RMA and RBAA approaches are compared in Table III with the spectroscopic data available. 18,21 Harada et al. 21 reported most of the Q-and R-branch lines including the splitting into several hyperfine components due to the spin angular momentum of the nitrogen nucleus (I = 1). As our calculation do not include the hyperfine structure and because the spin splitting is very small in comparison with the spacings of the rotational lines, we compare our results with those of Harada et al. 21 averaged over the hyperfine components. The agreement between our results and experiment is quite good, with a difference of less than 3.2% in the RMA in all cases, with the exception of the transition at 4604 MHz for which the error is about 13%. This is however better than the (∼30%) error obtained by Toczyłowski et al. 20 for this line while Harada et al. 21 did not mention it. The agreement between our results and experiment is even better when using the RBAA approach as the maximum error is now less than 2.4% again with the exception of the transition at 4604 MHz for which the error is about 18.6%. This is the only transition which for the error is increased when using the RBAA.

The transition ( j = 1←0)R(4) reported by Harada et al. 21 is missing in our comparison as it involves the upper state (j, l, J) = (1,4,5) which was not found to be bound using our PES.

For each couple (j = 1,l = n) with n ≥ 1, there are three levels (J = n-1,n,n+1) which are very close in energy. In the present calculation, the states (1,4,3) and (1,4,4) have the energies -0.532 and -0.096 cm -1 , respectively. Consequently, it is not possible for the third state (1,4,5), which is expected to be lying about ∼0.4 cm -1 above the state (1,4,4), to be bound. With a potential well deeper by about a half cm -1 or with a slightly more attractive long range interaction or less repulsive short-range interaction, the missing state (1,4,5) could appear in the calculations.

C. Inelastic cross sections

The inelastic cross sections were first calculated in the RMA in order to compare with the previous work. 11 Figure 5 shows the de-excitation cross sections for the first rotational levels. The shape, the positions, and the amplitudes of the resonances supported by the van der Waals well which appear on this figure are accurate fingerprints of the PES used in the calculations. We do not intend here to analyse the characteristics of these resonances which are typical of van der Waals systems and have been discussed in detail for similar systems by several authors. 38 We simply compare our results with those of Sarrasin et al., 11 obtained using the S01 PES. A very close similarity is observed between the latter cross sections and the ones presented in Fig. 5, indicating that the S01 PES and the present PES, restricted to the rigid linear HCN configuration, are very similar.

The propensity for even j transitions which can be observed in Fig. 5 was previously examined by Sarrasin et al. 11 This is related to the shape of the interaction potential. The radial even Legendre expansion coefficients are significantly larger than the odd ones.

Then we investigated the bending dependence of the cross section with the present PES, by computing the same rotational transitions using the RBAA approach. This approach does not include exactly the coupling between vibration and rotation which will be the object of a future work but allows checking significant variations of the PES as a function of the bending angle. These results are compared to those obtained using the RMA approach. The elastic cross sections, which are not represented, are almost unchanged. In Fig. 5, the inelastic cross sections for the | j| = 2 transitions are also very close for both approaches. But discrepancies can be observed for the | j| = 1 transitions below the collision energy of 1 cm -1 . There, the RMA cross sections are two or three times larger or smaller than the RBAA cross sections. Nevertheless, we can conclude that the RMA approach is quite satisfactory above 1 cm -1 to calculate rotational excitation cross sections for a linear triatomic molecule like HCN whose bending vibration frequency 39 is relatively small (711.98 cm -1 ) but still large compared to the rotational constant 35 (1.47822 cm -1 ).

V. CONCLUSION

We presented the first theoretical study of the collision of HCN with He including the bending vibration of HCN. We calculated a four dimensional analytical representation of the PES based on supermolecular ab initio calculations using a quadruple zeta basis set with mid-bond functions and BSSE correction. The van der Waals well was found to be 30.35 cm -1 deep and associated with the linear configuration (He-HCN) while a secondary minimum with a depth of 22.08 cm -1 , associated with a bent configuration, was also identified. Bound states calculations were performed using this PES. The results are in good agreement with the available experimental data. We checked that the restriction of the dynamics to the rigid linear configuration of HCN gives similar close coupling inelastic cross section than the previous theoretical works. We also presented a simple method (RBAA) of calculation of the rotational close coupling cross section which uses the average of the interaction potential over the bending wave functions of HCN. We found that taking into account the bending motion through the RBAA method does not change significantly the rotational excitation cross sections, while the agreement of the calculated bound state transition frequencies with the experiment is marginally improved. This first study shows in any case that the RMA approach is quite satisfactory for the computation of rotational excitation cross sections for a linear triatomic molecule like HCN. The same accuracy could also be expected for other rigid or semi-rigid triatomic (and larger) molecules discovered in the interstellar medium. This finding is particularly useful if we consider the calculations of rotational transitions of polyatomic molecules in collision with H 2 which are very computationally demanding. Nevertheless, this preliminary conclusion needs to be confirmed by a comparison of the RMA approach with accurate calculations using an Hamiltonian which includes the exact vibrotational coupling. Efforts in that direction are in progress.

I. INTRODUCTION

Because of their importance to model the chemistry of interstellar clouds, quantum inelastic scattering calculations involving small polyatomic molecules are the subject of many theoretical studies. 1 However most of them are limited to the use of the rigid rotor approximation as it is expected to be a quite accurate approach to calculate rotational transitions which are the most probable at the typical temperature of interstellar clouds. This is the case of most available studies of energy transfer collisions involving an atom and a linear 2 or bent triatomic molecule. 3 The neglect of the vibrational motion relies on the fact that the vibrational frequencies are commonly quite large compared to the rotational ones. The validity of this approximation was recently confirmed at low collision energy by the excellent agreement between experiment and calculation obtained for the H 2 -CO inelastic cross sections. 4 However, several authors pointed out that infrared transition among the molecular vibrational levels could significantly increase the intensities of the rotational transitions by populating the upper rotational levels. 5 The question of the validity of this approximation was then recently investigated by the group of Alexander and Dagdigian for CH 2 . 6 They calculated the bending levels of CH 2 and averaged the He-CH 2 potential over the bending angle using these functions in order to reduce the problem to a fictitious atom-rigid asymmetric top molecule collision. In a recent paper, 7 hereafter denoted Paper I, we used a similar approach which we called Rigid Bender Average Approximation (RBAA) for a collision involving this time He and a linear molecule: HCN. This last molecule and isocyanide (HNC) are among the most abuna) Author to whom correspondence should be addressed. Electronic mail: t.stoecklin@ism.u-bordeaux1.fr dant organic molecules in the interstellar medium. The averaging over the bending angle of the intermolecular potential reduced the calculation in our case to even simpler atom linear molecule calculations. The bending frequency of HCN is about a factor of 2 smaller than that of CH 2 suggesting that the exact treatment of the coupling between rotation and bending may be more important for this molecule. Therefore, the present study is dedicated to the development of a method including exactly the coupling between bending and rotation for a collision between an atom and a rigid bender. This study was furthermore motivated by recent astrophysical measurements of vibrationaly excited HCN 8,9 in the interstellar medium.

The authors assumed that the molecule is pumped to the excited bending level by infrared radiation and return to the vibrational ground state with rotational excitation. 10,5 What makes the exact treatment even more necessary is the recent detection of HCN l-type transitions in hot planetary nebula 11 as it involves nearly degenerate levels. Such transitions can be calculated using the RBAA approach but it seems important to check if the results given by this approximation are reliable by comparing with the exact results. In Sec. II of the present paper we give the close coupling equations for atom rigid bender collisions. In Sec. III, we compare our results to those given by the RBAA approach and for pure rotational transitions with atom-linear HCN close coupling results.

II. METHOD

We use the set of coordinates illustrated in Fig. 1 which was defined in Paper I dedicated to the calculation of the He-HCN rigid bender potential energy surface and we first calculate the bending levels of HCN alone for each value of its rotational angular momentum. While the problem of calculating the rovibrational states of a triatomic molecule in internal coordinates was solved long ago and can be performed using the MORBID code of Jensen 12 we decided to give here a short overview of the corresponding equations which we coded as the literature of that time is full of contradicting terms and misprints.

A. Calculation of the HCN rigid bender energies and wave functions in internal coordinates

We solve the variational problem using the rigid bender Hamiltonian of Sutcliffe 13 which for the Z molecular axis is along the intramolecular coordinate between H and the center of mass of CN. This effective Hamiltonian is obtained by fixing the CH and CN distances to their equilibrium values in the linear configuration of HCN:

H = K (1) V + K (2) V + K V R + V (γ ), (1)
where K (1) V and K (2) V are the contributions to the kinetic operator which are independent of rotation while K V R is the rotation-bending interaction term. The bending angle is denoted as γ and V(γ ) is the one-dimensional bending potential of HCN reported in Paper I,

K (1) V = - ¯2 2 1 μ 1 R 2 1 + 1 μ 2 R 2 2 ∂ 2 ∂γ 2 + cot γ ∂ ∂γ , (2)
K (2) V = ¯2 μ 12 R 1 R 2 cos γ ∂ 2 ∂γ 2 + cot γ ∂ ∂γ + sin γ ∂ ∂γ , (3)
K V R = 1 2 1 μ 1 R 2 1 cot 2 γ + 1 μ 2 R 2 2 1 sin 2 γ - 2 cos γ sin 2 γ 1 μ 12 R 1 R 2 ˆ 2 Z + 1 2μ 1 R 2 1 ˆ 2 x + ˆ 2 y + 1 2μ 1 R 2 1 cos γ sin γ - 1 2μ 12 R 1 R 2 1 sin γ [ ˆ X ˆ Z + ˆ Z ˆ X ] + ī 1 μ 1 R 2 1 ∂ ∂γ + 1 2 cot γ - 1 2μ 12 R 1 R 2 1 sin γ + 2 cos γ ∂ ∂γ ˆ Y , ( 4
)
with 1 μ1 = 1 MC + 1 MN , 1 μ2 = 1 MC + 1 MH , μ 12 = M C .
And where R 1 and R 2 are respectively the CN and CH interatomic distances while ˆ x , ˆ y , ˆ z are the projection over the molecule fixed axis of the rotational angular momentum of the HCN molecule.

We then follow Sutcliffe and Tennyson 14 and first take the matrix elements of H in a symmetric top basis set |jKM = 2J +1 4π D j * M,K (φ RHCN , θ RHCN , 0) where M and K are the projections of the rotational angular momentum j of HCN along the Z space fixed axis and along the Z molecular fixed axis, respectively, and where (φ RHCN , θ RHCN ) are the spherical coordinates of the vector R H CN joining H to the center of mass of CN in the space fixed frame:

jK M|K V R |jKM = δ K,K 1 2μ 1 R 2 1 [j (j + 1) -2K 2 ] + K 2 2 1 sin 2 γ 1 μ 1 R 2 1 + 1 μ 2 R 2 2 - 2 cos γ μ 12 R 1 R 2 + δ K+1,K C + jK 1 2μ 1 R 2 1 - ∂ ∂γ + K cot γ + δ K-1,K C - jK 1 2μ 1 R 2 1 ∂ ∂γ + K cot γ + δ K+1,K C + jK 1 2μ 12 R 1 R 2 cos γ ∂ ∂γ -K cot γ -K sin γ + δ K-1,K C - jK 1 2μ 12 R 1 R 2 -cos γ ∂ ∂γ + K cot γ -K sin γ , ( 5
)
with C ± J K = √ J (J + 1) -K(K ± 1). The second term of this expression cancels some other terms appearing in K V (1) + K V (2) . The terms remaining in K VR are split into the two following expressions: jK M|K (1) V

R |jKM = δ K,K ¯2 2μ 1 R 2 1 [j (j + 1) -2K 2 ] + δ K+1,K C + jK ¯2 2μ 1 R 2 1 - ∂ ∂γ + K cot γ + δ K-1,K C - jK ¯2 2μ 1 R 2 1 ∂ ∂γ + K cot γ (6)
and

jK M|K (2) V R |jKM = ¯2 2μ 12 R 1 R 2 δ K+1,K C + jK cos γ ∂ ∂γ -K cot γ -K sin γ + δ K-1,K C - jK -cos γ ∂ ∂γ + K cot γ + K sin γ . ( 7
)
A normalized associated legendre polynomial basis set P k l (cos γ ) as defined by Green 15 is then used to describe the bending vibration.

Y K l (γ, ϕ) = (-) K P K l (cos γ )e iKϕ ,
and we obtain for the different terms of the Hamiltonian (1) the following matrix elements:

jK M| P K l (cos γ ) K (1) V P K l (cos γ )|jKM = ¯2 2 δ K,K δ l,l 1 μ 1 R 2 1 + 1 μ 2 R 2 2 , ( 8
)
jK M P K l (cos γ ) K (2) V | P K l (cos γ )|jKM = -¯2 μ 12 R 1 R 2 δ K,K δ l ,l±1 (l > + m)(l > -m) (2l > + 1)(2l > -1) l 2 (l + 1) 2 . ( 9
)
In this last expression the top and the bottom terms enclosed in curly brackets are respectively associated with the values l = l + 1 and l = l -1 while l > is the maximum of l jK M| P K l (cos γ ) K (1) V R

P K l (cos γ )|jKM = ¯2 2μ 1 R 1 δ l ,l [δ K,K [j (j + 1) -2K 2 ] -δ K ,K+1 C + J K C + lK -δ K ,K-1 C - J K C - lK ], ( 10
)
jK M| P K l (cos γ ) K (2) V R P K l (cos γ )|jKM = ¯2 2μ 12 R 1 R 2 δ K ,K+1 C + jK × C + lK δ l ,l±1 (l > + K + 1)(l > -K -1) (2l > + 1)(2l > -1)
-K (a lK δ l ,l+1b lK δ l ,l-1 )

+ δ K ,K-1 C - jK C - lK δ l ,l±1 (l > -K +1)(l > + K -1) (2l > + 1)(2l > -1)
+ K(a l-K δ l ,l+1b l-K δ l ,l-1 ) , (11) where

a lK = (l + K + 1)(l + K + 2) (2l + 3)(2l + 1)
and

b lK = (l -K -1)(l -K) (2l -1)(2l + 1) .
The potential is expanded in Legendre polynomials:

V (γ ) = L C L P L (cos γ ), ( 12
)
and we obtain the following expression for the matrix elements of the potential:

jK M| P K l (cos γ )|V (γ )| P K l (cos γ )|jKM = δ K ,K L C L (-) K (2l + 1)(2l + 1) × l L l 0 0 0 l L l -K 0 K . ( 13
)
For each value of j, the diagonalisation of the resulting hamiltonian matrix gives the rigid bender energies ε νj and the corresponding space fixed rigid bender eigenfunctions

χ jM ν (γ ) = K n≥K C νj n,K P K n (cos γ )|jKM = K ν j,K (γ )|jKM , ( 14
)
where ν designates the bending quantum number. This equation shows that the rovibrational functions obtained from the rigid bender model can be put in the form of a product of an asymmetric top rotational function by a vibrational bending function ν j,K (γ ). This is in agreement with the intuitive representation that the bending vibration of a linear triatomic molecules makes it become an instantaneous asymmetric top, the square of the ν j,K (γ ) function giving the probability of a given bending angle for a specific values of the quantum numbers associated with the bending state ν the rotational state j and its projection along the molecular axis K. The calculation of the rigid bender energies and wave functions can be simplified by using the parity defined basis set suggested by Sutcliffe and Tennyson: 14

|jn KM, p = 1 2(1 + δ K0 ) [|jn KM + (-) p |jn -KM ] (15)
with

|jn KM = P K n (cos γ )|j KM . ( 16
)
Knowing that the action of the parity operator on a symmetric top function is given by |j KM = (-) (j -K) |j -KM and that P -K n (cos γ ) = (-) K P K n (cos γ ) (see, for example, Zare 16 ), it is easily seen that the parity of the function ( 14) is equal to (-) j+p with p = 0 or 1 and K is the absolute value of K. The use of this basis set gives parity selected rigid bender energies ε p νj and eigenfunctions

χ jM νp (γ ) = K≥0 n≥K Cνjp n, K |jn KM, p = K ˜ νp j, K (γ ) 2(1+δ K0 ) [|jn KM +(-) p+K |jn-KM ] (17)
with

˜ νp j, K (γ ) = n≥ K Cνjp n, K P K n (cos γ ). ( 18
)
The matrix elements of the triatomic rigid bender Hamiltonian in this symmetrised basis set are readily obtained from those given in expressions ( 9)-( 12) from the simple rule:

j n K M, q|H |jn KM, p = δ pq f (p, q) j K M|P K n (cos γ )|H | P K n (cos γ )|j KM (19)
with

f (p, q) = 1 forK = K = 0 or K , K > 0 and f (p, q) = √ 2 for K = 0, K > 0 or K > 0, K = 0.

B. Close coupling equations

The theory of the inelastic scattering of two rigid polyatomic molecules was developed long ago 17 but studies including both vibration and rotation of the fragments are seldom. Most of such studies use the pioneering approaches developed by Clary 18 or more recently by Bowman 19 which rely on the use of the Infinite order sudden approximation for the rotation and the Close Coupling for the vibration (VCC-IOS). We treat here the scattering of a rigid bender molecule colliding with an atom. This approach takes only into account the coupling between bending and rotation but could be easily generalized to include the stretching vibrations. We use the result of the previous paragraph to express the rovibrational wave function of the He-HCN complex in space fixed coordinates as

|νj lJ M = mj ml 2J + 1 4π jm j lm l |J M Y ml l ( R) × K ν j,K (γ )D j * mj ,K (φ RHCN , θ RHCN , 0), ( 20
)
where R are the spherical coordinates of the intermolecular vector R between He and the center of mass of HCN in the space fixed frame. Following Green 20 we now expand the intermolecular potential between the atom and the rigid bender molecule in the body fixed coordinates defined in Fig. 1 of Paper I like

V A-RB (R, γ, θ, ϕ) = λ μ v λμ (R, γ )Y μ λ (θ, ϕ), ( 21
)
which can be written:

V A-RB (R, γ, θ, ϕ) = λ μ≥0 v λμ (R, γ )[2 -δ μ0 ] P μ λ (cos θ ) cos(μϕ). ( 22
)
Expression (21) transformed by rotation in terms of space fixed angles reads:

V A-RB R, γ, φ RHCN , θ RHCN , R = λ μ v λμ (R, γ ) ν D λ * μν φ RHCN , θ RHCN , 0 Y ν λ ( R). (23)
The matrix elements of the interaction potential V A-RB between the atom and the rigid bender molecule in the space fixed coordinates are obtained from straightforward algebra:

[V A-RB ] J M νj l,ν j l (R) = νj lJ M|V A-RB |ν j l J M = π 0 dγ sin γ K K ν j,K (γ )W J M jKl,j K l (R, γ ) ν j ,K (γ ), (24)
where the W J M jkl,j k l are the atom symmetric top matrix elements in space fixed coordinates 20

W J M jKl,j K l (R, γ ) = λ μ V λμ (R, γ )X J M;λμ jKl,j K l , ( 25
) X J M;λμ jKl,j K l = δ J,J δ M,M (-1) [j +j +λ+J ] ×[(2j + 1)(2j + 1)(2l + 1)(2l + 1)(2λ + 1)] 1 2 × l j J j l λ l λ l 0 0 0 j λ j -K -μ K . ( 26
)
The Close Coupling equations to solve for an atom-triatomic rigid bender collision and for given values of J and M are then

d 2 dR 2 - l(l + 1) R 2 + k 2 νj (E) -[V A-RB ] J M νj l,ν j l (R) × G J M νj l,ν j l (R) = 0, ( 27
)
where k 2 νj (E) = 2μ H e-H CN [Eε nj ] and ε νj are the rigid bender energies calculated previously. From this expression we see that the atom-triatom rigid bender equations are similar in form to those of the atom linear molecule system, the only differences being the matrix elements of the potential and the rigid bender energies. Indeed if we put k = k = 0 and fix γ = 180 • we obtain the usual atom-linear molecule Close Coupling equations in space fixed coordinates. 21 The parity of the triatomic molecule can be easily introduced in these equations using instead of (20):

|νj lpJ M = mj ml 2J + 1 4π jm j lm l |J M Y ml l ( R) × K≥0 ˜ νp j, K (γ ) 2(1 + δ K0 ) D j * mj , K φ RHCN , θ RHCN , 0 + (-) p+ K D j * mj ,- K (φ RHCN , θ RHCN , 0) . ( 28
)
The global parity of this basis set is ε = (-) (j+l+p) while it is (-) (j+l+p+K) for a atom-symmetric top collision and (-) (j+l) for an atom colliding with a linear molecule. The calculations can then be split into two non-interacting parity blocks for each value of J. The matrix elements of the potential in this basis set are

[V A-RB ] J M νjpl,ν j ql (R) = νj lpJ M|V A-RB |ν j l qJ M = π 0 dγ sin γ K K ˜ νp j, K (γ ) W J M j Kpl,j K ql (R, γ ) ˜ ν q j , K (γ ) (29)
with

W J M j Kpl,j K ql (R, γ ) = λ μ
V λμ (R, γ ) XJM;λμ j Kpl,j K ql (30) and following Green: 20 XJM;λμ

j Kpl,j K ql = δ ε,ε [(1 + δ K0 )(1 + δ K 0 )] 1 2 × ωX J M;λ±μ j Kl,j K l + (-) p+ K δ μ,( K + K) X J M;λμ j -Kl,j K l , ( 31
)
where

ω = 1 if μ = K -K ≥ 0 (-) μ if μ = K -K ≥ 0 .
We have deliberately used the parity index ε in this expression in order to show that the symmetrised matrix elements of the potential are zero when the parity of the initial and final states is different. This expression shows that the results of the RBAA approach which uses the definition of parity for linear HCN cannot be completely equivalent to those given by the rigid bender inelastic collisions at the close coupling (RB-CC) method for transitions between different bending levels as discussed below.

III. CALCULATIONS AND RESULTS

We first give in Table I the rigid bender energies calculated using the approach presented in Sec. II A. These energies are compared to the experimental values reported by Harris et al. 22 As can be seen in this table the Rigid Bender Approximation gives energies which are about 2% accurate. The l doubling splitting is however over estimated by this approach as its absolute value is in average approximately multiplied by a factor five. This approach would then not be valid for spectroscopic purpose but is expected to be accurate enough to model the He-HCN collisional dynamics.

In the following we compare the results obtained with the atom-linear molecule close coupling (ALM-CC), RBAA, and RB-CC approaches and using our potential. 7 In all calculations we use 20 rotational levels of HCN. For the RB-CC and RBAA calculations we include in the basis set the three first bending levels for each rotational level of HCN. For each value of K (the projection of the rotational angular momentum of HCN along the Z molecular fixed axis), we use 30 associated Legendre functions to calculate the bending levels using the equations presented in Sec. II. The integral over the bending angle γ necessary to obtain the matrix elements of the inter molecular potential (28) was performed using a Gauss Legendre quadrature of 40 values of the bending angle.

We first consider the pure rotational transitions taking place inside the fundamental bending level ν = 0. We represented in Figures 234the elastic and de-excitation cross sections respectively of the levels (ν = 0, j = 1, 2, 3) of HCN. For each transition represented, we can observe a group of two curves which are almost identical and are associated with the ALM-CC and RBAA approaches while the third curve associated with the present RB-CC results is slightly different. The differences between the cross sections obtained from the three types of calculations are in any case negligible in the [1:1000] cm -1 energy interval needed for Astrochemistry. Only the very low energy range and the regions of the resonances are slightly different. This demonstrate that for these transitions the ALM-CC approach offers a level of accuracy equivalent to the RBAA approach, as already concluded in Paper I, and is, for the He-HCN system, almost equivalent to exact calculations. TABLE I. HCN rigid bender energies in cm -1 relative to the fundamental level (ν = 0, j = 0). The experimental HCN energies 22 are also given in the parentheses for comparison. As mentioned in the Introduction, l-type transitions of HCN have been detected recently in hot molecular gas, for example, in the proto-planetary nebula CRL 618. The double degeneracy of the bending mode of HCN is indeed lifted when the molecule is bending and rotating simultaneously giving rise to l-type doubling for j ≥ 1. For the first excited bending mode and for j ≥ 1, every rotational level is split into two sub-levels. Hereinafter we designate these two sub-levels by ν = 1 and ν = 2 and examine the l-type transitions between these sub-levels. In Figs. 5 and7, the elastic and the deexcitation cross sections starting from the (ν = 2, j = 1) level of HCN respectively calculated using the RBAA and RB-CC approaches are represented as a function of the collision energy. The upper panels of each of these two figures show all FIG. 2. Comparison of the elastic and de-excitation cross section of HCN (ν = 0, j = 1) in collisions with He as a function of collision energy calculated using the RB-CC, RBAA and ALM-CC approaches. The final level is indicated by one integer designating the rotational quantum number. these transitions while the lower panels show a blow up including only the transitions which are not labeled in the upper panel. The only open rotational channel of de-excitation towards the bending level ν = 1 is j = 1 as the (ν = 2, j = 1) and (ν = 1, j = 1) levels are almost degenerate while the first excited bending energy of the level j = 0, which is exempt from l-type doubling, is quite higher. The open rotational channel of de-excitation towards the bending level ν = 0 are the levels j = 0-19, j = 19 being the highest value of j considered in our calculations. As it can be seen in these figures, the RBAA gives at the best the right order of magnitude for some of the transitions but is definitely not accurate enough even for astrochemical purposes. The elastic and l-type transition cross sections (towards the levels (ν = 2, j = 1) and (ν = 1, FIG. 3. Comparison of the elastic and de-excitation cross section of HCN (ν = 0, j = 2) in collisions with He as a function of collision energy calculated using the RB-CC, RBAA, and ALM-CC approaches. The final level is indicated by one integer designating the rotational quantum number. that the RBAA approach fails to give an accurate estimate of the magnitude of the cross sections but gives at least the right ranking of the transitions cross sections. Not surprisingly, we find with both methods that the j = 0 transitions are favored and that the magnitude of the cross section decreases when j increases. Another interesting feature can be seen when comparing Figs. 5 and7. A dip in the RBAA cross sections appear around 700 cm -1 on many curves while it is absent from the RB-CC curves. This can be understood when looking at Figs. 8 and9 where the state selected (ν = 2, j = 1 → ν = 0, j = 17) transition cross sections are represented respectively at the RBAA and RB-CC levels. Clearly the dip appears at the RBAA level while it is absent of the RB-CC cross section. On the same figure the cross section associated with the opening of the (ν = 1, j ) and (ν = 2, j ) channels are represented with j = 15, 16, and 17. These levels are the most strongly coupled with the channel (ν = 0, j = 17) as they both are linked by the j ≤ 2 rule. We can see in Fig. 8 that the minimum of the dip of the RBAA curve corresponds approximately to the opening of the (ν = 0, j = 15) channel while the RB-CC (ν = 2, j = 1 → ν = 0, j = 17) cross section (Fig. 9) is reduced but shows no dip. We conclude from these different comparisons that the RBAA approach overestimates the coupling be- tween the bending levels. This analysis is confirmed when comparing the RBAA and RB-CC (ν = 2, j = 1 → ν = 0, j = 1) transition cross sections which are also represented in Figs. 8 and9. The RB-CC cross sections are effectively lower than those calculated using the RBAA approach. This is not really surprising as we use parity selected matrix elements of the potential both in the RBAA and RB-CC approaches while the definition of the parity for linear molecules used for the RBAA calculations differs from the one valid for a rigid bender which is used in the RB-CC method as detailed in Sec. II.

ν = 0 ν = 1 ν = 2 j HCN ε E ε E ε E + 0.

IV. CONCLUSION

We presented a new method for calculating exactly rovibrational cross sections for collisions between an atom and a rigid bender triatomic molecule. This approach was applied to the He-HCN collision and its results were compared to those obtained when using the rigid bender averaged approximation. We find for this system that the RBAA approach is almost equivalent to exact calculations for pure rotational transitions taking place inside the fundamental bending level ν = 0. On the contrary, for transitions involving two different bending levels, the RBAA approach fails to give an accurate estimate of the magnitude of the cross sections but gives most of the time at least the right ranking of the transition cross sections. We then conclude that l-type transitions cross sections have to be calculated at the RB-CC level for the He-HCN collision while pure rotational transitions cross sections may be calculated accurately at the RBAA level. This result should hold for other triatomic molecules which have a similar bending frequency but needs to be tested on other systems.

63 cm -1 ). Transitions involving several vibrational and rovibrational levels of C 3 have been observed. In some regions of the ISM, the LTE is not a good approximation and the collisional rate coefficients of C 3 with H 2 and He are needed to model its abundance.

Zhang et al [99] analysed the surfaces of the C 3 -He complex computed with two methods: CCSD(T) and MP4 for the linear configuration of C 3 . Also they compared these surfaces with one computed at CCSD(T) level when C 3 is bend, but only for a single bending angle (γ = 160 • ). However, to our knowledge, the collisional rate coefficients available for the C 3 molecule are only those calculated by Abdallah et al [100] for its collisions with Helium, within the RMA. However, due to the opening of the vibrational bending channel at very low temperature, their validity is limited to a low temperature interval and a new study including the bending motion is then necessary. Furthermore, the floppy character of this molecule makes of the study of its collisions with He a peculiarly interesting test of the RB-CC method. The building of a new PES including the bending of C 3 , necessary to perform these calculations, is first presented in the present chapter. The modifications of the RB-CC method needed to take advantage of the symmetry of C 3 are also presented as well as the computation of the bound levels of C 3 -He complex. We report the cross section computed at RMA and RBCC approach for the lower rotational levels and for the ground and first excited bending modes.

Rigid bender C 3 molecule

We cannot use for C 3 the same rovibrational Hamiltonian than the one introduced in our study of the HCN molecule [89] as it does not allow to take into account the C atoms exchange symmetry. We use instead the symmetrical Hamiltonian developed by Carter et al [101]. In this case, the x-axis is the bisector of the bending angle γ, while the molecule is lying in the xz plane. This system of coordinates is represented in Fig. 6.1. The kinetic operator in internal coordinates used in this case is described in details by Handy [102]. The Hamiltonian takes the form where H V and H RV are

H = H V + H RV + V, (6.1)
H V (γ) = - 2 2 1 µ 1 R 2 1 + 2 µ 2 R 2 2 ∂ 2 ∂γ 2 + cot γ ∂ ∂γ + 2 2m B R 1 R 2 cos γ ∂ 2 ∂γ 2 + cot γ ∂ ∂γ + ∂ 2 ∂γ 2 + cot γ ∂ ∂γ cos γ , (6.2)
H RV = 1 8 cos 2 (γ/2) 1 µ 1 R 2 1 + 1 µ 2 R 2 2 + 2 m B R 1 R 2 Π 2 z + 1 8 sin 2 (γ/2) 1 µ 1 R 2 1 + 1 µ 2 R 2 2 - 2 m B R 1 R 2 Π 2 x + 1 8 1 µ 1 R 2 1 + 1 µ 2 R 2 2 + 2 cos γ m B R 1 R 2 Π 2 y - 1 4 sin γ 1 µ 1 R 2 1 - 1 µ 2 R 2 2 Π z Π x + Π x Π z + 2i 1 µ 1 R 2 1 - 1 µ 2 R 2 2 cot γ 2 + ∂ ∂γ Π y . (6.3)
R 1 and R 2 are the rigid interatomic distances. Π is the total angular momentum operator of the molecule relative to the body-fixed axes. µ 1 and µ 2 are the reduced mass associated to each bond length as defined in the previous chapter, while m B is the mass of the central atom.

The rigid bender potential of C 3 is computed at the CCSD(T) level for the CC distance fixed at its experimental equilibrium value and is fitted using seven Legendre polynomials.

The last two terms in expression 6.3 vanish in for the C 3 molecule as R 1 = R 2 and µ 1 = µ 2 . In this case, the matrix elements of the H V and H RV in a symmetrised basis set, such as the one used previously for HCN, are

j K M, q| P K l |H V | P K l |j KM, p = δ p,q f (p, q)δ K , K 2 2 m b R 2 × l(l + 1)δ l,l - K 2 δ l ,l±2i (2l > + 1)(2l < + 1)(l > -K)!(l < + K)! (l > + K)!(l < -K)! - 2 2m B R 2 δ l ,l±1 l 2 > (l > + K)(l > -K) (2l > + 1)(2l > -1) + Kδ l ,l±(2i+1) (2l > + 1)(2l < + 1)(l > -K)!(l < + K)! (l < -K)!(l < + K)! , (6.4)
and

j K M, q| P K l |H RV | P K l |j KM, p = δ p,q f (p, q) 1 4m B R 2 δ K ,K × δ l ,l j(j + 1) -K 2 + δ l ,l±1 1 2 j(j + 1) -K 2 B l>K +δ l ,l±2i j(j + 1) + 5K 2 (l < + K)! K(l < -K)! N K l,l +δ l ,l±(2i+1) j(j + 1) -7K 2 (l < + K)! K(l < -K)! N K l,l +δ K ,K±2 A ± jK 2 P K l |(1 + 1 2 cos θ)| P K l + 1 2 P K l |( 1 sin 2 θ + cos θ sin 2 θ )| P K l , (6.5)
where i = 0, 1, 2, 3..., l < = min(l, l ), l > = max(l, l ) and

A ± jK = (j ± K + 1)(j ± K + 2)(j ∓ K)(j ∓ k + 1), (6.6)
B l>K = (l > + K)(l > -K) (2l > + 1)(2l > -1) , (6.7)
and

N K l,l = (2l + 1)(2l + 1)(l -K)!(l -K)! (l + K)!(l + K)! . ( 6.8)
The resulting matrix is diagonalized and this gives the rovibrational energies and the associated wave functions expanded over the symmetrised basis set 5.3.

Symmetry considerations

The C 3 molecule is made of three identical bosonic carbon atoms. The total wave function of rigid bender C 3 has then to be symmetrical under the exchange of the two terminal atoms [38] and the rovibrational wave function as well. The exchange of the two terminal atoms can be obtained by a rotation of 180 • around the x-axis and

C 2 (x)|jn KM, p = (-1) j+ K+p |jn KM, p , (6.9)
where |jn KM, p was defined in 5.3. Then the sum of the three quantum numbers j + K + p has to be even and one level out of two cannot exist as a result of the exchange symmetry. The resulting possible levels are represented in Fig. 6.2.

Rigid bender close coupling calculations

The interaction PES of C 3 with He has been computed at CCSD(T) levels using an aug-cc-pVQZ basis set with bond functions. These ab-initio energies have been fitted to an analytical function using a least square procedure similar to the one used for He-HCN. This four dimensional PES was used to solve the radial close coupling equations presented in the previous chapter. The bound levels of the C 3 -He complex have been calculated using two different methods: first, by considering C 3 to be a rigid linear molecule, i.e. at the RMA level, while in the second case, C 3 was considered to be a rigid bender using the RB-CC approach. The well depth of the He-C 3 complex is about 30 cm -1 while the bending frequency of C 3 is twice larger. This means that the bound states of the complex involve only the fundamental bending level of C 3 and it implies that the coupling between bending and rotation is not expected to be very important for these calculations. This is indeed what is found and detailed in the joined publication. The same com- parison is done for the calculations of the inelastic cross sections in order to show the effect of the coupling between the bending and the rotation of C 3 on to the dynamics. The RMA calculation are performed using a rotational basis set of 20 functions for C 3 and by taking the rotational constant to be equal to its experimental value B = 0.4305 cm -1 [97]. In the rigid bender calculations we included 20 rotational and 6 vibrational levels of C 3 . The de-excitation and excitation cross section of C 3 by collisions with He have been computed for the lower rotational levels and for the ground and first excited bending modes.

The variation of the cross sections of the pure rotational transitions (ν = 0) from j = 2, 4 as a function of collision energy are presented in Fig. 6.3. The RMA and RB-CC cross sections are found to be in excellent agreement at low collision energy (E < 100 cm -1 ) while at higher energies, the cross section computed at RB-CC level decreases faster than the RMA one. As

E(cm -1 ) RMA 0,4 → 0,0 RMA 0,4 → 0,2 RBCC 0,4 → 0,0 RBCC 0,4 → 0,2 RBCC 0,4 → 1,1 RBCC 0,4 → 1,3 RBCC 0,4 → 1,5 RBCC 0,4 → 1,7 100 10 RMA 0,2 → 0,0 RMA 0,2 → 0,2 RBCC 0,2 → 0,0 RBCC 0,2 → 0,2 E(cm -1 )
Figure 6.3: Comparison of the elastic and inelastic cross section using the RMA and the RB-CC approach. The levels are indicated by ν 0 , j 0 → ν f , j f . can be seen in this figure, this difference is due to the opening of the first excited bending channels of C 3 around 66 cm -1 . Such an effect could not be seen for the HCN-He system as the bending frequency of HCN is more than ten times larger. We conclude that for C 3 , even if one is interested only in pure rotational cross sections, the coupling between rotation and bending has to be taken into account for energies larger than 100 cm -1 as the RMA approach overestimates its value by 20% at 200 cm -1 (transition j = 4 → j = 0). The cross sections involving two different vibrational levels are also presented in Fig. 6.4 and 6.5 for the initial states ν = 1, j = 3 and ν = 1, j = 1 respectively. Their magnitudes are comparable to those of the pure rotational transitions cross sections showing that for this molecule bending and rotation have to be treated on the same footing.

Publication (submitted article)

Article submitted to J. Chem. Phys. on January 20, 2014.

I. INTRODUCTION

Since their first detection in dense intermolecular clouds in the seventies the vibrational ground state and the low-energy bending mode in stars forming cores 10 . In the dense interstellar medium (ISM) the presence of C 3 was proven using its mid-and farinfrared vibrational transitions 7,9,11 . C 3 was also identified in the diffuse interstellar medium, first by Maier et al. 8 and later in several other studies 12,13 along different lines of sight.

The detection of C 3 in the ISM has motivated many experimental [14][15][16][17][18] and theoretical [19][20][21][22] studies dedicated to the spectroscopy of this molecule. Contradicting measurements and theoretical calculations have alternatively reported a linear or non-linear equilibrium geometry for C 3 . Van Orden and Saykally 23 reviewed all previous works and proposed a linear structure with a very flat bending potential. They also discussed the barrier to linearity which is observed when the anti-symmetrical stretching is excited. A very small barrier to linearity (0.3 cm -1 ) is observed in the semi-empirical potential energy surface (PES) obtained by Špirko et al. 21 from a fit to experimental data. In contrast, the high-level CCSD(T) calculations of Mladenovic et al 20 resulted in a linear equilibrium geometry. While we can not still conclude whether C 3 is a linear or a quasi-linear molecule (i.e. a molecule with a barrier to linearity much smaller than the bending frequency), a good quantitative agreement has been obtained with the experimental energy levels 20,21 .

When the local thermal equilibrium (LTE) conditions do not apply, the determination of the abundance of C 3 require a study of its (de-)excitation resulting from collisions. In molecular clouds, the most abundant collider is usually H 2 , followed by He. Such calculations are usually performed within the rigid monomer approximation (RMA). However, it must be noted that the bending frequency of C 3 is small 24 (63.4 cm -1 ), and therefore a significant coupling of the vibrational bending with the intermonomer vibrational modes can be expected. The small bending frequency of C 3 brings a supplementary interest to the investigation of collision with He, as it could be seen as a case study of rovibrational energy transfer induced by collision.

We have recently developed the theory of rovibrational energy transfer for atoms colliding with a linear triatomic molecule and applied it to the HCN-He systems 25,26 . In this previous study 26 , denoted hereafter as Paper I, the coupling between rotation and vibrational bending was treated exactly within the rigid bender approximation at the Close Coupling level (RB-CC). Using this approach, we found that for the HCN-He collision, the pure rotational transitions could be computed accurately using a vibrationally averaged PES, similarly as was done in the investigations of CH 2 -He 27 , or even with the simple RMA. But for transitions involving two different bending levels, the use of the RB-CC approach was needed. In order to apply the latter approach to the C 3 -He inelastic collision, we will first focus, in the present work, on the development of a four dimensional PES for the C 3 -He system, which takes into account the bending motion. A first use of this PES to determine the bound states energies of the C 3 -He complex will also be presented.

The paper is organised as follows. In Section 2, we summarise the ab-initio calculations and detail the analytical form of the PES while the bound states calculations are presented in Section 3. Finally, we present and discuss the results of this study in section 4. Within the supermolecular approach, the potential energies of C 3 with He have been calculated with the coupled-cluster method with single and double excitations and a perturbative treatment of triple excitations (CCSD(T)). The interaction energy was corrected at all geometries for the basis set superposition error (BSSE) with the counterpoise procedure of Boys and Bernardi 28 . A comparison of the interaction energies calculated with basis sets of triple, quadruple, and quintuple-zeta quality is shown in Table I, with or without an additional set of bond functions centred at mid-distance between the He atom and the C 3 center of mass. The interaction energy, calculated at a configuration close to the equilibrium geometry, is quite stable in respect of the size of the basis set and the use of bond functions.

II. AB-INITIO CALCULATIONS AND ANALYTICAL FORM OF

For the largest basis set, it is safe to assume that the convergence of the one-electron basis is close to the complete basis set limit. Considering the computational cost associated with the various basis sets, we have chosen the quadruple zeta basis set with bond functions.

The potential energy has been computed for 12 values of the bending angle from 180 • to 80 The potential energy of the isolated C 3 molecule has been computed using the same abinitio method and basis set as we used for the C 3 -He system. We computed the ab-initio energy for 43 bending angles in the range from 50 • to 180 • . These energies were fitted to a linear combination of seven Legendre polynomials.

B. Analytical representation of the PES

The interaction potential is expanded as the sum of a short-range function V S and a long-range function V L with the expansion coefficients defined by the switching function S.

V int (R, θ, ϕ, γ) = S(R)V S (R, θ, ϕ, γ) + (1 -S(R))V L (R, θ, ϕ, γ) (1)
S(R) = 1 2 [1 -tanh (A 0 (R -R 0 ))] (2)
The V S and V L functions are in turn the sum of two terms. The first term represents the interaction energy in the case γ = π. Thus, it has no dependence in γ and ϕ. The second term is different from zero only for γ < π. It can be seen as the difference between the interaction energy at some γ < π and the interaction energy at γ = π, both at the same value of θ and R.

V

S (R, θ, ϕ, γ) = 12 l=0 F S l (R) Pl (cos θ) + 13 l=0 min (l,4) m=0 G S lm (R, γ) Plm (cos θ) cos (mϕ) (3) V L (R, θ, ϕ, γ) = 4 l=0 F L l (R) Pl (cos θ) + 4 l=0 min (l,2) m=0 G L lm (R, γ) Plm (cos θ) cos (mϕ) (4)
Due to the symmetry properties of the interaction potential (vide supra), the sum over l in the first term of right-hand side of eqs. 3 and 4 is restricted to even value and in the second term, the summation over l and m is restricted to even value of lm. Pl are normalized Legendre polynomials and Plm are normalized associated Legendre polynomials.

F S l (R) = e -α 1 R 11 n=0 R n A ln (5)
F L l (R) = 8 k=6 t k (β 1 R) R k C lk (6)
G S lm (R, γ) = e -α 2 R 5 j=0 Q m j (γ) 10 n=0 B lmnj R n (7) G L lm (R, γ) = 3 j=0 Q m j (γ) 8 k=6 D lmkj t k (β 2 R) R k (8)
where t k is the Tang-Toennies damping function:

t k (x) = 1 -e -x k i=0 x i i! (9)
The functions Q m j (γ) were chosen for all values of m, with the exception of m = 1, as

Q m j (γ) = 1 + P 2j+1 (cos γ) 2 (10)
and for m = 1,

Q 1 j (γ) = sin(jγ) (11)
In eq. 10, we used non-normalized odd Legendre polynomials.

The ab-initio points were fitted to this function in four steps, and in each step the least square procedure was used. In a first step, we fitted the linear case (γ = π), using the first term in equations 3 and 4. In a second step, the difference between the linear potential and the ab-initio points was fitted for each value of γ and R using the second term of 3 and 4.

These angular coefficients were interpolated using the functions 10 and 11. Finally we fitted the radial part to get the coefficients B lmnj and D lmkj .

The non-linear parameters R 0 , A, α 1 , β 1 , α 2 , and β 2 were set to the values R 0 = 5.2 Å,

A 0 = 5.6 Å-1 , α 1 = 4.8 Å-1 , β 1 = 2.20 Å-1 , α 2 = 2.
2 Å-1 and β 2 = 2.8 Å-1 . These values were determined by the trial and error method.

The ab-initio grid was computed only for γ ≥ 80 • and the extrapolation was done just until γ = 75 • . The rigid bender approximation used for C 3 is expected to be reliable only for the ground state and the first excited bending states. Considering that the potential energy of the C 3 molecule at γ = 75 • is about 5400 cm -1 and the vibrational frequency of the bending mode ν 2 is around 63 cm -1 , there is no need to represent the interaction energy for γ ≤ 75 • . Therefore we used this value as a cut-off limit. Beyond this limit, the interaction energy was set equal to its value at γ = 75 • . A similar procedure has been previously used 25 in the study of HCN in collision with He without any drawback.

III. BOUND STATES CALCULATIONS

Two kinds of calculations were performed, the first one considering C 3 as a linear rigid molecule i.e. using the RMA and the second one modelling C 3 by a rigid bender i.e. using the RB-CC method. In both cases we used the coupled-channel bound state method introduced long ago by Johnson 30 and adapted to the log-derivative and R-matrix propagators respectively by Huston 31 and Danby 32 . This method was already used in our study dedicated to the HCN-He system 25 . In the RMA approach, the rotational constant of the C 3

H = H v + H rv + V (γ) (12)
The vibrational term H v developed by Carter and Handy 34 is unchanged and its expression can be found in our previous work 25 . The rigid bender approximation of the rovibrational term H rv takes the following form

H rv = 3 4mR 2 1 1 cos 2 (γ/2) Π 2 z + 1 4mR 2 1 1 sin 2 (γ/2) Π 2 x + 1 2mR 2 1 1 + cos γ 2 Π 2 y ( 13
)
where R 1 is the C-C equilibrium distance and m is the mass of the carbon atom while Π x , Π y and Π z are the projections over the molecule-fixed axes of the total angular momentum

where M and K are the projections of j over the space-fixed z-axis and the body-fixed z-axis respectively. K is the absolute value of K, and p is equal to 0 or 1. Furthermore we have

|jn KM = P K n (cos γ)|jKM (15)
where

|jKM = 2j + 1 4π D j * M,K (α, β, γ) (16)
is a symmetric top wave function while P K n (cos γ) is a normalized associated Legendre polynomial describing the bending vibration. From the expression of H rv we see that a value of K can only be coupled with K = K, K±2. In other words, at each rovibrational state is associated a given parity of K. The calculation of the rovibrational wavefunctions are performed for a given value of j and p which means for a given parity as Π|jn KM, p = (-1) j+p |jn KM, p

We can further take advantage of the symmetry of the total wave function under the exchange of the bosonic C nuclei 35 which implies that the rovibrational wave function has to be symmetric too. The exchange of the carbon nuclei, as shown in Fig. 1, can be obtained by a rotation along the x-axis or equivalently by an inversion of the coordinates of the carbon nuclei followed by a rotation along the z-axis, C z 2 Π|j KM, p = (-) j+p+ K |j KM, p

This means that j + p + K has to be even. As each rovibrational state is associated with a given parity of K, we get immediately that for non-zero value of j one rovibrational state in two is missing, as illustrated for example by Gendriesch et al 18 .

Finally, the eigenvalues and eigenfunctions of H have been used in the calculations of the bound states of the He-C 3 complex with the RB-CC method. All the details of the RB-CC close coupling equations can be found in paper I. We included in the calculations six bending eigenfunctions of C 3 and ten rotational eigenfunctions for each bending functions.

The propagator step size was the same as for the RMA calculations.

IV. RESULTS AND DISCUSSION

A. Features of the potential energy surface

The quality of the fit was checked by using the root means square (RMS) of the differences between the ab initio and the interpolated energies. For this analysis, we have divided the grid of points in two regions. As we focus on the collision dynamics at low temperature, we can define the limit between the two regions by γ = 120 • . For such a bending angle, the energy of the isolated C 3 molecule is around 300 cm -1 . Thus the region defined by γ < 120 • may be represented with less accuracy. In the region γ ≥ 120 • , for the negative interaction energies E the RMS is 0.09 cm -1 , for 0 cm -1 ≤ E ≤ 100 cm -1 the RMS is 0.77 cm -1 and for 100 cm -1 ≤ E ≤ 5000 cm -1 the RMS is equal to 13.17 cm -1 . The second part of the grid (γ < 120 • ) has a RMS of 0.31, 2.68 and 35.46 cm -1 respectively. We have also checked the quality of the fit of the bending potential of the isolated C 3 molecule. The RMS for E < 500 cm -1 is 1.89 cm -1 , and for higher energies the relative error is less than 5.4%.

We found that the isolated C 3 molecule has a linear structure at equilibrium, with a flat potential energy which extends from γ = 180 • until about 160 • .

Figs. 2 and3 show contour plots of the interaction energy for selected values of γ with ϕ = 180 • and ϕ = 0 • respectively. For a linear configuration of C 3 , the minimum interaction energy with He is -26.73 cm -1 at R = 6.82 a 0 and θ = 81.4 • . By symmetry, there is a second minimum at θ = 98.6 • . The barrier to the θ = 90 • structure is 0.07cm -1 . This double minimum structure disappears when γ decreases as it can be observed in Fig. 2 for ϕ = 180 • . Conversely, the double minimum is more pronounced for ϕ = 0 • as seen in Fig. 3.

The potential well is 0.9 cm -1 deeper than the previous PES published by Abdallah et al 36 .

This small difference can arise from the use of a different basis set and a different distance between the carbon atoms. For R close to the equilibrium value, we can see in Figs. 2 and 3 that the interaction potential is strongly anisotropic. This is the consequence of the repulsion between He and the terminal carbon atoms, while the equilibrium structure is almost T-shaped. By comparing the different panels of Fig. 2, we observe that the long range energy is hardly changed while the short range energy varies strongly with γ. Consequently, the minimum of the interaction energy varies significantly with γ. Fig. 4 shows that the minimum of the interaction energy occurs when C 3 is bent. The global minimum of the interaction potential is -36,3 cm having the values found for the global minimum.

B. Bound states

In Table II, we report some of the lowest rovibrational energies obtained for the rigid bender C 3 molecule by diagonalising the Hamiltonian H (eq. 12), along with the corresponding vibrational quantum number ν, the total angular quantum number j and its projection K, and the value of the parameter p. As can be seen in this table, the rigid bender energies differ from the experimental values by less than 5%. This qualitative agreement indicates that the rigid bender model is a reliable approximation, at least for the lowest states.

The bound state energies of the He-C 3 complex, calculated with the RMA and RB-CC methods, are presented in Table III. The total angular momentum J and the parity ε are also reported. The energies computed using the RMA are found to be above those obtained using the RB-CC with only one exception at J = 4. This general behaviour can be explained by two factors. First the RMA rotational constant used for C 3 differs from the one obtained from the rigid bender calculations. Secondly, the minimum of the He-C 3 interaction potential is not in the linear configuration of C 3 as seen in Fig. 4. For J = 4, with the RB-CC method, we also found two more bound levels than with the RMA method.

The maximum value of the total angular momentum J leading to bound states is 6 using both methods. Since the ν = 0 → ν = 1 excitation energy of C 3 is larger than the well depth of the interaction potential, all the He-C 3 bound levels correspond to the bending quantum number of C 3 ν = 0. The dissociation energy computed using the RMA is 9.56 cm -1 while it is 9.73 cm -1 with the RB-CC method. The potential well supports 25 bound levels if C 3 is in the linear rigid configuration, and 27 is the bending motion is allowed. The relatively small differences between the results obtained with the RMA and RB-CC methods suggest that the RMA approach is a good approximation to compute the bound levels of an atom-triatom van der Waals complex, even for a very floppy molecule like C 3 .

V. CONCLUSION

We presented the first four dimensional analytical representation of the PES for the C 3 -He complex, including the bending motion of the C 3 molecule. This surface is based on supermolecular ab initio calculations using a quadruple zeta basis set with mid-bond functions and BSSE correction. The interaction potential has a minimum of -36,30 cm -1 at γ = 120 • , R = 6.57 a 0 , θ = 90 • and ϕ = 180 • . The bound states calculations for the C 3 -He complex were performed both at the RMA and at the RB-CC levels. The dissociation energy computed using the RMA is 9.56 cm -1 while it is 9.73 cm -1 at the RB-CC level.

The agreement between these two values shows that the RMA approach is valid for this system, mainly because the interaction potential is too shallow to couple excited bending states of C 3 , and also because the coupling between the inter-and intra-monomer motions is weak. More bound levels of lower energy were found however at the RB-CC level as a consequence of the non-linear geometry of C 3 associated with the minimum of the interaction potential. Future work dedicated to the rovibrational energy transfer in He-C 3 collisions will be presented shortly but we can however already expect that pure rotational transitions should be calculated correctly at the RMA level.

CHAPTER 7

General Conclusions

Conclusions

In this monograph the theoretical studies of four inelastic collisions of astrochemical interest were presented. The two first dedicated to the study of the collisions between two linear molecule were presented in the first part of this manucript while the two others were presented in the second part dealing with collisions between an atom and a rigid bender triatomic molecule. The conclusions reached in these two parts will now be reminded and the possible extensions of this work will also be briefly mentioned.

Inelastic collisions between two rigid linear molecules

In a first part, two studies of collisions of H 2 with a rigid linear molecules were respectively presented for the H 2 -CS and H 2 -HCN collisions. For these two systems a four dimensional PES was first developed. For the H 2 -CS system, the ab initio calculations were carried out at the coupled-cluster level with single and double excitations and a perturbative treatment of triple excitations, using a quadruple-zeta basis set and midbond functions, along with BSSE correction. This grid of points was fitted to an analytical function describing accurately both the region of the well and the long range asymptotic behaviour. The well depth of this PES is -173 cm -1 , and the equilibrium configuration is linear with the carbon atom pointing toward H 2 . All the rovibrational bound states of the H 2 -CS complex have been calculated. The dissociation energies of the ground states of the para and ortho species were found to be only 35.87 cm -1 and 49.88 cm -1 , respectively, indicating the energies of both ground states are a major part of the binding potential energy. We determined the cross sections for the first 16 rotational levels of CS, and the respective rate coefficients. In the case of the para-H 2 -CS system, a simple examination of some levels spacing indicate that the CS monomer behaves like a weakly hindered rotor. This is a consequence of the large zero point energy which allows the para-H 2 -CS system to move freely in almost the whole angular coordinates space. A set of rigid rotor close coupling cross sections for the inelastic collisions of CS with H 2 was also computed using this new 4D PES. For the collision of CS with para-H 2 , we observed a propensity to favor odd j CS over even j CS at low collision energy while at energy close to 1000 cm -1 , the propensity rule is reversed. The same tendency is observed for collisions involving ortho-H 2 . However, the inversion takes place at higher collision energy. This effect is related to the opening of the second channel of H 2 (j H 2 = 2 for para-H 2 or j H 2 = 3 for ortho-H 2 ). The cross sections for the rotational transition of CS in collision with para-H 2 were also compared with the scaled crosssections available for the collision of CS with He. The usual square root of the relative mass of the colliders, which is used by astronomers to obtain the cross sections between a molecule and H 2 from the data available with 4 He, is found to be a good qualitative approximation for this system.

The second study was dedicated to the HCN-H 2 collision considering HCN to be a rigid linear molecule. The 4D analytical PES developed for this system was based on a large grid of ab initio points computed at the CCSD(T)-F12a level and using an aug-cc-pVTZ basis set. The equilibrium structure of the HCN-H 2 complex was found to be linear with the nitrogen pointing towards H 2 . The corresponding well depth is 195.20 cm -1 . A secondary minimum was found only 11.61 cm -1 above the global minimum in which the H atom of HCN is pointing towards the center of mass of H 2 . As a first application, the rovibrational bound states were computed within the rigid-rotor approximation. The total number of bound states supported by our PES is 101 for HCN-para-H 2 and 330 for HCN-ortho-H 2 . The dissociation energies of the ground states of the para and ortho species are 37.79 cm -1 and 60.26 cm -1 , respectively. The calculated transitions frequencies are found to be in very good agreement with the experimental available data. This level of agreement suggests that our PES is accurate enough for computing accurate inelastic cross sections.

Inelastic collisions between an atom and a rigid bender triatomic molecule

The second part of this manuscript is dedicated to the development of new approaches allowing to include the bending-rotation interaction in the treatment of the inelastic collisions of a triatomic molecule with an atom. The two molecules considered are HCN and C 3 which both have been detected in several regions of the ISM and the latter is an archetypal example of a floppy molecule. These two studies required first the construction of an analytical model of PES including the bending vibration of the triatomic molecule. We started with the He-HCN system by building a four dimensional analytical representation of the PES based on supermolecular ab initio calculations using a quadruple zeta basis set with mid-bond functions and BSSE correction. The van der Waals well was found to be 30.35 cm -1 deep and associated with the linear configuration (He-HCN) while a secondary minimum of 22.08 cm -1 was also identified for a bent configuration. We checked first that the restriction of this PES to the rigid linear configuration of HCN gives similar close coupling inelastic cross section than the previous theoretical works and that the bound states calculated using this restricted PES are also in good agreement with the available experimental data. We also presented a simple method (RBAA) of calculation of the rotational close coupling cross section which uses the average of the interaction potential over the bending wave functions of HCN. We found that taking into account the bending motion through the RBAA method does not change significantly the rotational excitation cross sections, while the agreement of the calculated bound state transition frequencies with the experiment is marginally improved. This first study shows in any case that the RMA approach is quite satisfactory for the computation of rotational excitation cross sections for a linear triatomic molecule like HCN in its fundamental bending level. In a second publication we presented a new method for calculating exactly rovibrational cross sections for collisions between an atom and a rigid bender triatomic molecule. The results of this approach called RB-CC were compared to those obtained when using the RBAA approach.

ISM. As we have seen that l-doubling splitting is overestimated when using the rigid bender approximation, we will use the experimental values of the rovibrational energies of HCN which are available while using the rigid bender wave functions.

4. The calculation of the effect of dynamical switching, which is expected to be small for the lowest bending states of a linear molecule like HCN and was not taken into account in this first study, will also be investigated at higher collision energy for this system.

5. While there is no dynamical switching for the He-C 3 collision as we use a different rigid bender Hamiltonian for C 3 , the test of the use of experimental rovibrational energies of C 3 on the dynamics will also be of interest and the first determination of rovibrational rate coefficients needed by the astronomers will be performed shortly.

6. As the RB-CC method applies to any triatomic molecule whether linear or not, one can also think about many different use of the RB-CC method for collisions of He with other triatomic molecules or to model the collisions of these molecules by para or ortho-H 2 using available PES averaged over the H 2 rotational state. However as we know that this approximation sometime fails for the collision between H 2 and a rigid linear molecule, the extension of the RB-CC method to the case of a rigid linear molecule like H 2 interacting with a rigid bender triatomic molecule will be considered in a future work.

  Le premier chapitre introductif situe le contexte astrochimique des problèmes traités. Cette thèse est en effet financée par un programme européen d'astronomie: le programme Astronet CATS. Il s'agit donc de traiter les collisions inélastiques entre les molécules détectées et H 2 ou bien He qui sont les constituants les plus abondants des nuages interstellaires. Les modèles de chimie des nuages interstellaires se sont considérablement améliorés ces dernières années et un certain nombre de molécules triatomiques commencent à être prise en compte. Les systèmes traités dans cette monographie sont tout d'abord H 2 + CS et He + HCN linéaire. Ces deux systèmes, même s'ils nécessitent de calculer les surfaces d'énergie potentielle impliquées puis d'effectuer la dynamique, sont traités avec des méthodes existantes. En revanche la deuxième partie de la thèse est dédiée au développement de plusieurs méthodes originales permettant de décrire le couplage entre la rotation et le mouvement vibrationnel de pliage d'une molécule triatomique en collision avec un atome d'hélium. Ces modèles sont appliqués à deux molécules : HCN et C 3 . Ce chapitre introductif justifie donc la nécessité de prendre en compte le couplage entre rotation et pliage pour ces molécules dans le milieu interstellaire et mentionne les études existantes au moment du début de ce travail de thèse.

3. 1 4 . 1 Five coefficients ν m l 1

 1411 Set of coordinates used to describe the CS-H 2 system. The azimuthal angle ϕ is undefined when θ 1 or θ 2 is equal to 0 • or 180 • . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ,l 2 (R) of the potential energy expansion for the HCN-H 2 system . . . . . . . . . . . . . . . . . . . . . 5.1 Potential energy of the rigid bender model for the HCN molecule. 57 5.2 Coordinate systems for HCN-He. . . . . . . . . . . . . . . . 5.3 Rate coefficient for the collision of HCN with He for 4 ≤ T ≤ 100 K. The initial and final states are indicated as (ν 0 , j 0 ) → (ν f , j f ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Representation of the coordinates for the C 3 molecule. . . . 6.2 Representation of the allowed and not allowed states for the C 3 molecule. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Comparison of the elastic and inelastic cross section using the RMA and the RB-CC approach. The levels are indicated by ν 0 , j 0 → ν f , j f . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Elastic and inelastic RB-CC cross section of C 3 in collisions with He as a function of collision energy from the initial level ν = 1, j = 1. The final level is indicated by two integers designating the bending and the rotational quantum numbers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  are the spherical Hankel functions of the first and second kind.

Figure 3 . 1 :

 31 Figure 3.1: Set of coordinates used to describe the CS-H 2 system. The azimuthal angle ϕ is undefined when θ 1 or θ 2 is equal to 0 • or 180 • .

FIG. 1 .

 1 FIG.1. Set of body fixed coordinates used to describe the diatom-diatom system. The azimuthal angle ϕ is undefined when θ 1 or θ 2 is equal to 0 • or 180 • .

FIG. 2 .FIG. 3 .

 23 FIG. 2. Contour plot of the rigid rotor PES for θ 1 = 0 • . The contour lines are equally spaced by 10 cm -1 .

FIG. 4 . 60 FIG. 5 .

 4605 FIG. 4. Contour plot of the rigid rotor PES for θ 2 = 180 • . The contour lines are equally spaced by 10 cm -1 .

12 FIG. 6 .

 126 FIG.6. Cross second virial coefficient calculated with the present PES.

  FIG. 1. H 2 -CS interaction energy as a function of R for several attractive angular configurations. The values of θ 1 , θ 2 , and ϕ are indicated for each curve.

FIG. 2 .

 2 FIG.2. Rotational transition cross section of CS(j CS = 5) in collisions with para-H 2 (j H 2 = 0) as a function of the collision energy. The curves are labeled by (j H 2 , j CS ) → (j H 2 ,j CS ).

5 FIG. 5 .

 55 FIG.5. Rotational quenching cross section of CS(1 ≤ j CS ≤ 5) by collision with para-H 2 (j H 2 = 0) as a function of the collision energy. The initial rotational quantum number j CS is reported on each curve.

5 FIG. 6 .FIG. 7 .

 567 FIG.6. Rotational quenching cross section of CS(1 ≤ j CS ≤ 5) by collision with ortho-H 2 (j H 2 = 1) as a function of the collision energy. The initial rotational quantum number j CS is reported on each curve.

FIG. 8 .

 8 FIG.8. Comparison of the rotational quenching cross section of CS by collision with4 He (green lines) and with para-H 2 (j H 2 = 0) (red lines) as a function of the collision energy. The rotational transition j CS → j CS is reported in each panel.

1 FIG. 9 .

 19 FIG.9. De-excitation rate coefficients of CS(j CS = 15) by collisions with H 2 at 100 K as a function of j CS = j CSj CS .

Figure 4 . 1 :

 41 Figure 4.1: Five coefficients ν m l 1 ,l 2 (R) of the potential energy expansion for the HCN-H 2 system

FIG. 2 .

 2 FIG. 2. Potential energy cuts of 4D PES for several angular orientations.Energy is in cm -1 .

FIG. 3 .

 3 FIG. 3. Contour plot of the PES for θ = 0 • . Contour lines are equally spaced by 20 cm -1 . Red contour lines represent the repulsive interaction energies.

Figure 5 (

 5 FIG. 4. Contour plot of the PES for ϕ = 0 • and R relax in the range [7.0, 8.0] a 0 . Contour lines are equally spaced by 10 cm -1 . Red contour lines represent the repulsive interaction energies.

Figure 5 . 1 :

 51 Figure 5.1: Potential energy of the rigid bender model for the HCN molecule.

Figure 5 . 2 :

 52 Figure 5.2: Coordinate systems for HCN-He.

Figure 5 . 3 :

 53 Figure 5.3: Rate coefficient for the collision of HCN with He for 4 ≤ T ≤ 100 K. The initial and final states are indicated as (ν 0 , j 0 ) → (ν f , j f ).

FIG. 1 .

 1 FIG. 1. Definition of the body-fixed coordinate system for the He-HCN system. The planar configuration represented here corresponds to ϕ = 180 • . The angle ϕ is not defined for γ or θ equal to 0 • or 180 • .

FIG. 2 .

 2 FIG. 2. Contour plot of the PES for selected values of γ and for ϕ = 0 • . Negative contour lines (blue) are equally spaced by 3 cm -1 .

FIG. 3 .

 3 FIG. 3. Contour plot of the PES for selected values of γ and for ϕ = 180 • . Negative contour lines (blue) are equally spaced by 3 cm -1 .

FIG. 4 .

 4 FIG. 4. Interaction energy between He and l-HCN calculated with the present PES (solid lines) and with the S02 PES (dashed lines), for θ = 0 • (black), θ = 90 • (blue), and θ = 180 • (red).

FIG. 5 .

 5 FIG. 5. Comparison of the rotational transition cross sections of HCN in collisions withHe calculated using the RMA and the RBAA approaches for even j (top panel) and for odd j (bottom panel).

2 FIG. 1 .

 21 FIG. 1. Definition of the body-fixed coordinate system for the He-HCN system. The planar configuration represented here corresponds to ϕ = 180 • . The angle ϕ is not defined for γ or θ equal to 0 • or 180 • .

FIG. 4 .FIG. 6 .

 46 FIG.4. Comparison of the elastic and de-excitation cross section of HCN (ν = 0, j = 3) in collisions with He as a function of collision energy calculated using the RB-CC, RBAA, and ALM-CC approaches. The final level is indicated by one integer designating the rotational quantum number.

FIG. 7 .

 7 FIG.7. Elastic and de-excitation RB-CC cross section of HCN(ν = 2, j = 1) in collisions with He as a function of collision energy. The final level is indicated by two integers designating the bending and the rotational quantum numbers.

FIG. 9 .

 9 FIG. 9. Comparison of the inelastic RB-CC cross section of HCN(ν = 2, j = 1 → ν , j = 17) in collisions with He as a function of collision energy. The final level is indicated by two integers designating the bending and the rotational quantum numbers. The cross section associated with the opening of the first excited bending level of the rotational state j = 0 of HCN is also represented.

Figure 6 . 1 :

 61 Figure 6.1: Representation of the coordinates for the C 3 molecule.

Figure 6 . 2 :

 62 Figure 6.2: Representation of the allowed and not allowed states for the C 3 molecule.

Figure 6 . 4 :

 64 Figure 6.4: Elastic and inelastic RB-CC cross section of C 3 in collisions with He as a function of collision energy from the initial level ν = 1, j = 1. The final level is indicated by two integers designating the bending and the rotational quantum numbers.

Figure 6 . 5 :

 65 Figure 6.5: Elastic and inelastic RB-CC cross section of C 3 in collisions with He as a function of collision energy from the initial level ν = 1, j = 3. The final level is indicated by two integers designating the bending and the rotational quantum numbers.

Figure 1 .

 1 Figure 1. Body-fixed coordinates. The linear C 3 molecule is along the z-axis, and the bent C 3 molecule is in the plane xz. The angle ϕ is indefinite when γ = 180 • .

molecule 23

 23 was taken equal to 0.43 cm -1 . A modified version of our scattering code based on the log-derivative propagator, was used following the recommendations of Hutson 31 . The calculations were performed for two values of the propagator step size (0.05 a 0 and 0.1 a 0 ), and the values of the bound state energies were obtained from a Richardson extrapolation31 .We tested the convergence of the bound state energies of the He-C 3 complex as a function of the size of the rotational basis set of C 3 . As the nuclear spin of the carbon atom is zero, ten even values of the rotational quantum number were included in the basis set describing C 3 and the maximum propagation distance was set to 50 a 0 .The RB-CC approach presented in paper I had to be adapted to take into account the exchange symmetry of the C atoms inside C 3 as the rigid bender Hamiltonian used in this previous work cannot be used for symmetric triatomics. We use instead the symmetric Hamiltonian developed by Carter et al33 . Fig.1shows the coordinates used in this study for the C 3 molecule. The x-axis is the bisector of the bending angle γ, while the molecule is lying in the xz plane.The symmetric form of the Hamiltonian developed by Carter et al33 is written as usual as a sum of a vibrational, a rovibrational and a potential term

j of C 3 .

 3 As in our previous work, we use the following symmetrised rovibrational basis set to diagonalise the Hamiltonian H |jn KM, p = |jn KM + (-1) p |jn -KM

Figure 2 .

 2 Figure 2. Contour plot of the interaction energy for selected values of γ and for ϕ = 180 • . Negative contour lines (blue) are equally spaced by 4 cm -1 . The lowest positive contour line (red) show the 4 cm -1 energy and there is a factor of 2 in energy between successive positive contour lines (red).

Figure 3 .

 3 Figure 3. Contour plot of the interaction energy for selected values of γ and for ϕ = 0 • . The contour lines have the same spacing as in Fig. 2.

1 Figure 4 .

 14 Figure 4. Interaction energy versus the bending angle of C 3 . All other coordinates R, θ and ϕ are relaxed.

Figure 5 .

 5 Fig 5 shows the contour plot of the interaction energy in function of θ and ϕ, with R and γ

TABLE II .

 II Calculated energies (cm -1 ) of the rovibrational bound states of oH 2 -CS for J ≤ 2, along with the associated quantum numbers (J, ε).

		State			State				State			State					State
	J	ε		Energy	J	ε	Energy	J	ε		Energy	J	ε	Energy		J	ε	Energy
	0	+		78.1927	1	+	105.7680	1	-		110.1782	2	+	102.9420	2	-	89.4881
	0	+		84.9060	1	+	108.5430	1	-		112.9009	2	+	104.1448	2	-	91.2534
	0	+		92.1872	1	+	110.7159	1	-		114.0169	2	+	106.1486	2	-	92.6863
	0	+		98.4315	1	+	113.4028	1	-		114.5219	2	+	106.8613	2	-	95.1564
	0	+		101.7227	1	+	113.8182	1	-		116.2404	2	+	107.3817	2	-	96.0826
	0	+		104.1509	1	+	115.1382	1	-		117.3911	2	+	109.7498	2	-	97.7320
	0	+		111.9908	1	+	117.7081	1	-		117.9314	2	+	110.4607	2	-	99.6486
	0	+		112.2215	1	+	118.9015	1	-		118.4700	2	+	112.2532	2	-	101.3505
	0	+		114.1382	1	+	120.3939	1	-		120.9221	2	+	113.1311	2	-	102.3588
	0	+		117.1608	1	+	120.7320	2	+		71.8257	2	+	113.5232	2	-	103.8251
	0	+		120.6376	1	-	73.9313	2	+		75.1621	2	+	114.3119	2	-	105.6296
	0	-		87.2179	1	-	78.9144	2	+		78.8875	2	+	115.2841	2	-	106.4514
	0	-		94.7724	1	-	82.4785	2	+		80.3857	2	+	115.7166	2	-	108.0734
	0	-		104.9616	1	-	84.2987	2	+		84.2431	2	+	117.4254	2	-	110.1003
	0	-		115.3423	1	-	87.1817	2	+		84.6973	2	+	117.7393	2	-	111.3348
	0	-		120.0817	1	-	89.8010	2	+		87.5395	2	+	118.3457	2	-	111.7665
	1	+		74.0346	1	-	91.7632	2	+		89.1242	2	+	119.0900	2	-	113.8275
	1	+		82.5703	1	-	95.2903	2	+		90.4075	2	+	119.5686	2	-	115.4673
	1	+		85.2856	1	-	98.5494	2	+		92.2985	2	+	120.4822	2	-	116.4560
	1	+		88.0804	1	-	99.2229	2	+		94.2299	2	+	120.9077	2	-	117.0772
	1	+		90.9907	1	-	99.7993	2	+		95.7337	2	+	121.1620	2	-	118.1388
	1	+		92.4167	1	-	100.9420	2	+		97.1789	2	-	75.4563	2	-	118.6224
	1	+		95.8193	1	-	104.9230	2	+		99.2674	2	-	78.9411	2	-	119.1144
	1	+		99.5506	1	-	105.6870	2	+		100.7164	2	-	84.5982	2	-	120.3551
	1	+		100.5843	1	-	108.0265	2	+		101.3734	2	-	85.5051	2	-	121.0159
	1	+		102.0414	1	-	109.8234	2	+		101.8032	2	-	88.2248
						.8716	1	-	3	2	8 8	-21.7297	2	+	1	1	5	-2.0243
	0	+	1	1	-32.0920	1	-	3	4	9 6	-15.9052			1	3	5
	0	+	2	2	-25.7775	1	-	4	3	9 6	-13.0722			2	0	6
	0	+	3	3	-19.2932	1	-	4	5	9 7	-5.6069			2	2	13
	0	+	4	4	-9.5563	1	-	2	1	8	-3.1860			2	4	36
	0	+	1	1	-3.6394			2	3	23				3	5	6
			2	2	74				5	4	58				4	6	7
			5	5			1	-	5	4	9 2	-1.9849			5	3	15
	0	+	0	0	-1.0338	1	-	1	0	1 2	-0.1433	2	+	1	3	5	-1.5577
			1	1	33				1	2	12				3	1	10
			2	2	16				2	3	10				4	6	70
			3	3	41				3	4	36		2	+	4	6	9 2	-0.5862
	1	+	1	1	-33.8144			4	5	6		2	-	1	2	9 6	-31.8691
	1	+	2	2	-27.6721			5	4	16		2	-	2	1	9 6	-30.5297
	1	+	3	3	-19.2149	2	+	1	1	8 0	-33.2079	2	-	2	3	8 9	-25.2144
	1	+	4	4	-9.8959	2	+	0	2	7 6	-33.0614	2	-	3	2	8 9	-22.7255
	1	+	1	1	-1.5300	2	+	2	0	8 5	-30.9901	2	-	3	4	9 6	-16.0049
			2	2			2	+	2	2	5 0	-28.6911	2	-	4	3	9 6	-13.1797
			3	3			2	+	1	3	4 5	-26.9339	2	-	4	5	9 8	-5.6317
			4	4			2	+	3	1	7 9	-23.8148	2	-	5	4	9 7	-2.1704
			5	5			2	+	2	4	7 8	-20.6995	2	-	2	3	7	-0.4070
	1	-	0	1	-34.9375	2	+	3	3	7 9	-19.8274			3	2	11
	1	-	1	0	-34.3308	2	+	4	2	9 4	-15.6287			3	4	36
	1	-	1	2	-30.9609	2	+	3	5	9 7	-11.7101			4	5	34
	1	-	2	1	-28.5378	2	+	4	4	9 8	-9.6357
	1	-	2	3	-23.8440	2	+	5	3	9 7	-5.3502

TABLE III .

 III Energy spacing (cm -1 ) between pH 2 -CS bound states associated with two successive values of j 2 for L = 0.

				Rotational spacing
			Energy spacing	of the free	Relative
	j 2	j 2	between levels	CS diatom	difference (%)
	1	0	1.54	1.64	6.1
	2	1	3.34	3.28	1.8
	3	2	4.65	4.92	5.5
	4	3	6.99	6.56	6.6
	5	4	8.26	8.20	0.8
	6	5	9.66	9.84	1.9

TABLE IV .

 IV Energy spacing (cm -1 ) between para bound levels associated with successive values of L. B rel is the rotational constant of the two-body system H 2 + CS calculated from the spacing.

				Energy spacing
	j 2	J = L	J = L	between levels	B rel
	0	1	0	0.93	0.47
		2	1	1.88	0.47
		3	2	2.52	0.42
		4	3	4.28	0.54
		5	4	4.98	0.50
		6	5	5.13	0.43
		7	6	6.34	0.45
		8	7	6.97	0.44
	Average				0.46
	1	2	1	1.95	0.49
		3	2	2.84	0.47
		4	3	4.12	0.51
		5	4	4.97	0.50
		6	5	5.43	0.45
		7	6	6.64	0.48
		8	7	7.22	0.45
	Average				0.48

rotator, i.e., follow a B rel L(L+1) law, then the value of B rel can be calculated for each energy spacing. Table

4 )

 4 FIG.4. Inelastic cross section for the opening of the first excited rotational levels of para-and ortho-H 2 as a function of the collision energy. The curves are labeled by (j H 2 , j CS ) → (j H 2 ,j CS ).

  2,b) 1 Université de Bordeaux, ISM, CNRS UMR 5255, 33405 Talence Cedex, France 2 Departamento de Física, Universidad de Matanzas, Matanzas 40100, Cuba 3 LOMC -UMR 6294, CNRS-Université du Havre, 25 rue Philippe Lebon, BP 540, 76058, Le Havre, France 4 Department of Optics and Spectroscopy, Tomsk State University, 36 Lenin av., Tomsk 634050, Russia 5 Instituto Superior de Tecnologías y Ciencias Aplicadas, Quinta de Los Molinos, Plaza, La Habana 10600, Cuba (Received 4 October 2013; accepted 12 November 2013; published online 9 December 2013)

TABLE I .

 I Molecular properties of monomers HCN and H 2 calculated at CCSD(T)/aV5Z level of theory using finite-field method. Literature values are in parenthesis. All values are in a.u.

	Param.	Definition	HCN a	H 2
	μ z	Dipole moment	-1.188 (-1.180 b )	0
	zz	Quadrupole moment	1.703 (1.646 b )	0.4823
	zzz	Octupole moment	-9.822 (-9.762 b )	0
	zzzz	Hexadecapole moment	23.274 (22.45 b )	0.3177
	zzzzz	5th-order multipole moment	-110.120	0
	zzzzzz	6th-order multipole moment	352.377	0.1702
	α xx	Dipole polarizability	13.933	4.7284
	α zz		22.248	6.7168
	A z, zz	Dipole-quadrupole polarizability -11.714	0
	A x, zx		-1.489	0
	E x, xxx	Dipole-octupole polarizability	-34.229	-1.7800
	E z, zzz		100.156	4.4462
	C xx, xx	Quadrupole polarizability	33.996	4.8331
	C xz, xz		39.886	4.4469
	C zz, zz		67.184	6.3610

TABLE III .

 III Lower bound levels (E < 77 cm -1 ) of the HCN-ortho-H 2 van der Waals complex. The approximate quantum numbers j 2 , j 12 and L are also given.

			State	Energy			State	Energy
	J	ε	(j 2 , j 12 , L)	(cm -1 )	J	ε	(j 2 , j 12 , L)	(cm -1 )
	0	+	(0,1,1)	58.38	2	-	(2,2,2)	72.97
	0	+	(1,2,2)	68.60	3	-	(0,1,2)	63.48
	0	+	(2,3,3)	75.26	3	-	(1,2,1)	67.74
	1	-	(0,1,2)	59.24	3	-	(1,2,5)	69.42
	1	-	(1,1,1)	63.95	3	-	(2,1,2)	74.41
	1	-	(1,2,1)	65.30	3	-	(2,3,4)	76.65
	1	-	(2,2,2)	70.11	3	+	(1,1,2)	67.91
	1	-	(2,3,2)	74.00	3	+	(1,1,4)	70.33
	1	+	(1,1,0)	63.93	3	+	(2,2,1)	75.85
	1	+	(1,1,2)	65.63	4	+	(0,1,3)	66.83
	1	+	(2,3,3)	76.52	4	+	(1,2,2)	70.61
	2	+	(0,1,1)	60.94	4	+	(1,2,6)	73.03
	2	+	(1,2,0)	65.56	4	-	(1,2,3)	70.94
	2	+	(1,1,2)	66.76	4	-	(2,2,2)	74.22
	2	+	(2,2,1)	72.46	5	-	(0,1,4)	70.97
	2	+	(3,4,2)	73.03	5	-	(1,2,3)	74.25
	2	+	(3,4,4)	74.93	5	+	(1,1,4)	74.69
	2	-	(1,1,1)	65.58	6	+	(0,1,5)	75.88
	2	-	(1,1,3)	67.45
									-	(1,0)	-27.89
								0	+	(1,1)	-31.46	2	-	(2,1)	-24.03
								0	+	(2,2)	-23.87	3	-	(0,3)	-33.16
								1	-	(0,1)	-37.00	3	-	(1,2)	-26.90
								1	-	(1,0)	-31.27	3	-	(1,4)	-24.63
								1	-	(1,2)	-29.09	3	-	(2,1)	-21.92
								1	-	(2,1)	-25.39	3	+	(1,3)	-24.69
								1	-	(2,3)	-20.68	3	+	(2,2)	-21.81
								1	+	(1,1)	-30.08	4	+	(0,4)	-30.14
								1	+	(2,2)	-23.28	4	+	(1,3)	-22.95
								2	+	(0,2)	-35.45	4	+	(1,5)	-22.19
								2	+	(1,1)	-29.69	4	-	(1,4)	-20.63
								2	+	(1,3)	-26.30	5	+	(0,5)	-26.40
								2	+	(2,0)	-25.51	6	+	(0,6)	-21.97
								2	+	(2,2)	-22.08

TABLE IV .

 IV Comparison of observed and calculated transition frequencies in MHz.

	Transition	Frequencies
	J -J	Calculated	Observed	% difference
	HCN-para-H 2
	1-0	23 649	23 665 a	0.07
	3-2	68 823	69 155 b	0.48
	4-3	90 616	90 933 b	0.35
	5-4	112 040	112 149 b	0.10
	6-5	133 024	133 575 b	0.41
	HCN-ortho-H 2
	1-0	25 723	25 768 a	0.18
	3-2	76 103	76 209 b	0.14
	4-3	100 446	100 590 b	0.14
	5-4	124 179	124 370 b	0.15
	6-5	147 368	147 614 b	0.17

a Average of hyperfine components from Ref. 23. b Average of hyperfine components from Ref. 20.

1

  Université de Bordeaux, ISM, CNRS UMR 5255, 33405 Talence Cedex, France 2 Departamento de Física, Universidad de Matanzas, Matanzas 40100, Cuba 3 Université Pierre et Marie Curie, LPMAA, UMR CNRS 7092, 75252 Paris, France 4 Observatoire de Paris, LUTH, UMR CNRS 8102, 92195 Meudon, France (Received 13 April 2013; accepted 24 June 2013; published online 16 July 2013)

TABLE II .

 II Bound levels of the He-HCN van der Waals complex.

TABLE III .

 III Comparison of observed and calculated transition frequencies in MHz.

			RMA		RBAA
	Transition	Observed	Calculated	% error	Calculated	% error
	j = 1 ← 0
	P(1)	97034 a	97696	-0.7	97198	0,2
	P(2)	96756 a	98411	-1.7	97823	1,1
	P(3)	98149 a	100188	-2.1	99477	1,4
	P(4)	101559 a	103801	-2.2	102900	1,3
	Q(1)	98132 b	101236	-3.2	100534	2,4
	Q(2)	101191 b	104373	-3.1	103545	2,3
	Q(3)	106244 b	109482	-3.0	108453	2,1
	Q(4)	113737 b	116863	-2.7	115565	1,6
	R(0)	98696 b	101006	-2.3	100328	1,7
	R(1)	101432 c	103782	-2.3	102962	1,5
	R(2)	105795 c	108350	-2.4	107295	1,4
	R(3)	112782 a	115460	-2.4	114074	1,1
	R(4)	122944 a	. . .	. . .	. . .	. . .
	j = 0 ← 0
	R(0)	15894 c	15674	1.4	15669	1,4
	R(1)	31325 c	30895	1.4	30892	1,4
	j = 1 ← 1
	R(2)	4604 c	3976	13,6	3750	18,6

a Reference 21. b Average of hyperfine components from Ref. 21. c Reference 18.

  1 , the carbon chains are expected to play an important role in the chemistry of interstellar and circumstellar clouds. It is especially interesting to estimate the abundance of the small pure carbon molecules, such as C 2 , C 3 and C 5 , because they are supposed to be the building blocks for other more complex interstellar molecules. Pure carbon chains have no permanent dipole moment but can be detected in dense and cold interstellar clouds through their infrared active low-energy bending vibrations as suggested byVan Orden 2 . In the present study we focus our interest on the carbon trimer C 3 . Its emission spectrum near 4050 Å was observed3,4 first in comets as early as 1881, and it was identified 5 in the laboratory in 1951.Since then, C 3 has been observed in several regions of the interstellar medium (ISM): in the atmospheres of cool stars 6 , in circumstellar shell7 and in diffuse interstellar clouds 8 . Recently, Cernicharo et al. 9 have detected nine lines of the bending mode of C 3 in Sagittarius B2 and IRC +10216. Mookerjea et al. observed several rovibrational transitions between

Table I .

 I CCSD(T) interaction energy of the C 3 -He system at R = 6.8 a 0 , θ = 90 • and γ = 180 • . The use of bond functions is denoted by +bf. were spaced uniformly in steps of 10 • for θ in the range [0,90] • and steps of 30 • for ϕ in the range [0,180] • . The total number of points was 26526. All calculations have been carried out with the MOLPRO package 29 .

	Basis	Energy (cm -1 )	relative time cost
	aug-cc-pVTZ+bf	-26,53	1
	aug-cc-pVQZ	-24,62	2,7
	aug-cc-pVQZ+bf	-26,64	4,8
	aug-cc-pV5Z	-25,62	16,5
	aug-cc-pV5Z+bf	-26,72	28,3

• . The radial grid included 34 points ranging from 4.1 a 0 to 22.7 a 0 . The angular grids

Table II .

 II The rovibrational energies of the rigid bender C 3 molecule. The experimental values are shown in parenthesis.

	ν	K	j	p	Energy (cm -1 )
	0	0	0	0	0.00
	0	0	2	0	2.65
	0	0	4	0	8.83
	0	0	6	0	18.53
	0	0	8	0	31.76
	1	1	1	0	66.41 (63.42) 24
	1	1	2	1	68.26
	1	1	3	0	71.02
	1	1	4	1	74.71
	1	1	5	0	79.31
	1	1	6	1	84.84
	1	1	7	0	91.28
	1	1	8	1	98.64
	1	1	9	0	106.92
	2	2	2	0	141.43 (133.07) 14
	2	2	3	1	144.30
	2	2	4	0	148.11
	2	2	5	1	152.89
	2	2	6	0	158.58
	2	2	7	1	165.31
	2	2	8	0	172.82
	2	2	9	1	181.52

Table III .

 III Bound levels of the C 3 -He complex.

	State	RMA	RB-CC	State	RMA	RC-CC
	J	ε	Energy (cm -1 ) Energy (cm -1 ) J	ε	Energy (cm -1 ) Energy (cm -1 )
	0	+	-9,56	-9,73	3	-	-6,95	-7,09
	0	+	-2,89	-3,00	3	-	-5,68	-5,81
	1	+	-3,39	-3,63	3	-	-2,50	-2,68
	1	-	-9,09	-9,26	4	+	-5,44	-5,55
	1	-	-3,85	-4,08	4	+	-3,79	-3,92
	1	-	-2,11	-2,22	4	+	-	-1,26
	2	+	-8,19	-8,35	4	+	-	-1,02
	2	+	-7,07	-7,18	4	-	-4,13	-4,23
	2	+	-3,40	-3,61	4	-	-1,07	-1,01
	2	+	-0,89	-0,98	5	+	-2,05	-2,13
	2	-	-7,10	-7,20	5	-	-3,65	-3,74
	2	-	-2,32	-2,54	5	-	-1,48	-1,62
	3	+	-5,82	-5,92	6	+	-1,58	-1,64
	3	+	-0,77	-0,95

A. A. Penzias, P. M. Solomon, R. W. Wilson, and K. B. Jefferts, Astrophys. J. 168, L53 (1971).

ACKNOWLEDGMENTS

We acknowledge Jose Cernicharo for stimulating this studies and for fruitful discussions. This work has been supported by the Agence Nationale de la Recherche (ANR-HYDRIDES), contract ANR-12-BS05-0011-01 and by the CNRS national program "Physique et Chimie du Milieu Interstellaire." M.H.V. acknowledge the french embassy of Cuba for financial support. F.L. thank the CPER Haute-Normandie/CNRT/Energie, Electronique, Matériaux. Ab initio calculations were performed using High Performance Computing resources of SKIF-Cyberia (Tomsk State University).

ACKNOWLEDGMENTS

Computer time for this study was provided by the Mésocentre de Calcul Intensif Aquitain (MCIA), computing facilities of the Université de Bordeaux, and Université de Pau et des Pays de l'Adour.

CHAPTER 4

Bound states of the HCN-H 2 complex

Introduction to the study

From the astrochemical point of view, HCN and HNC are very interesting molecules. They are among the most abundant organic molecules in the ISM. These molecules have been detected in several molecular clouds [69-73], in planetary atmospheres [74], in comets [75], carbon star atmospheres [76] and circumstellar masers [77].

In this chapter, we focus on HCN. Several investigations of the collision of HCN with He and H 2 have been done both theoretically [78-81] and experimentally [82-85]. However, in all the theoretical studies of the collision with H 2 , the PES have been simplified.

Introduction to the study 6.1.1 Background

The C 3 molecule has been observed in the different regions of the ISM [21, 22, 92-96] and a summary of these detections can be found in several reviews [97, 98]. The study of this molecule in astrochemical conditions is important as it is expected to be one of the building blocks of more complex organic molecules such as carbon chains. It is furthermore a very floppy molecule exhibiting one of the smallest bending frequency (around Rovibrational energy transfer in the He-C 3 collision: potential energy surface and bound states

Otoniel Denis-Alpizar, 1, 2, a) Thierry Stoecklin, 1, b) and Philippe Halvick We present a four-dimensional potential energy surface (PES) for the collision of C 3 with He. Ab initio calculations were carried out at the coupled-cluster level with single and double excitations and a perturbative treatment of triple excitations, using a quadruple-zeta basis set and mid-bond functions. The minimum of the interaction potential is found to be -36.30 cm -1 and corresponds to an almost T-shaped structure of the Van der Waals complex along with a bent configuration of C 3 . This PES is used to determine the rovibrational energy levels of the He-C 3 complex using the rigid monomer approximation (RMA) and the recently developed atom-rigid bender approach at the Close Coupling level (RB-CC). The calculated dissociation energies are -9.56 cm -1 and -9.73 cm -1 respectively at the RMA and RB-CC levels. This is the first theoretical prediction of the bound levels of the He-C 3 complex.

a) Electronic mail: otonieldenisalpizar@gmail.com b) Electronic mail: thierry.stoecklin@u-bordeaux1.fr

We found for this system that the RBAA approach is almost equivalent to exact calculations for pure rotational transitions taking place inside the fundamental bending level ν = 0. On the contrary, for transitions involving two different bending levels, the RBAA approach fails to give an accurate estimate of the magnitude of the cross sections but gives most of the time at least the right ranking of the transition cross sections. We then conclude that l-type transitions cross sections have to be calculated at the RB-CC level for the He-HCN collision while pure rotational transitions cross sections may be calculated accurately at the RBAA level.

The last study presented in this manuscript was dedicated to the C 3 -He collision. The first four dimensional analytical representation of this PES which includes the bending motion of the C 3 molecule was developed by fitting a grid of CCSD(T) ab initio points calculated using a quadruple zeta basis set with mid-bond functions and BSSE correction. The resulting interaction potential has a minimum of -36,30 cm -1 at γ = 120 • , R = 6.57 a 0 , α = 90 • and ϕ = 180 • . The bound states calculations for the C 3 -He complex were performed both at the RMA and at the RB-CC level. The dissociation energy computed using the RMA is 9.56 cm -1 while it is 9.73 cm -1 at the RB-CC level. The agreement between these two values shows that the RMA approach is valid for this system, mainly because the interaction potential is too shallow to couple excited bending states of C 3 . More bound levels of lower energy were found however at the RB-CC level as a consequence of the non-linear geometry of C 3 associated with the minimum of the interaction potential. Perspectives 1. A first important use of the results of this work will be found in the input of the new rate coefficients calculated for H 2 -CS in the astrophysical models of the chemistry of dense and diffuse molecular clouds.

2. The cross section and the rate coefficients for the HCN-H 2 collision using the new PES determined in this work will have to be performed before being also implemented in the same models.

3. An improvement of the RB-CC method will also be considered in the calculations of the l-doubling transitions cross sections detected in the