# Étude théorique de la formation catalytique de petites molécules sur des modèles de grains interstellaires

Ichraf Oueslati

<figure><img src="image_1.png" /></figure>

To cite this version:

Ichraf Oueslati. Étude théorique de la formation catalytique de petites molécules sur des modèles de grains interstellaires. Physique [physics]. Université Pierre et Marie Curie; Faculté des Sciences de Tunis, 2015. Français. NNT: 2015PA066074 . tel-04312209

HAL Id: tel-04312209

https://theses.hal.science/tel-04312209

Submitted on 28 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Thèse présentée pour obtenir le grade de
**Docteur de l'Université Pierre et Marie Curie**

Spécialité: Physique

École doctorale Physique en Ile-de-France

**Docteur de la Faculté des Sciences de Tunis**

Spécialité : Chimie

Par

**Ichraf OUESLATI**

Étude théorique de la formation catalytique de
petites molécules sur des modèles de grains
interstellaires

Thèse soutenue le 13 Avril 2015 devant le jury composé de :

<table><tr><td>M. Olivier Parisel</td><td>Directeur de recherches au CNRS, Université Pierre et Marie Curie</td><td>Président</td></tr><tr><td>Mme Salima Boughdiri</td><td>Maître de conférences, FST, Université Tunis El Manar</td><td>Rapporteur</td></tr><tr><td>Mme Céline Léonard</td><td>Professeur, Université Paris-Est Marne-la-Vallée</td><td>Rapporteur</td></tr><tr><td>M. Mohamed Taieb Ben Dhia</td><td>Maître de conférences, FST, Université Tunis El Manar</td><td>Examinateur</td></tr><tr><td>Mme Boutheïna Kerkeni</td><td>Professeur, ISAMM, Université de la Manouba</td><td>Directrice de thèse</td></tr><tr><td>Mme Lydia Tchang-Brillet</td><td>Professeur, Université Pierre et Marie Curie</td><td>Directrice de thèse</td></tr><tr><td>M. Stefan T. Bromley</td><td>Professeur, Université de Barcelone</td><td>Invité</td></tr></table>
Observatoire de Paris - Meudon
LERMA

Thèse en vue de l'obtention du diplôme de
docteur de l'Université Pierre & Marie Curie
& de la faculté des Sciences de Tunis

Directeurs de thèse : Boutheina KERKENI et Lydia TCHANG-BRILLET
Co-encadreur: Annie SPIELFIEDEL

Etude théorique de la formation
catalytique de petites molécules sur des
modèles de grains interstellaires

Ichraf OUESLATI

<ichraf.oueslati@obspm.fr>
<figure><img src="image_2.png" /></figure>
Remerciements

Mes premiers remerciements vont à mes directrices de thèse, Lydia Tchang-Brillet et Boutheïna Kerkeni, qui ont dirigé cette thèse avec beaucoup de rigueur et de compétences et qui m'ont accordé soutien et disponibilité de tous les instants. Qu'elles trouvent ici toute ma reconnaissance pour leur confiance, leur attention et leurs nombreux conseils qu'elles m'ont fourni tout au long de ces trois années.

Je tiens à remercier particulièrement Nicole Feautrier et Annie Spielfiedel de m'avoir fait partager leur enthousiasme pour la recherche et ainsi de m'avoir donné l'envie de poursuivre dans cette direction. J'ai réellement apprécié de travailler et de discuter avec elles.

Je remercie également Jean Hare, directeur de l'Ecole doctorale de physique de l'Université Pierre et Marie Curie pour m'avoir accueillie avec bienveillance et Laurent Pagani qui m'a parrainée auprès de l'Ecole doctorale et qui a suivi le déroulement de ma thèse avec intérêt.

J'adresse toute ma gratitude à Céline Léonard et à Salima Boughdiri pour avoir accepté de juger ce travail de thèse en tant que rapporteurs. Je les remercie pour la lecture attentive et critique du manuscrit et pour les commentaires pertinents qu'elles m'ont apportés. Je remercie Mohamed Taieb Ben Dhia pour avoir examiné et évalué ce travail et Olivier Parisel de m'avoir fait l'honneur de présider le jury de thèse.

Un grand merci à Jean Louis Lemaire qui m'a offert l'occasion de participer au programme LASSIE FP7 qui m'a permis d'effectuer des stages académiques et industriels, de participer à de nombreux congrès et d'y rencontrer des groupes de chercheurs de toutes les disciplines et de tous les coins du monde.

Mes remerciements vont également à Stefan T. Bromley pour sa grande contribution à la réalisation d'une partie de ce travail, pour l'accueil chaleureux et pour toutes les discussions enrichissantes qu'on a eues ensemble. Je souhaite également remercier toute l'équipe de l'université de Barcelona et en particulier Leny Laksmipriya.

Je remercie également toutes les personnes du laboratoire LERMA et en particulier Fabrice Dayou, Christian Balanca, Benoît Albert, Dominique Lopes, Uddhab Chaulagain, Lionel de Sà et Lisseth Gavilain pour leur accueil, pour leur soutien, pour leur gentillesse et pour l'ambiance chaleureuse qui y règne.

Je remercie aussi la directrice du laboratoire LPMC, Najoua Kammoun, de m'avoir
accueilli dans son laboratoire ainsi que toute l'équipe et en particulier mes amis : Amira
Siai, Ghofrane Ouerfelli, Malek Boukallaba et Mariem Hechmi pour le soutient scientifique
et psychologique qu'ils m'ont apporté.

Enfin, je remercie mes parents et tous les membres de ma famille, qui m'ont soutenu
au cours de cette thèse et tout au long de mes études. Leur présence si chaleureuse et
l'intérêt témoigné à mon égard m'ont été d'un grand secours. Je remercie aussi mon fiancé
Aymen Ezzine qui a été toujours présent pour m'encourager et me soutenir. Ses sacrifices,
son soutien moral, sa gentillesse sans égal et son profond attachement m'ont permis de
réussir ces trois années de thèse. Je termine par remercier tous mes amis.
# Résumé

L'hydrogène moléculaire est la molécule la plus abondante dans l'Univers. Il est reconnu que cette molécule ne peut se former que par un processus catalytique mettant en jeu la surface des grains de poussière interstellaires. Les études théoriques ont été menées principalement sur des surfaces de graphite et très récemment sur la forstérite [010] cristalline et amorphe. C'est dans ce cadre que se situe notre étude qui porte sur les réactions de formation de l'hydrogène moléculaire sur de nouveaux modèles de surfaces silicées et silicatées.

Dans un premier temps, nous avons étudié la réaction d'abstraction par l'hydrogène atomique à partir du tétramethylsilane (un prototype de grains silicés) en phase gazeuse. On a utilisé des méthodes précises de la chimie quantique basées sur la théorie des perturbations à l'ordre 2 ainsi que sur la méthode des clusters couplés. Les effets isotopiques primaire (KIE) et secondaire (SKIE) ont été mis en évidence. La théorie canonique de l'état de transition généralisée (CVT) incluant la correction de l'effet tunnel à faible courbure (SCT) a été appliquée pour étudier la cinétique des réactions dans un intervalle de températures allant de 180 K à 2000 K. Cette étude a montré que la réaction est colinéaire.

Etant donnée l'importance de l'effet tunnel à basse température, nous avons entrepris des calculs de dynamique quantique. Pour ce faire, une surface d'énergie potentielle en coordonnées hyper-sphériques $V (ρ, δ)$ a été construite. Le calcul quantique a été réalisé en utilisant une approche à dimensionnalité réduite à deux dimensions appliquée à ce problème de collision réactive colinéaire. Ce modèle traite explicitement les liaisons qui se rompent et qui se forment lors du processus réactionnel. Les modes spectateurs sont introduits sous forme de l'énergie de point zéro. La résolution des équations couplées permet d'extraire la matrice S de diffusion conduisant aux probabilités d'état à état ainsi qu'aux taux de réactions quantiques. Les résultats montrent que $H_2$ est principalement formé dans son état vibrationnel fondamental. La comparaison avec les résultats expérimentaux réalisés dans le domaine de températures 425-570 K montre un accord satisfaisant. La comparaison des résultats quantiques avec les résultats semi-classiques précédemment obtenus confirme l'importance de l'effet tunnel aux basses températures.

Par ailleurs, une étude de formation de $H_2$ sur des nanosilicates, prototypes des surfaces silicatées, a été menée. Les méthodes de la théorie de la fonctionnelle de la densité (DFT) ont permis d'identifier les sites de physisorption et de chimisorption d'atomes d'hy-
drogène, d'étudier la diffusion de l'atome H physisorbé sur le cluster, de déterminer les caractéristiques énergétiques de ces sites et les énergies d'activation pour la désorption et la recombinaison de H₂. Différents modes d'adsorption ont été identifiés et plusieurs chemins exothermiques et endothermiques pour la recombinaison des deux atomes H ont été obtenus pour ces nanoclusters. Les résultats trouvés ont montré que les grains amorphes/poreux de la composition de la forstérite ont tendance à dissocier H₂ et que les grains de silicate pauvre en Mg (e.g. de la composition de l'enstatite) ou les grains silicatés cristallins/compacts ont tendance à catalyser la formation de H₂.

Ces travaux ont été publiés dans trois articles [1, 2, 3].
Abstract

Molecular hydrogen is the most abundant molecule in the Universe. It was recognized long ago that the formation of molecular hydrogen most likely occurs on dust grains. Theoretical studies have focused on H₂ formation mainly on model graphite surfaces and very recently on polycrystalline and amorphous foresterite [010]. It is within this framework that lies our study addressing the molecular hydrogen formation on new models of siliceous and silicateous surfaces.

First, we studied the abstraction reaction by atomic hydrogen from tetramethylsilane (a prototype of siliceous grains) in the gas phase. We used accurate methods of quantum chemistry based on the second-order perturbation theory and on the coupled clusters method. The primary and secondary kinetic isotope effects (KIEs and SKIEs, respectively) were highlighted. The canonical variational transition state theory (CVT) including the small curvature tunneling (SCT) correction was applied to compute the reaction kinetics for a wide temperature range (180-2000 K). This study showed that the reaction progresses along a collinear pathway with a moderate tunneling effect.

Given the importance of the quantum tunneling effects at low temperatures, we investigated state-selected dynamics using quantum dynamics calculations. In order to achieve this purpose, a two-dimensional potential energy surface in the hyperspherical coordinates representation (ρ, δ) was built. Quantum calculations were performed using a reduced dimensionality approach applied to this collinear reactive collision problem. This model treats explicitly the bonds being broken and formed. Spectator modes are treated adiabatically by inclusion of their projected zero-point energy corrections in the effective potential. Solving coupled equations allows us to extract the scattering S-matrix which leads to state-to-state probabilities as well as the thermal rate constants. The results show that H₂ is mainly formed in its ground vibrational state. The comparison with the experimental results, performed within a temperature range between 425 and 570 K, shows a reasonable agreement. The comparison of quantum and semi-classical results confirms the importance of quantum effects at low temperatures.

Furthermore, the molecular hydrogen formation on nanosilicates, prototypes of silicate surfaces, was investigated. Using Density functional theory (DFT), we identified physisorption and chemisorption sites for hydrogen atoms, then studied the diffusion of physisorbed hydrogen on the nanoclusters and calculated the energy properties and the
activation energies for H₂ recombination and desorption into the gas phase. Different types of adsorption sites were identified. Several endothermic and exothermic reaction pathways with respect to the recombination of two chemisorbed H atoms were obtained for these nanoclusters. The highlighted results showed in general, that amorphous/porous grains with forsteritic composition tend to dissociate H₂ and that the more crystalline/compact silicate grains would tend to catalyse H₂ formation.

The results obtained through the thesis gave rise to three papers published in inter-
national journals [1, 2, 3].
# Table des matières

Introduction 21

I Physico-Chimie moléculaire 25

1 Formation de l'hydrogène moléculaire sur les grains interstellaires 27

1.1 Introduction 27

1.2 Milieu interstellaire, nuages sombres, nuages diffus 27

1.2.1 Milieu Neutre 27

1.2.2 Milieu ionisé 28

1.3 La chimie du milieu interstellaire 29

1.3.1 Réactions en phase gaz 29

1.3.2 Réactions sur les surfaces des grains 30

1.4 Les grains interstellaires 33

1.5 Formation de H₂ dans le milieu interstellaire 35

1.5.1 Formation de H₂ en phase gazeuse 36

1.5.2 Formation de H₂ sur les surfaces de grains interstellaires 36

2 Méthodes de chimie quantique 43

2.1 Introduction 43

2.2 Equation de Schrödinger 43

2.3 Approximation de Born-Oppenheimer 44

2.3.1 Approximation adiabatique 46

2.4 Fonction d'onde électronique 46

2.4.1 Spin-orbitales 47

2.4.2 Déterminants de Slater 48

2.5 Propriétés de l'énergie électronique 49

2.5.1 Principe variationnel 49

2.6 Méthode Hartree Fock 49

2.6.1 Energie Hartree-Fock 49

2.6.2 Méthode du "champ autocohérent", SCF 50
<table>
  <tbody>
    <tr>
      <td>2.6.3</td>
      <td>Les limites du modèle Hartree Fock</td>
      <td>53</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>Méthodes d'interaction de configurations (CI)</td>
      <td>54</td>
    </tr>
    <tr>
      <td>2.7.1</td>
      <td>Méthode Multi Configuration Self-Consistent-Field (MCSCF)</td>
      <td>57</td>
    </tr>
    <tr>
      <td>2.7.2</td>
      <td>Méthode Complete Active Space SCF (CASSCF)</td>
      <td>58</td>
    </tr>
    <tr>
      <td>2.7.3</td>
      <td>Méthode Multi Reference CI (MRCI)</td>
      <td>58</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>Méthode des Clusters Couplés (CC)</td>
      <td>59</td>
    </tr>
    <tr>
      <td>2.8.1</td>
      <td>Méthode des Clusters Couplés CCSD(T)</td>
      <td>59</td>
    </tr>
    <tr>
      <td>2.8.2</td>
      <td>Méthode des Clusters Couplés explicitement corrélés : CCSD(T)-F12</td>
      <td>61</td>
    </tr>
    <tr>
      <td>2.9</td>
      <td>Théorie de la perturbation de Møller-Plesset</td>
      <td>62</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>La théorie de la fonctionnelle de la densité</td>
      <td>64</td>
    </tr>
    <tr>
      <td>2.10.1</td>
      <td>Densité électronique</td>
      <td>65</td>
    </tr>
    <tr>
      <td>2.10.2</td>
      <td>Théorème de Hohenberg et Kohn</td>
      <td>65</td>
    </tr>
    <tr>
      <td>2.10.3</td>
      <td>Théorème de Kohn et Sham</td>
      <td>67</td>
    </tr>
    <tr>
      <td>2.11</td>
      <td>Bases atomiques</td>
      <td>69</td>
    </tr>
    <tr>
      <td>2.11.1</td>
      <td>Bases split valence</td>
      <td>70</td>
    </tr>
    <tr>
      <td>2.11.2</td>
      <td>Fonctions de polarisation</td>
      <td>70</td>
    </tr>
    <tr>
      <td>2.11.3</td>
      <td>Fonctions diffuses</td>
      <td>70</td>
    </tr>
    <tr>
      <td>2.11.4</td>
      <td>Bases explicitement corrélées</td>
      <td>71</td>
    </tr>
    <tr>
      <td>2.11.5</td>
      <td>Bases complètes : CBS</td>
      <td>71</td>
    </tr>
    <tr>
      <td>2.11.6</td>
      <td>Érreur de superposition de base (Basis set superposition error BSSE)</td>
      <td>71</td>
    </tr>
    <tr>
      <td>2.12</td>
      <td>Conclusion</td>
      <td>72</td>
    </tr>
    <tr>
      <td><b>3</b></td>
      <td><b>Dynamique moléculaire</b></td>
      <td><b>73</b></td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Introduction</td>
      <td>73</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Surface d'énergie potentielle (PES)</td>
      <td>75</td>
    </tr>
    <tr>
      <td>3.2.1</td>
      <td>Modes normaux de vibration</td>
      <td>75</td>
    </tr>
    <tr>
      <td>3.2.2</td>
      <td>Chemin d'énergie minimum (MEP)</td>
      <td>77</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Probabilités semi-classiques</td>
      <td>78</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Méthodes de dynamique quantique : Modèle à dimensionnalité réduite (2D)</td>
      <td>80</td>
    </tr>
    <tr>
      <td>3.4.1</td>
      <td>Cas d'une réaction triatomique : A + BC→AB + C</td>
      <td>80</td>
    </tr>
    <tr>
      <td>3.4.2</td>
      <td>Cas d'un système polyatomique quelconque</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.4.3</td>
      <td>Construction de la PES à dimensionnalité réduite (RD)</td>
      <td>86</td>
    </tr>
    <tr>
      <td>3.4.4</td>
      <td>Equations couplées</td>
      <td>90</td>
    </tr>
    <tr>
      <td>3.4.5</td>
      <td>Méthode de propagation de la R-matrice</td>
      <td>92</td>
    </tr>
    <tr>
      <td>3.4.6</td>
      <td>La matrice de diffusion (S)</td>
      <td>94</td>
    </tr>
    <tr>
      <td>3.4.7</td>
      <td>Sections efficaces</td>
      <td>95</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Conclusion</td>
      <td>96</td>
    </tr>
    <tr>
      <td><b>4</b></td>
      <td><b>Cinétique réactionnelle</b></td>
      <td><b>97</b></td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Introduction</td>
      <td>97</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Propriétés thermodynamiques</td>
      <td>97</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>Fonction de partition électronique</td>
      <td>99</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Fonction de partition translationnelle</td>
      <td>100</td>
    </tr>
  </tbody>
</table>
<table><tr><td>4.2.3</td><td>Fonction de partition rotationnelle</td><td>100</td></tr><tr><td>4.2.4</td><td>Fonction de partition vibrationnelle</td><td>101</td></tr><tr><td>4.3</td><td>Approche semi-classique</td><td>101</td></tr><tr><td>4.3.1</td><td>Théorie de l'état de transition</td><td>101</td></tr><tr><td>4.3.2</td><td>Limites de la théorie de l'état de transition</td><td>107</td></tr><tr><td>4.3.3</td><td>Théorie variationnelle canonique de l'état de transition (CVT)</td><td>108</td></tr><tr><td>4.3.4</td><td>Prise en compte de l'effet tunnel</td><td>109</td></tr><tr><td>4.4</td><td>Approche quantique</td><td>112</td></tr><tr><td>4.5</td><td>Conclusion</td><td>115</td></tr></table>

## II Résultats 117

### 5 Réaction d'abstraction d'un hydrogène à partir du tétraméthylsilane : Approche semi-classique 119

<table><tr><td>5.1</td><td>Introduction</td><td>119</td></tr><tr><td>5.2</td><td>Méthodologie</td><td>121</td></tr><tr><td>5.2.1</td><td>Calculs de structure électronique</td><td>121</td></tr><tr><td>5.2.2</td><td>Probabilités de réaction semi-classiques</td><td>122</td></tr><tr><td>5.2.3</td><td>Calcul des constantes de vitesse</td><td>123</td></tr><tr><td>5.3</td><td>Résultats de structure électronique</td><td>126</td></tr><tr><td>5.3.1</td><td>Points stationnaires</td><td>126</td></tr><tr><td>5.3.2</td><td>Analyse vibrationnelle</td><td>128</td></tr><tr><td>5.3.3</td><td>Propriétés énergétiques de réaction</td><td>128</td></tr><tr><td>5.3.4</td><td>Détermination du chemin de réaction</td><td>129</td></tr><tr><td>5.4</td><td>Constantes de vitesse</td><td>129</td></tr><tr><td>5.4.1</td><td>Réaction Si(CH<sub>3</sub>)<sub>4</sub> + H</td><td>131</td></tr><tr><td>5.4.2</td><td>Effets cinétiques isotopiques (KIEs)</td><td>132</td></tr><tr><td>5.5</td><td>Conclusion</td><td>138</td></tr></table>

### 6 Réaction d'abstraction d'un atome H à partir du tétraméthylsilane : Approche quantique 141

<table><tr><td>6.1</td><td>Introduction</td><td>141</td></tr><tr><td>6.2</td><td>Surface d'énergie potentielle ab initio</td><td>142</td></tr><tr><td>6.3</td><td>Dynamique quantique</td><td>146</td></tr><tr><td>6.4</td><td>Résultats et discussion</td><td>149</td></tr><tr><td>6.4.1</td><td>Probabilités de réaction</td><td>149</td></tr><tr><td>6.4.2</td><td>Constantes de vitesse</td><td>151</td></tr><tr><td>6.5</td><td>Conclusion</td><td>154</td></tr></table>

### 7 Adsorption et réactivité de l'hydrogène sur des nanoclusters de silicate 157

<table><tr><td>7.1</td><td>Introduction</td><td>157</td></tr><tr><td>7.2</td><td>Méthodologie</td><td>159</td></tr><tr><td>7.3</td><td>Résultats et discussion</td><td>161</td></tr><tr><td></td><td></td><td>11</td></tr></table>
<table><tbody><tr><td>7.3.1</td><td>Physisorption d'atomes d'hydrogène</td><td>161</td></tr><tr><td>7.3.2</td><td>Chimisorption d'atomes d'hydrogène</td><td>163</td></tr><tr><td>7.3.3</td><td>Adsorption d'un deuxième atome d'hydrogène</td><td>168</td></tr><tr><td>7.3.4</td><td>Formation de H<sub>2</sub></td><td>172</td></tr><tr><td>7.3.5</td><td>Discussion des résultats</td><td>175</td></tr><tr><td>7.3.6</td><td>Conséquence sur la chimie des grains</td><td>178</td></tr><tr><td>7.4</td><td>Conclusion</td><td>178</td></tr><tr><td colspan="2">Conclusion et perspectives</td><td>183</td></tr><tr><td colspan="2">Bibliographie</td><td>185</td></tr><tr><td colspan="2">Annexes</td><td>201</td></tr><tr><td>A</td><td>La théorie de la fonctionnelle de la densité</td><td>203</td></tr><tr><td>A.1</td><td>Approximation de densité locale (LDA)</td><td>203</td></tr><tr><td>A.2</td><td>Approximation du gradient généralisé (GGA)</td><td>204</td></tr><tr><td>A.3</td><td>Functionnelles meta-GGA</td><td>206</td></tr><tr><td>A.4</td><td>Functionnelles hybrides</td><td>206</td></tr><tr><td>A.4.1</td><td>Functionnelle hybride B3LYP</td><td>206</td></tr><tr><td>A.4.2</td><td>Functionnelles hybrides méta-GGA</td><td>206</td></tr><tr><td>B</td><td>Théorie élémentaire de diffusion</td><td>209</td></tr><tr><td>B.1</td><td>Sections efficaces de diffusion</td><td>210</td></tr><tr><td>B.2</td><td>Amplitude de diffusion</td><td>210</td></tr><tr><td>B.3</td><td>Ondes partielles et matrice S</td><td>211</td></tr><tr><td>C</td><td>Discrete Variable Representation (DVR)</td><td>213</td></tr><tr><td>D</td><td>Traitement des rotateurs internes</td><td>217</td></tr><tr><td>D.1</td><td>Rotation libre (free rotor)</td><td>217</td></tr><tr><td>D.1.1</td><td>Rotation contrariée</td><td>219</td></tr><tr><td>E</td><td>Analyse vibrationnelle pour la réaction d'abstraction d'un hydrogène à partir du TMS</td><td>221</td></tr><tr><td>F</td><td>Réaction d'abstraction d'un hydrogène à partir du tétraméthylsilane par un radical hydroxyle : Approche semi-classique</td><td>225</td></tr><tr><td>F.1</td><td>Optimisation des géométries</td><td>226</td></tr><tr><td>F.2</td><td>Analyse vibrationnelle</td><td>229</td></tr><tr><td>F.3</td><td>Energétique de réaction</td><td>229</td></tr><tr><td>F.3.1</td><td>Détermination du chemin de réaction</td><td>230</td></tr><tr><td>F.4</td><td>Constantes de vitesse</td><td>230</td></tr></tbody></table>
<table><tr><td>G Annexe du chapitre 7</td><td>235</td></tr><tr><td>G.1 La relation de Brønsted-Evans-Polanyi</td><td>235</td></tr><tr><td>G.2 Estimation de la BSSE</td><td>235</td></tr><tr><td>G.3 Analyse vibrationnelle des états de transition de la réaction 2H<sub>chim</sub> → H<sub>2,phys</sub></td><td>237</td></tr></table>
Table des figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Diagramme énergétique du phénomène d'adsorption (Lynch, 2001).</td>
      <td>31</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Les trois mécanismes de la chimie de surface. E<sub>b</sub> correspond à la barrière de diffusion d'un site à un autre, et E<sub>D</sub> est l'énergie nécessaire pour désorber la molécule à partir de la surface (Herbst, 2014)</td>
      <td>33</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Cycle de vie de la poussière interstellaire (Bromley et al., 2014).</td>
      <td>34</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Surface d'énergie potentielle du système A + BC → AB + C.</td>
      <td>76</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Barrière de potentiel à une dimension.</td>
      <td>78</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Système de coordonnées de Jacobi pour la réaction colinéaire A + BC → AB + C. « × » dénote l'emplacement des centres de masse de AB et BC.</td>
      <td>82</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Les coordonnées hypersphériques (ρ, δ) définies à partir des coordonnées de Jacobi des réactants (R<sub>r</sub>, r<sub>r</sub>) et des produits (R<sub>p</sub>, r<sub>p</sub>) pour la réaction A + BC → AB + C. Deux « contours plots » (en bleu) de la surface d'énergie potentielle aux faibles valeurs de ρ (région entourant l'état de transition) et aux grandes valeurs de ρ (puits de potentiel des réactants A + BC et des produits AB + C).</td>
      <td>84</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Système de coordonnées colinéaires pour la réaction H<sub>a</sub> + H<sub>b</sub>YZ → H<sub>a</sub>H<sub>b</sub> + YZ.</td>
      <td>86</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Représentation de V(ρ, δ) pour trois valeurs de ρ dans l'espace défini pour la propagation de la R-matrice.</td>
      <td>91</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Diagramme énergétique d'un exemple de réaction dans la théorie de l'état de transition</td>
      <td>102</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Différences d'énergie du point zéro et différences d'énergie d'activation au cours de la rupture des liaisons R-H et R-D analogues.</td>
      <td>107</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>(a) : Surface d'énergie potentielle montrant la MEP (nommé QQ') et le phénomène de « corner cutting » (nommé PP'). (b) : Représentation schématique de la courbure du chemin réactionnel (Wolfsberg, Hook et Paneth, 2009).</td>
      <td>110</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>5.1</td>
      <td>Géométrie optimisée de l'état de transition TS(1.a) au niveau MP2/cc-pVTZ (atome Si : jaune, atomes C : bleu, atomes H : gris)</td>
      <td>126</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Energie potentielle classique (V<sub>MEP</sub>), énergie potentielle vibrationnellement adiabatique (V<sub>a</sub><sup>G</sup>) et ZPE en fonction de s (amu<sup>1/2</sup> a<sub>0</sub>) au niveau de la méthode CCSD(T)/cc-pVTZ//MP2/cc-pVTZ pour la réaction Si(CH<sub>3</sub>)<sub>4</sub> + H.</td>
      <td>131</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Comparaison des constantes de vitesse calculées : k<sub>1.a</sub>(TST), k<sub>1.a</sub>(TST/ZCT), k<sub>1.a</sub>(CVT) et k<sub>1.a</sub>(CVT/SCT) avec les valeurs expérimentales pour la réaction Si(CH<sub>3</sub>)<sub>4</sub>+H.</td>
      <td>132</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Les KIEs calculés par les méthodes TST (a) et TST/ZCT (b) pour Si(CH<sub>3</sub>)<sub>4</sub> + H (k<sub>1.a</sub>/k<sub>1.b</sub>, bleu), Si(CD<sub>3</sub>)<sub>4</sub> + D (k<sub>1.a</sub>/k<sub>1.c</sub>, rouge), Si(CD<sub>3</sub>)<sub>4</sub> + H (k<sub>1.a</sub>/k<sub>1.d</sub>, vert) et Si(CH<sub>3</sub>)<sub>3</sub>CH<sub>2</sub>D +H (k<sub>1.a</sub>/k<sub>1.e</sub>, noir), en fonction de la température.</td>
      <td>135</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Les KIEs calculés par la méthode CVT/SCT pour Si(CH<sub>3</sub>)<sub>4</sub> + D (k<sub>1.a</sub>/k<sub>1.b</sub>, Si(CD<sub>3</sub>)<sub>4</sub> + D (k<sub>1.a</sub>/k<sub>1.c</sub>, Si(CD<sub>3</sub>)<sub>4</sub> + H (k<sub>1.a</sub>/k<sub>1.d</sub>) et Si(CH<sub>3</sub>)<sub>3</sub>CH<sub>2</sub>D +H (k<sub>1.a</sub>/k<sub>1.e</sub>), en fonction de la température.</td>
      <td>136</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Constantes de vitesse en fonction de 1000/T pour la réaction Si(CH<sub>3</sub>)<sub>4</sub> + H et les réactions substituées calculées par les méthodes TST/ZCT (a) et CVT/SCT (b) : k<sub>1a</sub> bleu, k<sub>1c</sub> rouge, k<sub>1d</sub> noir, et 12 × k<sub>1e</sub> vert.</td>
      <td>137</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Coordonnées utilisées pour le système H+Si(CH<sub>3</sub>)<sub>4</sub>. atomes Si : jaune, atomes C : bleu et atome H : gris.</td>
      <td>143</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>(a) Courbe de l'énergie potentielle V (en Hartree) en fonction de δ (rad) obtenue par ajustement et reliant les points ab initio calculés pour ρ=3.31 a.u. et ρ=10.18 a.u.. (b) Contours de V (ρ, δ). Les énergies des contours et des points ab initio (en astérisque) varient entre 0.122 a.u. et 0.158 a.u. Les positions du TS et des complexes de van der Waals sont indiquées avec des flèches.</td>
      <td>144</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Probabilité cumulative de réaction calculée J = 0 en fonction de l'énergie translationnelle (définie comme la différence entre l'énergie totale E et l'énergie vibrationnelle de l'état fondamental E(0)) : (a) pour les états initiaux s' = 0 et s' = 1 du TMS, (b) sommée sur tous les états initiaux et finaux.</td>
      <td>149</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Probabilités d'état-à-état pour un moment magnétique total (J=0) pour la réaction H+Si(CH<sub>3</sub>)<sub>4</sub> (s') → H<sub>2</sub> (s)+Si(CH<sub>3</sub>)CH<sub>2</sub> et différentes valeurs de s et s' en fonction de l'énergie translationnelle E - E(s' = 0). (a) : s' = 0; (b) : s' = 1.</td>
      <td>150</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>Comme Fig. 6.4 pour (c) :s' = 2; (d) :s = 2.</td>
      <td>150</td>
    </tr>
    <tr>
      <td>6.6</td>
      <td>Adiabats hypersphériques, ϵ(ρ<sub>i</sub>) pour les premiers états des réactants et des produits en fonction de ρ. Les courbes pleines correspondent aux états des réactants et les courbes en pointillé aux états des produits.</td>
      <td>152</td>
    </tr>
  </tbody>
</table>
Table des figures

<table>
  <tbody>
    <tr>
      <td>D.1</td>
      <td>Angle de torsion φ entre deux groupements : (CH<sub>3</sub>)H et Si(CH<sub>3</sub>)<sub>3</sub>.</td>
      <td>218</td>
    </tr>
    <tr>
      <td>F.1</td>
      <td>Géométrie optimisée de l'état de transition TS(2) au niveau MP2/cc-pVTZ (atomes Si : jaune, C : bleu, H : gris et O : rouge)</td>
      <td>228</td>
    </tr>
    <tr>
      <td>F.2</td>
      <td>Energie potentielle classique (V<sub>MEP</sub>), énergie potentielle vibrationnelle-ment adiabatique (V<sub>a</sub><sup>G</sup>) et ZPE en fonction de s (amu<sup>1/2</sup> a<sub>0</sub>) calculées au niveau CCSD(T)/cc-pVTZ//MP2/cc-pVTZ pour la réaction Si(CH<sub>3</sub>)<sub>4</sub> + OH.</td>
      <td>231</td>
    </tr>
    <tr>
      <td>F.3</td>
      <td>Comparaison des constantes de vitesse calculées avec les valeurs expéri-mentales et théoriques pour la réaction Si(CH<sub>3</sub>)<sub>4</sub>+OH.</td>
      <td>232</td>
    </tr>
    <tr>
      <td>G.1</td>
      <td>Géométries optimisées au niveau MPWB1K des états de transition corres-pondant aux labels « TS » dans la figure (7.10). Les indices se rapportent aux figures montrant les configurations initiales pour deux atomes d'hy-drogène chimisorbés sur les clusters de nano-pyroxène et de nano-olivine<sub>a</sub>. Atomes Si : jaune, C : bleu, H : gris et O : rouge.</td>
      <td>239</td>
    </tr>
  </tbody>
</table>
Introduction

Le milieu interstellaire (MIS) entre les étoiles et les galaxies est un lieu privilégié où des
agrégats, des poussières et plus de 160 espèces moléculaires évoluent dans des conditions
physiques (température et densité) hors équilibre thermodynamique [4, 5, 6, 7], mettant
en jeu des processus physico-chimiques variés. La spectroscopie fine dans les domaines
infrarouge et submillimétrique a permis d'accéder à la physique et à la chimie de ces
constituants, qu'ils soient en phase gazeuse, solide, ou à l'interface de ces deux phases. En
dépit d'avancées significatives, la caractérisation des processus clés pour la chimie du MIS
reste cependant largement incomplète. La compréhension de ces mécanismes complexes à
l'origine de cette évolution représente donc un défi majeur. La prédiction des propriétés
physico-chimiques requises pour une modélisation fiable de la chimie de ce milieu est
aussi cruciale.

L'enchaînement des réactions chimiques en phase gazeuse conduit à la formation de
molécules interstellaires, même les plus complexes. Cette production requiert cependant
la présence de l'hydrogène moléculaire. Son abondance joue également un rôle essentiel
sur l'évolution dynamique du MIS, notamment sur les processus de formation d'étoiles
dans les milieux diffus. Dans les conditions extrêmes de ce milieu, H₂ ne peut se former en
phase gazeuse. Le mécanisme le plus efficace pour la formation de cette molécule est celui
de l'association d'atomes d'hydrogène sur les grains de poussière interstellaire [8, 9, 10]
qui agissent comme un catalyseur absorbant une partie de l'énergie libérée (434 kJ/mol).
La compréhension de ce processus reste jusqu'à aujourd'hui à élucider et notamment
les mécanismes mettant en jeu les surfaces silicatées. Ces dernières présentent l'inconvé-
nient d'avoir des surfaces assez complexes, d'avoir des défauts cristallins et d'avoir un
taux de cristallinité différent selon les régions du MIS considérées. Ces caractéristiques
structurales sont à l'origine de certaines propriétés intéressantes qui peuvent affecter la
formation de l'hydrogène moléculaire. Un autre élément manque à la compréhension de
cette réaction c'est la répartition de l'énergie libérée lors de sa formation, sous forme
d'énergie interne (excitation ro-vibrationnelle), énergie cinétique ou échauffement du
grain [11, 12, 13, 14].

De nombreux travaux expérimentaux ont déjà été réalisés sur les surfaces carbonées
[15, 16], sur les silicates amorphes [17, 18, 19] et cristallins [20], ainsi que sur les surfaces recouvertes de glaces [21, 22, 23]. Toutefois, plusieurs questions demeurent sans réponse en ce qui concerne l'analyse des résultats obtenus. La modélisation théorique dans ses différentes formes est considérée comme un appui important à cet effort puisqu'elle permet de donner un aperçu plus détaillé de la chimie de la poussière interstellaire (la structure, la composition, les barrières de réaction). Fournir ces détails qui sont difficiles à déterminer avec précision à partir des observations et de l'expérience représente l'enjeu majeur pour les études théoriques.

Compte tenu de la relative simplicité de sa structure, plusieurs travaux théoriques ont été consacrés à l'étude des interactions H avec des modèles de grains de poussières carbonées [24, 25, 26, 27], en particulier les surfaces de graphite ou graphène ou hydrocarbure aromatique polycyclique (PAH). Selon Garcia-Gil et al. [28], les résultats obtenus avec les surfaces de graphite montrent que ces dernières ne fournissent pas les conditions requise pour une formation efficace de $H_2$. Par conséquent, il est utile de considérer les silicates, qui sont les autres constituants majeurs des grains de poussière interstellaire. Les mécanismes réactionnels mettant en jeu des surfaces silicatées ont été peu étudiés théoriquement jusqu'ici. Des études très récentes sur les surfaces de forstérite cristalline (010) [29, 14, 28, 30] et sur des nanoclusters qui sont à l'origine la formation des silicates plus gros par voie de nucléation [31, 32], mettent en évidence l'importance d'une voie de formation par chimisorption sans barrière. Ce résultat pourrait avoir un grand impact sur l'astrochimie. C'est dans ce contexte que notre travail s'inscrit. Il concerne l'étude des mécanismes de formation de l'hydrogène moléculaire $H_2$ à la surface des clusters de nano-olivine, nanopyroxène et tétraméthylsilane (TMS) qui sont des prototypes de grains secs silicatés et silicés, respectivement.

Dans ce mémoire, nous commençons par détailler les processus physico-chimiques utiles à la compréhension actuelle des phénomènes ainsi que les prérequis astrophysiques montrant l'importance du processus de formation de $H_2$ notamment sur les surfaces des grains interstellaires. Le second chapitre décrit des méthodes sophistiquées de calcul de chimie quantique, consacrées à la résolution de l'équation de Schrödinger électronique pour les systèmes moléculaires polyélectroniques, ainsi que les bases utilisées pour décrire les orbitales atomiques. La première partie du troisième chapitre est dédiée à la description d'une surface d'énergie potentielle (PES), solution de l'équation électronique et impliquée dans une collision réactionnelle, de ses points caractéristiques et notamment le chemin réactionnel intrinsèque et son rôle dans le calcul des probabilités de réaction dans la théorie semi-classique. La deuxième partie de ce chapitre est consacrée à l'étude de la dynamique de l'équation de Schrödinger nucléaire, qui permettra d'introduire la surface d'énergie potentielle en fonction des coordonnées hypersphériques pour un système à dimensionnalité réduite traitant explicitement deux degrés de liberté (2D) en se basant sur des méthodes ab initio. L'étude de la dynamique quantique fondée sur la théorie de la diffusion y est aussi présentée. L'étude de la cinétique réactionnelle fait l'objet du quatrième chapitre où est décrit le calcul des constantes de vitesse de réaction soit par
les méthodes fondées sur la théorie de l'état de transition incluant les approximations semi-classiques, soit par la méthode quantique.

Avant d'aborder l'étude des réactions sur les surfaces silicatées, nous avons considéré la réaction d'abstraction d'un atome d'hydrogène à partir du tétraméthysilane (TMS). Ce travail peut être considéré comme une étape vers l'étude de réactions impliquant des molécules de plus en plus complexes. Le cinquième chapitre décrit donc les méthodes de calcul employées pour l'étude de l'abstraction à partir du TMS ainsi que la discussion des résultats obtenus concernant l'optimisation des géométries des fragments, les calculs des propriétés énergétiques, les fréquences de vibration et le calcul des constantes de vitesse par des approches semi-classiques [1]. Une étude isotopique a aussi été menée. La construction d'une surface d'énergie potentielle, obtenue par des calculs *ab initio* précis suivie par une méthode d'ajustement analytique (fit), a permis l'interprétation des résultats obtenus par l'étude de dynamique quantique (2D), tels que le calcul des probabilités de réaction cumulatives et des probabilités d'état-à-état, ainsi que des constantes de vitesse (taux de réaction)[2].

Le dernier chapitre expose les études DFT réalisées afin d'identifier les structures géométriques et les chemins énergétiques dans le processus d'adsorption, diffusion, recombinaison et désorption de H₂ sur des grains de poussière amorphe de silicate. Récemment, plusieurs études ont été menées sur la surface cristalline de forstérite et sur un cluster de nanosilicate possèvant stoechiométrie de la forstérite afin de caractériser les interactions avec des atomes d'hydrogène et d'étudier la formation de l'hydrogène moléculaire. Ces travaux ont montré que la chimisorption des atomes d'hydrogène se fait sans barrière et que le mécanisme le plus probable suit le mécanisme de Langmuir-Hinshelwood où deux atomes d'hydrogène chimisorbés se recombinent pour former la molécule H₂. Grâce à leurs faibles tailles, les nano-clusters offrent plus de voies à la formation de H₂. Les barrières à la recombinaison obtenues à partir des différents sites identifiés montrent une relation de Bronsted-Evans-Polanyi (BEP) en fonction de l'énergie de liaison de deux atomes d'hydrogènes chimisorbés sur le grain indépendemment de la taille ou de la cristallinité/amorphicité de ce dernier [3]. La surface cristalline et les clusters ont été utilisés pour modéliser les grains interstellaires.
# Première partie

## Physico-Chimie moléculaire
<figure><img src="image_3.png" /></figure>
Chapitre 1

# Formation de l'hydrogène moléculaire sur les grains interstellaires

## 1.1 Introduction

Le milieu interstellaire, MIS, est le milieu qui sépare les étoiles dans la voie lactée. Il possède une structure complexe, composée de poussières interstellaires (1%) [4] baignant dans un gaz constitué principalement d'hydrogène qui représente 90% du gaz total, de 9% d'hélium et d'environ 1% d'éléments lourds comme C, N, O, Si, Fe, Mg et S. Les grains jouent un rôle important dans la physique et la chimie du MIS puisqu'ils peuvent absorber et diffuser la lumière des étoiles, catalyser des réactions qui ne peuvent avoir lieu dans la phase gaz, entrer en collision avec des atomes et des ions du gaz ou d'autres grains et libérer des éléments lourds.

Ce chapitre a pour but d'expliquer les motivations qui ont amené ce sujet de recherche. Il est composé de quatre parties. La première comporte une description du milieu inter-stellaire, sa composition ainsi qu'une description des différentes zones qui le constituent. La seconde partie porte sur la chimie du milieu interstellaire. La troisième partie donne un aperçu sur la formation et l'évolution des grains interstellaires. Enfin la dernière partie a pour objectif de présenter les différents travaux expérimentaux et théoriques ayant porté sur la molécule H₂.

## 1.2 Milieu interstellaire, nuages sombres, nuages diffus

Le MIS comprend de nombreux environnements différents, avec des propriétés physiques différentes. Ils sont classés en deux grandes catégories : neutre et ionisé.

### 1.2.1 Milieu Neutre

Le milieu interstellaire neutre ne contient pas d'ions hydrogène, mais uniquement de l'hydrogène sous forme atomique ou moléculaire. Ce milieu peut être séparé en deux
parties présentant des caractéristiques très différentes.

Nuages diffus

Entourés de régions plus denses et plus froides, les nuages diffus sont composés principalement de l'hydrogène atomique neutre H noté HI. C'est un milieu atomique neutre froid (température varie entre 50 et 200 K). Les principales molécules diatomiques qui ont été découvertes dans ce milieu sont principalement H₂, CO, SiO, CS, CH, CH⁺, CN et OH mais aussi des molécules plus complexes telles que HCO⁺, NH₂⁺ et H₂CO.
[33] Ces molécules sont aussi présentes dans les comètes et dans des galaxies extérieures (notamment la molécule CO). D'autres grosses molécules résistantes aux UV, telles que les hydrocarbures polycycliques aromatiques (PAHs) ou d'autres clusters de carbone, ont été observées. [34] En dehors de la formation de H₂, la chimie sur les grains ne joue pas un rôle important. On a plutôt une chimie en phase gaz.

Nuage moléculaire

Ils se forment lors de la contraction des nuages HI. Ce milieu moléculaire neutre est une des régions les plus froides du milieu interstellaire (la température peut descendre plus bas que 10 K). Une grande proportion des poussières interstellaires se trouve dans ces nuages denses qui absorbent les rayonnements UV. De ce fait, on peut généralement considérer ces nuages comme des zones opaques et sombres dans le milieu interstellaire. Dans les nuages moléculaires, la réaction la plus importante est celle de la formation de H₂ à la surface des grains de poussière par recombinaison de deux atomes d'hydrogène. Cette réaction est fortement exothermique (-434 kJ/mol soit -4.48 eV). Cette énergie qui représente l'énergie de liaison se répartit entre l'excitation vibrationnelle et rotationnelle de H₂, le chauffage du grain et l'énergie cinétique de la molécule quittant le grain.
[35, 11, 12, 14]

1.2.2 Milieu ionisé

Dans le MIS ionisé, le gaz interstellaire est ionisé par les rayonnements UV lointains des étoiles chaudes, par des mécanismes comme l'ionisation collisionnelle lors des chocs, par des rayons X et par des particules chargées de haute énergie. L'espèce majoritaire de ce milieu est H⁺. Ce milieu est constitué de trois phases dont la première est la région H II (nébuleuses gazeuses), la seconde est le milieu interstellaire ionisé diffus et la dernière est le milieu interstellaire chaud issu de la dispersion des restes de supernova et des bulles.

Les régions HII :

Idéalement, cette région se présente comme un grand nuage d'hydrogène entourant une étoile chaude ou un amas d'étoiles chaudes qui émettent des UV, provoquant ainsi l'ioni-sation du gaz contenu dans ce milieu. Ceci explique l'absence d'hydrogène atomique. On trouve ces nébuleuses à proximité des nuages moléculaires. Une des régions HII les plus connues est la nébuleuse d'Orion.
Les régions HII émettent aussi un continu infrarouge. En effet, les poussières qui s'y trouvent sont chauffées par le rayonnement des étoiles excitatrices, et réémettent de l'infrarouge.

**Milieu ionisé diffus :**

Le MIS contient du gaz ionisé en dehors des régions HII. Ce gaz provient des fuites hors de ces régions ou de l'ionisation de l'hydrogène par des rayonnements des étoiles chaudes. Ce dernier mode d'ionisation est relativement très faible. Ce milieu est de faible densité et la température est estimée à 8000 K. L'hydrogène y est entièrement ionisé.

**Gaz chaud :**

Des observations ont montré l'existence d'un gaz diffus à des températures de quelques $10^5$ à $10^7$ K. Il remplit une grande partie du milieu interstellaire et il est peu dense. Il provient des restes de supernova et des bulles. Ce gaz peut se répandre à de grandes distances du disque galactique formant un halo (on peut isoler trois zones distinctes dans une galaxie, le bulbe au centre, le disque galactique et le halo).

## 1.3 La chimie du milieu interstellaire

Les conditions physiques du MIS ont des implications profondes sur les réactions chimiques possibles. Les réactions exothermiques peuvent avoir lieu aux faibles températures des nuages moléculaires et du milieu atomique froid. Certaines réactions devant franchir une barrière de potentiel ne sont pas favorisées dans le cas où les températures sont très basses. Les réactions ion-molécule forment des ions moléculaires de plus en plus complexes qui réagissent à leur tour pour former des molécules neutres. Les molécules peuvent aussi se photodissocier en présence de photons.

### 1.3.1 Réactions en phase gaz

Les conditions extrêmes dans les nuages interstellaires font que le nombre de processus efficaces est limité. Les grandes familles de réactions qui peuvent intervenir en phase gaz sont énumérées.

-Les réactions radiatives : ce sont une combinaison entre deux entités ionisées ou neutres avec formation d'une molécule et évacuation d'énergie en excès par émission d'un photon.

-Les réactions ion-molécule : en se rapprochant de la molécule neutre, les ions créent un dipôle électrique induit. La force attractive engendrée favorise ainsi la réaction. Ces réactions sont possibles dans les milieux à très basses températures à condition qu'elles soient exothermiques et de faible barrière d'activation.

-La recombinaison dissociative : des ions moléculaires se recombinent avec un électron pour former des molécules neutres dans un état dissociatif.

-Les réactions neutre/neutre : elles sont dues à l'interaction atome/atome ou molé-
cule/molécule neutres. Elles ont lieu dans le milieu interstellaire neutre et jouent un rôle important dans les nuages moléculaires.

- La photodissociation : ce mécanisme conduit à la destruction des molécules dans le MIS diffus et dans les parties externes des nuages moléculaires où le rayonnement UV peuvent pénétrer. Néanmoins, des photons ultraviolets peuvent y être crées par l'excitation de H₂ ou de He par les rayons cosmiques ou par collisions avec des électrons secondaires.

## 1.3.2 Réactions sur les surfaces des grains

La chimie de surface joue un rôle crucial dans la formation et la destruction des molécules dans le MIS. Elle est à l'origine de plusieurs réactions qui ne peuvent pas avoir lieu dans la phase gaz. Un exemple qui illustre ce phénomène est la formation de H₂ et de CO₂. Ce dernier a été identifié dans les manteaux de glace recouvrant les grains dans les parties profondes des nuages moléculaires.

De façon générale, les réactions sur les surfaces (substrats) impliquent une étape d'adsorption d'un ou plusieurs réactants. Deux types d'adsorption existent suivant que la particule incidente (adsorbat) subisse une interaction de nature physique (physisorption) ou bien qu'une liaison chimique ait été formée (chimisorption). La réaction met en jeu également des processus de diffusion et de désorption.

### La physisorption

La faible interaction du cortège électronique de la particule incidente et des atomes de la surface suffit à permettre la cohésion de l'ensemble. Les forces mises en jeu dans l'adsorption physique sont du type Van der Waals (interaction dipôle-dipôle). Le terme attractif entre la particule libre et le solide résulte de la somme des interactions de cette particule avec un grand nombre d'atomes ou de molécules constitutifs du solide au voisinage de la surface. Cette attraction conduit à un minimum de potentiel qui représente l'énergie d'adsorption exprimée par rapport à l'état de référence où l'interaction est nulle (distance infinie). Comme l'énergie augmente avec la distance, la particule va spontanément se rapprocher du solide et se stabiliser à la distance correspondant au minimum : c'est l'état physisorbé. A cette physisorption peut succéder une désorption lorsque la particule possède assez d'énergie pour « remonter » le potentiel et quitter la surface.

### La chimisorption

Dans le cas de la chimisorption, une liaison chimique se forme entre l'espèce adsorbée et un ou plusieurs atomes superficiels de la surface : c'est une véritable réaction chimique. L'existence d'une telle réaction impose des conditions de localisation de l'adsorption à l'endroit où la formation de la liaison est la plus favorable : c'est un site d'adsorption. Les diagrammes de Lennard-Jones (Figure 1.1) [36] représentent clairement les différents domaines correspondant à l'adsorption physique et chimique. L'énergie potentielle du système est portée en fonction de la distance adsorbat-surface. L'adsorption chimique est
<figure><img src="image_4.png" /><figcaption>FIGURE 1.1 – Diagramme énergétique du phénomène d'adsorption (Lynch, 2001).</figcaption></figure>

qualifiée de dissociative (la molécule est adsorbée sous forme dissociée) ou associative (les atomes se combinent à la surface pour former la molécule).

La diffusion

Ce mouvement est caractéristique d'un processus thermique. La fréquence de ce phénomène augmente avec la température. Il dépend aussi de la barrière à la diffusion. Cette mobilité résulte principalement de l'agitation thermique et de l'effet tunnel qui rend possible le passage d'un site de physisorption à un site de chimisorption. En effet, les atomes ou molécules adsorbés sur la surface peuvent se déplacer pour atteindre d'autres sites d'adsorption leur permettant d'augmenter leur temps de résidence sur les grains. Ils peuvent aussi rencontrer d'autres atomes ou molécules pour réagir. Trois mécanismes (Figure 1.2) [7] ont été identifiés :

-Le processus Langmuir-Hinshelwood (LH) est un mécanisme de diffusion d'abord suggéré par Langmuir en 1921 et développé par Hinshelwood en 1926 [37]. Dans ce mécanisme, les réactifs sont initialement adsorbés (physisorbés ou chimisorbés) sur la
surface et en équilibre thermique avec le substrat. Le mécanisme de ce processus peut être résumé selon cette équation :

$$A_{ads} + B_{ads} \xrightarrow{diffusion} AB_{ads} \xrightarrow{desorption} AB_{g} \quad (1.1)$$

C'est le mécanisme le plus étudié pour la formation de H₂ dans l'hypothèse où le nombre d'atomes d'hydrogène adsorbés sur un grain est faible. [7] Les atomes d'hydrogène diffusent entre les différents sites actifs en surmontant une barrière de potentiel ou par effet tunnel. Les réactions peuvent se produire lorsque deux espèces atteignent le même site actif. Des expériences de TPD (Thermal Program Desorption) menées à faible température sur le carbone amorphe et sur les silicates [15] ont montré la faiblesse de la liaison d'adsorption des atomes d'hydrogène ce qui laisse penser qu'il y a une barrière à la chemisorption.

En plus du flux d'atomes d'hydrogène, la formation de H₂ via le mécanisme LH dépend aussi de la nature des grains, qu'ils soient carbonés ou silicatés et de la température du milieu [7]. En effet, à des faibles températures, les atomes d'hydrogène ne peuvent plus se déplacer alors qu'à des températures plus élevées, ils désorbent avant qu'ils puissent réagir ensemble. Pour que le mécanisme LH soit opérationnel, le processus doit donc se faire à des températures non faibles [7]. Par exemple, 50% de l'abondance de l'hydrogène moléculaire se trouvent dans les nuages diffus dont la température varie de 20 à 25 K [7]. En plus, pour que le mécanisme de LH soit efficace, les surfaces des grains doivent présenter des sites d'adsorption avec une variété d'énergies de liaison [7].

**Le processus Eley-Rideal** est un mécanisme direct, proposé par Eley et Rideal en 1940 [38]. Seul l'un des réactifs est fortement adsorbé (chimisorbé) sur la surface du grain. La réaction a lieu quand un réactif en phase gazeuse entre en collision directement avec celui déjà adsorbé sur la surface. A des faibles températures du MIS (≈ 5 K), la vitesse de diffusion des atomes d'hydrogène est assez faible pour qu'un grand nombre d'entre eux restent relativement immobiles sur la surface [7]. C'est dans ces conditions que le mécanisme d'Eley-Rideal, dans lequel Hg heurte directement un Hads et réagit avec lui, devient dominant dans la mesure où il n'existe pas de barrière à la réaction. Dans ces régions, cependant, la surface est susceptible d'être recouverte de glace d'eau.

**Le processus Harris-Kasemo (HK)** c'est un mécanisme pour lequel une partie de l'énergie cinétique d'un atome incident lui sert à diffuser sur la surface du grain et ainsi rencontrer un autre atome adsorbé. Ce processus, nommé *Hot-Atom*, se distingue du mécanisme LH par l'origine non thermique de la mobilité du réactant. Les processus LH et HK peuvent se produire en même temps pour former H₂ aux faibles températures [11, 7]. Les réactions suivant le mécanisme HK peuvent se faire même si l'atome d'hydrogène non thermalisé doit parcourir une longue distance pour réagir avec un Hphys [7]. Cependant cette distance dépend de la température de la phase gaz et des propriétés de la surface (énergie d'adsorption, morphologie) [7]. D'après une étude récente de Lemaire et al [11], la formation de D₂ sur les silicates ne peut se faire que par le processus HK pour une température de 12 K.
<figure><img src="image_5.png" /><figcaption>FIGURE 1.2 – Les trois mécanismes de la chimie de surface. $E_b$ correspond à la barrière de diffusion d'un site à un autre, et $E_D$ est l'énergie nécessaire pour désorber la molécule à partir de la surface (Herbst, 2014)</figcaption></figure>

### La désorption

Elle limite le temps de résidence des adsorbats ce qui les empêche de diffuser et de réagir avec d'autres atomes adsorbés. Comme l'étape d'adsorption, ce processus est considéré comme l'étape déterminante de la cinétique chimique : les énergies d'activation des deux phénomènes sont en effet relativement élevées.

## 1.4 Les grains interstellaires

Les grains de poussières sont formés dans les enveloppes d'étoiles en fin de vie et ils sont intimement mélangés au gaz interstellaire. Une représentation schématique du cycle de vie de la poussière interstellaire est donnée dans la figure (1.3) [39].

Ils représentent 1% de la matière interstellaire et possèdent une taille typique de 0.1 µm [4]. Ils sont caractérisés par les propriétés suivantes : [40]

- Ils absorbent une partie des photons émis par les étoiles et ils réémettent l'énergie absorbée par rayonnement thermique dans l'infrarouge moyen et lointain.

- Les atomes ou les molécules qui interagissent avec les poussières interstellaires peuvent former un manteau de glace ou réagir pour former d'autres molécules. Les grains jouent ainsi le rôle de catalyseur de réactions qui sont impossibles dans la phase gaz.

- Ils libèrent des électrons après absorption de photons UV. Ces électrons constituent une source de chauffage du MIS.

- Après destruction des grains (fragmentés par évaporation et par pulvérisation cathodique [41] suite aux collisions des grains dans les ondes de choc de supernova et dans la phase gaz respectivement), des éléments lourds sont libérés dans le gaz interstellaire. Ces atomes peuvent se condenser sur des grains formant ainsi des manteaux réfractaires.

La nature de l'étoile en fin de vie à partir de laquelle sont formés les grains influe sur la composition de ceux-ci. En effet, les étoiles en fin d'évolution sont riches soit en oxygène soit en carbone. Le spectre des étoiles carbonées présente principalement
<figure><img src="image_6.png" /><figcaption>FIGURE 1.3 – Cycle de vie de la poussière interstellaire (Bromley et al., 2014).</figcaption></figure>

une bande due aux grains de carbure de silicium (SiC) et une autre due au carbone
amorphe hydrogéné. Le spectre des enveloppes d'étoiles riches en oxygène a montré des
bandes caractéristiques des silicates amorphes et des silicates cristallins. Les données ISO
(Infrared Space Observatory) ont révélé que (i) la plupart de ces silicates formés sont
riches en magnésium [42] (ii) une fraction de 10% est crystalline [42] et (iii) l'enstatite
(MgSiO₃) est plus abondant que la forstérite (Mg₂SiO₄) [42, 43].

Dans le MIS, les amas formés à partir des étoiles en fin de vie correspondent aux noyaux de condensation des grains de poussière. Comparés à ceux des enveloppes des étoiles, leurs spectres présentent des modifications dont les plus importantes sont les bandes des carbones amorphes hydrogénés partiellement aromatiques et des silicates amorphes (principalement des olivines) [44, 45].

Les silicates sont vraisemblablement amorphisés durant leur séjour dans le milieu
interstellaire par des ondes de choc qui brisent les grains par pulsation ionique[44, 45]. Ce
processus change les propriétés physiques d'un solide et notamment la structure cristalline
qui est détruite. Ces mêmes poussières sont transformées pendant leur incorporation dans
les étoiles jeunes (Young Stellar Objects) [31] puisque l'on retrouve des signatures cristal-
lines dans les disques protoplanétaires [46]. Ces différentes observations montrent que les
silicates évoluent à chaque étape de leur cycle de vie. Pourtant, les poussières silicatées
sont le plus souvent présentes dans des milieux froids dans lesquels peu d'évolutions sont à
attendre par des processus thermiquement activés. De petites tailles, elles sont néanmoins
fortement exposées aux rayonnements et aux flux de particules qui balayent constamment
le milieu interstellaire. Ainsi, les phénomènes d'irradiation pourraient être à l'origine des
évolutions de la poussière silicatée [47]. Néanmoins dans les régions de HII, une quantité
notable de silicates cristallins est identifiée [48]. Ils sont formés suite au chauffage des
grains par le rayonnement intense des étoiles excitatrices. Le taux de cristallinité a été
estimé à 2,2 % [44] contre 10 à 40 % autour des étoiles en fin de vie. La poussière silicatée
amorphe est constituée d'environ 85 % de silicates du type olivine probablement enrichis
en fer et de 15 % de silicates du type pyroxène probablement enrichis en magnésium
[44]. Dans les nuages froids, dans lesquels la poussière n'est pas exposée aux chocs, les
silicates ont une composition du type olivine riche en fer et riche en magnésium [49]. La
composition des silicates est contrainte par les abondances cosmiques. Les plus simples
sont :

- L'olivine est constitué de minéraux de formule générale $XYSiO_4$ comprise entre deux pôles : la forestérite ($Mg_2SiO_4$) et la fayalite ($Fe_2SiO_4$). Le X et Y sont les cations $Mg^{+2}$, $Fe^{+2}$ ou $Ca^{+2}$. Ils sont répérés également par des bandes d'extinction.

- Le pyroxène est un silicate naturel de calcium, de magnésium et de fer cristallisant dans les systèmes orthorhombique et monoclinique. La structure des pyroxènes est caractérisée par un assemblage de chaînes simples de complexes tétraédriques $SiO_4$. Ils forment deux grandes familles : les orthopyroxènes ($FeMg)_2Si_2O_6$, dont la plus abondante dans le MIS est l'enstatite $MgSiO_3$, et les clinopyroxènes ($CaFeMg)_2Si_2O_6$.

Les grains carbonés subissent des modifications importantes dans le MIS. Des indications montrent l'existence d'émetteurs de bandes IR aromatiques formées à partir des grains de carbone amorphe hydrogéné tels que les hydrocarbures aromatiques polycycliques (PAH). Ce sont des molécules contenant des cycles aromatiques carbonés, typiquement des cycles benzéniques. Ils possèdent des bandes aromatiques observées dans le MIS [50] et peuvent constituer jusqu'à 20% de la matière carbonée interstellaire.

Les différentes sortes de grains peuvent être recouverts de manteau de glace lorsqu'ils sont dans les nuages denses. De nombreuses composantes de ces glaces ont été identifiées telles que $H_2O$, $CO$, $CH_2OH$, $CH_3OH$, $NH_3$. Les grains recouverts de glace ont des températures très faibles de l'ordre de 10-30K. Leur morphologie n'est pas encore correctement connue. Cependant en comparant les spectres IR de films de glace d'eau obtenus en laboratoire et ceux du milieu interstellaire, il y a correspondance avec les films de glaces d'eau amorphes dites de "haute densité" qui présentent un caractère poreux.

1.5 Formation de H₂ dans le milieu interstellaire

L'enchaînement de réactions chimiques en phase gazeuse peut conduire à la formation de molécules interstellaires, même les plus complexes. Cette production de molécules re-quiert cependant la présence de la molécule $H_2$, nécessaire au démarrage de toute cette chimie interstellaire active et complexe. Cette molécule est l'espèce la plus abondante dans l'univers et l'étude du mécanisme de sa formation encore mal caractérisé sera d'une
grande importance puisque cette molécule se forme dans des régions où les conditions de température de gaz ainsi que la nature et la température des grains sont très diverses.

### 1.5.1 Formation de H₂ en phase gazeuse

De nombreuses études ont démontré que la formation de H₂ en phase gazeuse par association radiative ou ionique est impossible [51, 9]. Lors d'une collision, deux atomes d'hydrogène dans l'état fondamental (1s) ont une énergie cinétique qui place le système à une énergie de 4.5 eV supérieure à l'état fondamental $ ^1\Sigma_g^+ $ de la molécule H₂. Il est au seuil de dissociation de la molécule et peut redonner deux atomes d'hydrogène. Pour évacuer cette énergie en excès, le système doit émettre un photon en effectuant une transition radiative vers un niveau d'énergie de rotation et de vibration stable (E<4.5 eV). Cependant, cette transition est trop improbable pour avoir lieu durant le temps de vie très court du système. D'autre part, si la molécule d'hydrogène est formée dans un état excité, elle pourrait libérer un photon en revenant à l'état fondamental. Cependant, le premier état excité de H₂ est $ ^3\Sigma_u^+ $, la transition vers le fondamental est donc interdite par le spin.

L'hydrogène moléculaire ne peut se former par un processus radiatif que si pendant que les deux premiers atomes d'hydrogène sont très proches, un troisième H est présent pour emporter l'énergie en effectuant une collision. Ceci correspond à une réaction à trois corps qui requiert une très forte densité qu'il est rare de trouver dans le MIS [52].

En phase gaz, H₂ peut se former selon les réactions [53] :

$$H + e \rightarrow H^{-} + h\nu$$

$$H^{-} + H \rightarrow H_{2} + e$$

La vitesse de la formation de H₂ dépend de la vitesse de la première réaction, qui est très faible par rapport à la deuxième, et de la concentration de $H^{-}$ formés qui peuvent soit former H₂, soit réagir avec $H^{+}$ ($H^{-} + H^{+} \rightarrow 2H$) ou avec un photon ($H^{-} + h\nu \rightarrow H + e$). L'hydrogène moléculaire peut aussi se former via le processus réactionnel suivant [53] :

$$H + H^{+} \rightarrow H_{2}^{+} + h\nu$$

$$H_{2}^{+} + H \rightarrow H_{2} + H^{+}$$

Cependant, le taux d'ionisation dans le MIS est trop faible pour que ces voies soient efficaces [53].

### 1.5.2 Formation de H₂ sur les surfaces de grains interstellaires

Le principe de formation de H₂ paraît très simple. Deux atomes d'hydrogène collés sur un grain se rencontrent et forment une molécule de H₂. En fait, le mécanisme est assez complexe et mal connu. De façon idéale, les atomes d'hydrogène qui se déplacent dans le MIS du fait de l'agitation thermique, entrent en collision avec un grain de poussière et se
collent à sa surface : c'est le phénomène d'adsorption. La probabilité de s'y coller dépend de la température du gaz, de celle du grain et de la nature chimique du grain qui influe sur l'énergie d'adsorption. Les deux atomes d'hydrogène ne restent pas immobiles mais peuvent être animés d'une grande mobilité sur cette surface, se déplaçant très rapidement d'un site à un autre. Plusieurs études expérimentales et théoriques ont été dédiées à la compréhension du mécanisme de formation de l'hydrogène moléculaire sur différents substrats, la distribution de l'énergie libérée et le niveau ro-vibrationnel dans lequel se forme la molécule H₂.

## Formation de H₂ sur la glace d'eau

Sur les surfaces de glaces d'eau typiques de celles qui recouvrent les poussières du milieu interstellaire dense et froid, plusieurs expériences [54, 55, 56, 19] ont été menées afin d'étudier la désorption de l'hydrogène moléculaire qui s'est révélée indispensable à l'interprétation des expériences de formation de l'hydrogène moléculaire sur les glaces d'eau. Le groupe de l'université de Cergy-pontoise [57, 54, 58, 59, 60, 61] a étudié cette réaction à l'aide de l'expérience FORMOLISM. Les distributions d'énergie d'adsorption des molécules H₂, HD et D₂ sur différents types de glaces d'eau ont été déterminées [57, 54]. Dans le cas des surfaces non poreuses, la totalité des molécules formées est rapidement libérée dans la phase gazeuse. Elles ont été détectées dans des états rovibrationnels excités jusqu'à $v=7$. Ces réactions ne sont efficaces que sur une faible gamme de température de surface. L'efficacité de formation chute en effet rapidement à 13 K, le temps de résidence des atomes sur la surface étant insuffisant au dessus de cette température. Pour les surfaces poreuses, les molécules formées ne sont plus observées directement dans la phase gazeuse comme dans le premier cas, mais uniquement par TPD (Temperature Programmed Desorption). La grande hétérogénéité de la surface permet difficilement aux molécules de désorber sans percuter à nouveau la surface.

La formation de l'hydrogène moléculaire sur la glace d'eau amorphe a été étudiée sur le plan théorique par des méthodes de dynamique moléculaire [62, 63]. Des méthodes de trajectoires classiques ont été aussi utilisées afin d'étudier l'adsorption et la diffusion de l'hydrogène sur la glace cristalline [64] ou amorphe [65]. Les probabilités d'adsorption calculées par ces différentes approches ne sont pas comparables. Ceci est lié à la morphologie de la surface de glace utilisée [65]. Les calculs montrent que l'atome d'hydrogène peut diffuser sur la surface de la glace avant de s'immobiliser [63, 65]. Les distances que peut parcourir cet atome, dépendent de l'énergie de l'atome incident. A partir d'une analyse des données expérimentales, Matar *et al.* [59] ont développé une expression analytique reliant l'énergie d'adsorption de H, D, H₂ et D₂ à l'énergie cinétique de l'hydrogène incident.

## Formation de H₂ sur le graphite, le carbone amorphe et les hydrocarbures aromatiques polycycliques (PAHs)

Expérimentalement, il a été longtemps considéré qu'en dehors de la physisorption à très basse température ($< 30$ K), les atomes d'hydrogène ne peuvent pas être chimisorbés sur une surface de graphite sans défauts cristallins [66]. En effet, la conversion de l'état
physisorbé à l'état chimisorbé nécessite un apport d'énergie estimée à 0.2 eV selon les travaux théoriques de Jeloaica et Sidis [67], de Ferro et al. [68] ainsi ceux de Ivanovskaya et al. [69]. Récemment, ceci a été confirmé expérimentalement [70]. Toutefois, les expériences de Zecho et al. [71, 72] ont montré que les atomes d'hydrogène et de deutérium ayant une énergie cinétique estimée à 2000 K (0.17 eV) peuvent être chimisorbés. Ainsi, plusieurs expériences [73, 74, 75, 76, 77] utilisent un faisceau d'atomes d'hydrogène à des températures allant de 1600 à 2200 K (0.14-0.19 eV) pour étudier la chimisorption des atomes d'hydrogène sur le graphite et le graphène. Dans le cas du carbone amorphe hydrogéné [78] et du HOPG (Highly Oriented Pyrolytic Graphite) [71], la molécule HD a été formée par le mécanisme Eley-Rideal (ER) dans lequel les atomes d'hydrogène du faisceau irradiant la surface réagissent directement avec les atomes de deutérium déjà adsorbés sur la surface. L'hydrogénation du graphène a aussi été étudiée [79, 80]. Les résultats obtenus ont mis en évidence le changement des paramètres de maille et des propriétés électroniques de la surface. Dans une étude d'interaction de deutérium avec le coronène (un hydrocarbure aromatique polycyclique) [81], les changements enregistrés dans les spectres IR ont été attribués au processus d'abstraction et de formation de HD et $D_2$. Une étude de superhydrogenation a également été menée [82]. L'irradiation avec les atomes d'hydrogène a conduit à leur adsorption sur le bord extérieur du coronène et ensuite à la formation de $H_2$ par un mécanisme ER. L'amélioration des méthodes DFT destinées à traiter les interactions atome-surface a conduit à leurs utilisations dans plusieurs études théoriques de formation de $H_2$ sur des surfaces de graphite [83, 84, 85, 86, 87]. Comme mentionné précédemment, il est crucial de connaître l'état rovibrationnel de la molécule désorbée. Plusieurs études ont été menées sur les états rovibrationnels de $H_2$ formés par le mécanisme Eley-Rideal toute en relaxant toute la surface de graphite [88, 89, 27] et en variant l'énergie de l'atome incident [90]. Dans le cas où deux atomes d'hydrogène sont physisorbés sur la surface du graphène, leurs recombinaison conduit à la formation de $H_2$ dans un état vibrationnel élevé ($v=8-12$) [88, 89, 27].

## Formation de $H_2$ sur des grains silicés

Des émissions en infra rouge ont été observées entre 11 et 11,5 µm et elles ont été attribuées aux grains de carbure de silicium (SiC) formés dans l'enveloppe circumstellaire de l'IRC +10216 [91, 92] et détectés dans le MIS [93]. C'est un des précurseurs chimiques des grains interstellaires qui proviennent majoritairement des résidus d'étoiles mortes ou en fin de vie.

Selon Walmsley et al. (1999) [94], les atomes de silicium réagissent avec le méthyle (CH₃) et/ou le carbène (CH₂) pour former un produit qui subira une série de photodissociations pour former finalement le carbure de silicium dans les régions de photodissociation (PDR). Cependant, ce scénario ne tient pas compte du rôle joué par les clusters hydrogénés de carbone et de silicium (SiCHₓ, x=6-3) qui ont été observés dans l'enveloppe circumstellaire de l'IRC +10216 [95, 96, 97, 98, 99, 100].

Kaiser et Osumara (2005) [95] ont proposé un autre modèle qui a été approuvé par la détection de ces mêmes clusters. La photolyse des molécules précurseurs chimiques du
silane et du méthane fournit les radicaux (CH₃) et (SiH₃). Ces derniers peuvent se com-
biner pour former le cluster SiCH₆ dans une réaction à trois corps où ce dernier emporte
l'énergie cinétique en excès. Si la zone où la réaction se déroule est éloignée du centre de
l'étoile, la molécule SiCH₆ ne peut pas être stabilisée par la collision à trois corps et elle
se décompose en SiCH₅ + H. Une succession de photolyses donnera à la fin le carbure de
silicium par élimination de l'hydrogène atomique et moléculaire. Le mécanisme proposé
par Kaiser et Osumara [95] met l'accent sur le rôle important que jouent ces clusters dans
la formation de SiC.

Dans ce travail, nous avons choisi d'étudier la réaction d'abstraction d'hydrogène à partir du tétraméthylsilane (TMS) par des atomes H, en tant que prototype pour l'étude de la formation de l'hydrogène moléculaire sur les grains silico-carbonés [101, 102, 103] des enveloppes d'étoiles riches en carbone. La molécule Si(CH₃)₄ pourrait être un précurseur des grains SiC. Elle a été utilisée pour modéliser des polytypes de structure cubique (βSiC) [104].

Formation de H₂ sur des grains silicatés

Les grains de poussière silicatés, en particulier, sont considérés comme de bons candidats pour catalyser la réaction H + H → H₂ dans le MIS. Les premières expériences de formation de H₂ et de ses analogues deutérés sur des olivines polycristallines ont été effectuées par Pirronello et ses collaborateurs [20, 105]. La formation de l'hydrogène moléculaire sur les silicates amorphes a été largement étudiée par le groupe de Vidali de l'université de Syracuse [17, 18, 106, 107] et par le groupe de Lemaire à l'université de Cergy-Pontoise [11, 108, 109, 110].

Jusqu'à aujourd'hui les caractéristiques des grains interstellaires sont encore mal connues. Néanmoins, on sait que les grains silicatés se trouvent sous la forme amorphe [111]. Le groupe de Vidali [17, 18, 106, 107] a étudié la formation d'hydrogène moléculaire dans la gamme de température (10-30 K) sur des surfaces de silicates amorphes ((Mg₁-x Fex)₂ SiO₄) de diverses compositions : le rapport Mg/Fe varie de 0 à 75%. Il s'est avéré que le processus de formation de la molécule H₂ est plus efficace sur les surfaces de silicates amorphes que sur les surfaces polycristallines. Ces mêmes travaux expérimentaux ont aussi trouvé que (i) les atomes d'hydrogène et de deutérium sont physisorbés sur les surfaces froides (≈10 K), (ii) quand la surface est réchauffée, une désorption de molécules HD est observée à une température qui dépend du type de la surface et (iii) le mécanisme impliqué est celui de Langmuir-Hinshelwood pour la gamme de températures allant de 10 à 20 K. Concernant le second point, il faudrait savoir que pour les expériences TPD, il est nécessaire de chauffer la surface pour rendre mobiles les atomes ou les molécules et activer ainsi la formation juste avant désorption. A partir des données expérimentales, un modèle [15, 17] permettant de remonter aux paramètres clés (énergie d'adsorption, barrière de mobilité, rétention des molécules formées par la surface) a été développé. Pour le silicate polycristallin, l'énergie de désorption de l'hydrogène est estimée à 32 meV [15] alors que pour la surface amorphe, elle est égale à 44 meV [17].

Les travaux de Lemaire et ses collaborateurs ont montré que la molécule D₂ formée sur
les silicates amorphes se trouve dans un état rovibrationnel élevé (v= 4, J= 2) [11]. Ils ont aussi mis en évidence deux mécanismes, Langmuir-Hinshelwood et Harris-Kasemo, qui sont en compétition à basses températures (<18 K) : une partie des molécules, formées (environ 33%) dans un état excité (v= 4, J= 2) par le mécanisme HK, se désorbe immédiatement dans la phase gazeuse alors que l'autre partie, formée par le mécanisme LH, ne parvient pas à quitter la surface avant le TPD.

L'effet de la présence des molécules de deutérium déjà adsorbées (pre-coverage) sur une surface de silicate amorphe sur la formation de D₂ a été étudié par Gavilan et al. [19]. Contrairement à ce qu'on attendait, ils ont trouvé que l'efficacité est augmentée. En effet, la diffusion des atomes D sur une surface avec des molécules D₂ déjà présentes est améliorée puisque les sites nécessitant une énergie d'adsorption élevée (et ainsi une plus grande barrière de diffusion) sont déjà occupés par D₂.

Cependant, tous ces résultats expérimentaux donnent seulement des tendances. Les principales causes sont les suivantes [111] : (i) les flux d'atomes d'hydrogène sont inévita-blement plus grands que dans le MIS, (ii) le coefficient de collage est une quantité qui est difficile à mesurer pour l'hydrogène atomique, (iii) dans le MIS, il y a une distribution des tailles de grains et (iv) la formation de l'hydrogène moléculaire dans l'espace se produit dans des conditions de régime approximativement stationnaire (par rapport à l'échelle de temps de plusieurs processus qui se produisent sur des surfaces, telles que le collage ou la désorption, qui sont beaucoup plus rapides), tandis que ce n'est pas toujours le cas dans les conditions du laboratoire.

De par leur structure complexe, les mécanismes réactionnels mettant en jeu ce type de surfaces ont été peu étudiés théoriquement jusqu'ici. La surface silicatée a été modé-lisée principalement par des nano-clusters [112, 113] ou par des approches périodiques [29, 28, 30, 14] utilisant une surface cristalline ayant la stoechiométrie de la forstérite (010). Une étude théorique récente [29] a mis en évidence l'existence d'une possible voie de formation de H₂ par chimisorption sans barrière sur une surface de pristine olivine Mg₂SiO₄ [forstérite (010)], susceptible de produire des molécules ro-vibrationnellement excitées. Ce résultat est intéressant, puisqu'il a permis d'établir un itinéraire à la forma-tion de H₂ à partir de deux hydrogènes chemisorbés à température relativement élevée où H<sub>phys</sub> est instable (nuage diffus, T=100K). Seul l'aspect énergétique a été considéré dans cette étude. Garcia-Gil et al. [28] ont également exploré les différents sites d'adsorption de l'hydrogène sur la même surface de forstérite (010) au moyen de calculs périodiques moyennant des méthodes DFT (PBE) et des bases de valences double et triple zêta implémentées dans le programme SIESTA [114]. Selon ces deux études, le site de la physi-sorption le plus favorable est sur l'atome de Mg et celui de la chimisorption est sur l'atome d'oxygène voisin. Ces mêmes résultats ont été trouvés dans le travail de Navarro-Ruiz et al. [14] qui a été mené sur la même surface. Ils ont trouvé que la diffusion de l'hydrogène est limitée puisque les constantes de vitesse de diffusion calculées avec la théorie de l'état de transition (TST) corrigées par l'approche de Wigner, sont négligeables à 10 K. Elles ne deviennent significatives qu'à partir de 100 K, température à laquelle les atomes seront désorbés. Ils ont conclu que l'hydrogène ne peut se déplacer sur le grain qu'à travers
l'adsorption et la désorption. La chemisorption est aussi entravée à 10 K à cause de la
barrière à la chemisorption directe de l'hydrogène. La formation de la molécule H₂ se
fait selon le mécanisme de Langmuir-Hinshelwood à partir de deux atomes d'hydrogène
initialement physisorbés avec une barrière cinétique quasiment négligeable, ou à partir de
deux atomes d'hydrogène initialement chimisorbés avec une barrière surmontable à des
températures bien supérieures à 10 K.

En plus de la forstérite, un autre travail a été mené sur la fayalite Fe₂SiO₄ (010) en
utilisant le programme VASP [30]. Les calculs ont montré qu'il n'y a pas de barrière
d'activation pour l'adsorption de l'hydrogène atomique sur ces surfaces. Cependant il y
a formation soit d'un hydrure (adsorption sur le magnésium) soit d'un proton suite à la
formation des groupes hydroxyle (OH).

Des études de formation de H₂ [112, 113] ont été menées sur des grains de silicate
modélisés respectivement par des petits clusters d'olivine ayant la stœchiométrie de la
forstérite (MgO)₆(SiO₂)₃ et par des petits clusters de pyroxène ayant la stœchiométrie
de l'enstatite (MgO)₄(SiO₂)₄. Ces nanoclusters sont des espèces non recouvertes par la
glace d'eau ce qui rend ces modèles plus pertinents pour modéliser les grains de poussière
dans les nuages diffus. Ces grains comprennent environ 10% de nanoparticules amorphes
et cristallines avec des diamètres ≤ 15 Å [115]. Ces nanosilicates sont des intermédiaires
fondamentaux dans la formation de poussières de silicates des sorties stellaires, et sont
omniprésents partout dans le MIS. Les travaux de Kerkeni et Bromley [113] ont permis
d'identifier des sites de physisorption et de chimisorption d'atomes d'hydrogène, d'étu-
dier la diffusion de l'atome H physisorbé sur un cluster d'olivine de type foresterite. Les
caractéristiques énergétiques de ces sites (les énergies de chimisorption, de physisorption
et d'activation pour la désorption et la recombinaison de H₂) ont été déterminées. Le
chemin réactionnel le plus énergétiquement favorable pour la formation de H₂ à partir de
deux H<sub>chem</sub> possède une barrière à la recombinaison de 6 kJ/mol et une exothermicité
de -86 kJ/mol. D'autres profils énergétiques sont favorables pour la dissociation de H₂.
Cependant, chacun de ces chemins possède une assez grande barrière (50-55 kJ/mol) qui
aurait tendance à inhiber le processus de dissociation. Le chemin réactionnel (exother-
mique) trouvé sur la forstérite (010) par Goumans et al. [29] est localisé entre les voies
exothermiques et endothermiques pour la formation de H₂ sur les nanosilicates. Compa-
rés aux surfaces cristallines, les nanosilicates offrent une plus grande variété de profils
énergétiques pour la formation et la dissociation de H₂.
# Chapitre 2

## Méthodes de chimie quantique

### 2.1 Introduction

La chimie quantique consiste en l'application de la mécanique quantique à des systèmes polyélectroniques afin d'en extraire les différentes propriétés structurales (géométries), spectroscopiques (fréquence de vibration), énergétiques (énergies de liaison, énergies d'excitation etc...) et réactionnelles (barrière d'activation, liaisons chimiques etc...). A partir de ces quantités, des surfaces d'énergie potentielle précises sont construites. Ces surfaces servent à modéliser et à étudier en particulier la dynamique des réactions chimiques.

Dans ce chapitre, les méthodes théoriques nécessaires au traitement des problèmes électroniques sont présentées. Après avoir introduit l'approximation de Born-Oppenheimer qui permet de réduire la complexité de la résolution de l'équation de Schrödinger, nous décrivons en détails la méthode Hartree-Fock qui est à la base des méthodes utilisées dans ce travail. Nous décrivons ensuite les méthodes *ab initio* d'interaction de configuration (MCSCF, MRCI et Cluster Couplés) de plus en plus sophistiquées, qui permettent de tenir compte de l'énergie de corrélation et de bien décrire les énergies électroniques. Nous introduisons ensuite la méthode MP2 basée sur une approximation perturbative. La théorie de la fonctionnelle de la densité (DFT) et les bases atomiques sont aussi présentées.

Nous avons employé ces méthodes pour décrire la réaction d'abstraction de l'hydrogène atomique du tétraméthylsilane par collision avec un atome d'hydrogène et un radical hydroxyle ainsi que l'interaction de l'hydrogène atomique avec les nanograins de silicates et la formation de l'hydrogène moléculaire sur ces nanograins.

### 2.2 Equation de Schrödinger

L'équation de Schrödinger indépendante du temps s'écrit :

$$ \hat{H}\Psi(R, r) = E\Psi(R, r) \qquad (2.1) $$

où $\hat{H}$ est l'opérateur hamiltonien associé au système, E est l'énergie de l'état stationnaire du système. r et R désignent respectivement, l'ensemble des coordonnées des électrons et
des noyaux dans le référentiel du centre de masse des noyaux.

On considère un système moléculaire de n électrons (indexés par i) de masse $m_e$ et de N noyaux (indexés par $\alpha$) de charge $Z_\alpha$ et de masse $M_\alpha$. L'hamiltonien non relativiste s'écrit comme suit :

$$ \hat{H} = \hat{T}_e + \hat{T}_N + \hat{V}_{ee} + \hat{V}_{eN} + \hat{V}_{NN} \quad (2.2) $$

Les opérateurs d'énergie cinétique et potentielle s'écrivent :

$$ \hat{T}_e = - \sum_{i=1}^{n} \frac{\hbar^2}{2m_e} \nabla_i^2 \quad (2.3) $$

$$ \hat{T}_N = - \sum_{\alpha=1}^{N} \frac{\hbar^2}{2M_\alpha} \nabla_\alpha^2 \quad (2.4) $$

$$ \hat{V}_{eN} = - \sum_{i=1}^{n} \sum_{\alpha=1}^{N} \frac{Z_\alpha e^2}{4\pi\epsilon_0 |r_i - R_\alpha|} \quad (2.5) $$

$$ \hat{V}_{ee} = \sum_{i=1}^{n} \sum_{j>i}^{n} \frac{e^2}{4\pi\epsilon_0 |r_i - r_j|} \quad (2.6) $$

$$ \hat{V}_{NN} = \sum_{\alpha=1}^{N} \sum_{\beta>\alpha}^{N} \frac{Z_{\alpha} Z_{\beta} e^2}{4\pi\epsilon_0 |R_{\alpha} - R_{\beta}|} \quad (2.7) $$

Les deux premiers termes de l'hamiltonien sont des opérateurs d'énergie cinétique des n électrons et des N noyaux alors que les trois derniers termes sont respectivement les potentiels d'interaction électron-noyau, électron-électron et noyau-noyau.

## 2.3 Approximation de Born-Oppenheimer

L'équation de Schrödinger (eq. 2.1) devient très complexe ce qui rend sa résolution analytique pour un système polyélectronique (n>2) impossible. L'approximation de Born Oppenheimer [116] est utilisée dans le cadre moléculaire pour simplifier cette équation en traitant les mouvements électroniques et nucléaires séparément. En effet, la masse d'un noyau est beaucoup plus grande que la masse d'un électron (mp/me ≈ 1836). Il s'en suit que le mouvement des électrons est beaucoup plus rapide que celui des noyaux. On peut ainsi en première approximation, pour chaque configuration fixe des noyaux, calculer une fonction d'onde purement électronique, décrivant le mouvement des électrons.

Suite à cette approximation :

- l'hamiltonien $\hat{H}$ se décompose en une partie purement nucléaire et une partie incluant toutes les contributions électroniques :
$$ \hat{H} = \hat{T}_N + \hat{H}_e \qquad (2.8) $$

$$ \hat{H}_e = \hat{T}_e + \hat{V}_{ee} + \hat{V}_{eN} + \hat{V}_{NN} \qquad (2.9) $$

Puisque l'hamiltonien électronique $\hat{H}_e$ s'exprime en fonction de $\hat{V}_{eN}$, il dépend paramétriequement des coordonnées nucléaires $R_α$. Les vecteurs propres $\psi_{ei}$ de $\hat{H}_e$ en dépendent également. L'équation de Schrödinger électronique s'écrit :

$$ \hat{H}_e \psi_{ei}(R_\alpha, r_i) = E_{ei}(R_\alpha) \psi_{ei}(R_\alpha, r_i) $$

Pour simplifier l'écriture des équations, les unités atomiques sont utilisées :

$$ \left[ -\sum_{i=1}^{n} \frac{1}{2} \nabla_i^2 - \sum_{i=1}^{n} \sum_{\alpha=1}^{N} \frac{Z_\alpha}{|r_i - R_\alpha|} + \sum_{i=1}^{n} \sum_{j>i} \frac{1}{|r_i - r_j|} + V_{NN} \right] \psi_{ei}(R_\alpha, r_i) = E_{ei}(R_\alpha) \psi_{ei}(R_\alpha, r_i) \qquad (2.10) $$

$E_{ei}(R_α)$ est l'énergie associée à l'état électronique i pour la configuration $R_α$. Le terme d'interaction entre les noyaux $V_{NN}$ est constant puisqu'il ne dépend pas de la position des électrons.

- La fonction d'onde totale $\Psi$ du système s'écrit sous la forme d'une combinaison linéaire de produits de fonctions d'onde nucléaires $\psi_{Ni}$ et de fonctions d'onde électroniques $\psi_{ei}$ dépendant paramétriquement des positions des noyaux :

$$ \Psi(R_\alpha, r_i) = \sum_i \psi_{Ni}(R_\alpha) \psi_{ei}(R_\alpha, r_i) \qquad (2.11) $$

En remplaçant dans l'équation de Schrödinger (2.1) l'expression 2.11 de $\Psi$, on obtient :

$$ \begin{aligned} \hat{H} \Psi(R_\alpha, r_i) = & E \Psi(R_\alpha, r_i) \\ & \sum_i (\hat{T}_N + \hat{H}_e(R_\alpha)) \psi_{Ni}(R_\alpha) \psi_{ei}(R_\alpha, r_i) = E \sum_i \psi_{Ni}(R_\alpha) \psi_{ei}(R_\alpha, r_i) \\ & \sum_i \left[ \sum_\alpha \left( -\frac{1}{2M_\alpha} \nabla_\alpha^2 + E_{ei}(R_\alpha) \right) \psi_{Ni}(R_\alpha) \psi_{ei}(R_\alpha, r_i) \right. \\ & \left. + \sum_i \psi_{Ni}(R_\alpha) \sum_\alpha \left( -\frac{1}{2M_\alpha} \nabla_\alpha^2 \psi_{ei}(R_\alpha, r_i) \right) \right. \\ & \left. - \sum_i \sum_\alpha \left( -\frac{1}{2M_\alpha} \left[ \nabla_\alpha \psi_{Ni}(R_\alpha) \right] \left[ \nabla_\alpha \psi_{ei}(R_\alpha, r_i) \right] \right) \right] = E \sum_i \psi_{Ni}(R_\alpha) \psi_{ei}(R_\alpha, r_i) \end{aligned} \qquad (2.12) $$

En utilisant l'orthonormalité des fonctions d'onde électroniques et en multipliant par $\psi_{ek}^*$
et en intégrant sur les coordonnées électroniques, on a :

$$ \left(\sum_{\alpha} -\frac{1}{2M_{\alpha}}\nabla_{\alpha}^{2} + E_{ei}(R_{\alpha})\right)\psi_{Ni}(R_{\alpha}) \right. \\ \left. + \sum_{i}\sum_{\alpha}\langle\psi_{ek}(R_{\alpha}, r_{i})| - \frac{1}{2M_{\alpha}}\nabla_{\alpha}^{2}|\psi_{ei}(R_{\alpha}, r_{i})\rangle\psi_{Ni}(R_{\alpha}) \right. \\ \left. - \sum_{i}\sum_{\alpha} -\frac{1}{2M_{\alpha}}\nabla_{\alpha}\psi_{Ni}(R_{\alpha})\langle\psi_{ek}(R_{\alpha}, r_{i})|\nabla_{\alpha}|\psi_{ei}(R_{\alpha}, r_{i})\rangle = E\psi_{Ni}(R_{\alpha}) \right) \quad (2.13) $$

Les différents états électroniques sont couplés dynamiquement via la matrice des dérivées premières du couplage cinétique ($\langle\psi_{ek}(R_{\alpha}, r_{i})|\nabla_{\alpha}|\psi_{ei}(R_{\alpha}, r_{i})\rangle$) et la matrice des dérivées secondes du couplage cinétique ($\langle\psi_{ek}(R_{\alpha}, r_{i})| - \frac{1}{2M_{\alpha}}\nabla_{\alpha}^{2}|\psi_{ei}(R_{\alpha}, r_{i})\rangle$).

### 2.3.1 Approximation adiabatique

Dans l'approximation adiabatique, les termes de couplage dynamique non diagonaux sont négligeables puisqu'ils sont proportionnels à $\frac{1}{M_{\alpha}}$. De plus, on néglige également le terme diagonal $\langle\psi_{ek}(R_{\alpha}, r_{i})|\nabla_{\alpha}|\psi_{ei}(R_{\alpha}, r_{i})\rangle$ qui est très petit et totalement négligeable pour des problèmes de collision. Les électrons très mobiles suivent adiabatiquement le mouvement des noyaux. Les fonctions d'onde $\psi_{ei}$ varient donc peu avec $R_{\alpha}$ et par conséquent, les dérivées des états propres de $\hat{H}_{e}$ par rapport à $R_{\alpha}$ sont négligeables. Pour l'état électronique i, l'équation (2.13) se réduit alors à :

$$ [\hat{T}_N + E_{ei}(R_\alpha)]\psi_{Ni}(R_\alpha) = E\psi_{Ni}(R_\alpha) \quad (2.14) $$

On remarque que la valeur propre de l'Hamiltonien électronique s'interprète comme l'énergie potentielle du mouvement nucléaire. Les noyaux se meuvent dans le potentiel $E_{ei}(R_{\alpha})$ créé par les électrons et les noyaux.

L'approximation adiabatique n'est valable que pour des états électroniques énergétiquement bien séparés. Elle ne l'est plus dans le cas d'un croisement évité de deux surfaces d'énergie potentielle où les deux états électroniques se mélangent et les fonctions d'onde électroniques varient fortement avec les distances internucléaires.

Une fois l'équation de Schrödinger électronique (éq. 2.10) résolue, il est possible de résoudre le problème nucléaire. En calculant $E_{ej}(R_{\alpha})$ pour un ensemble de configurations nucléaires ($R_{\alpha}$), une surface d'énergie potentielle (SEP) est obtenue. La résolution de l'équation de Schrödinger nucléaire (équation 2.14) permet de déterminer les mouvements de vibration et de rotation d'une molécule ou les mouvements des noyaux lors d'une collision.

## 2.4 Fonction d'onde électronique

Afin de résoudre l'équation de Schrödinger électronique, nous avons recours à des approximations qui rendent le calcul plus simple. Mais tout d'abord, la forme générale
de la fonction d'onde électronique doit être décrite. Elle doit inclure le spin de chaque électron sous la forme de spin-orbitales et elle doit satisfaire le principe de Pauli pour un système de n électrons.

### 2.4.1 Spin-orbitales

Un électron est décrit au moyen d'une orbitale spatiale $\phi(\vec{r})$, fonction du vecteur position $\vec{r}$. La probabilité de trouver l'électron dans le volume élémentaire $d\vec{r}$ est $|\phi(\vec{r})|^2 d\vec{r}$. Les orbitales spatiales forment un ensemble orthonormé tel que :

$$ \int dr \phi_i^*(\vec{r}) \phi_j(\vec{r}) = \delta_{ij} \qquad (2.15) $$

Cependant, pour décrire complètement un électron, nous avons besoin de préciser également son spin. Deux fonctions de spin, $|\alpha\rangle$ et $|\beta\rangle$, orthogonales entre elles, sont alors introduites :

$$ \begin{aligned} \langle\alpha|\alpha\rangle &= \langle\beta|\beta\rangle = 1 \\ \langle\alpha|\beta\rangle &= \langle\beta|\alpha\rangle = 0 \end{aligned} \qquad (2.16) $$

Ainsi, il existe deux spin-orbitales associées à une orbitale spatiale :

$$ \Phi(\vec{q}_i) = \phi(\vec{r}_i)\alpha \qquad (2.17) $$

$$ \Phi(\vec{q}_i) = \phi(\vec{r}_i)\beta \qquad (2.18) $$

où $\vec{q}_i$ représente les coordonnées spatiales et de spin de l'électron i. L'ensemble des spin-orbitales forme une base orthonormée :

$$ \int dr \Phi_i^*(\vec{q}_i) \Phi_j(\vec{q}_i) = \langle \Phi_i | \Phi_j \rangle = \delta_{ij} \qquad (2.19) $$

La partie spatiale des orbitales moléculaires (OM) peut être prise sous la forme d'une combinaison linéaire d'orbitales atomiques (LCAO). Cette approximation consiste à développer les orbitales moléculaires sur la base des orbitales atomiques $\chi_\mu$ des différents atomes formant la molécule étudiée :

$$ \Phi_i = \sum_{\mu=1}^{M} C_{\mu i} \chi_{\mu} \qquad (2.20) $$

où M est le nombre d'orbitales atomiques utilisé pour représenter une orbitale moléculaire. $C_{\mu i}$ sont les coefficients de développement de l'orbitale moléculaire i sur la base orbitales atomiques $\chi_\mu$.
Dans le cadre de l'approximation orbitale, la fonction d'onde d'un système à n électrons est représentée sous la forme d'un produit de spin-orbitales :

$$ \psi_e(\vec{q}_1, \vec{q}_2, \dots, \vec{q}_n) = \Phi_1(\vec{q}_1)\Phi_2(\vec{q}_2)\dots\Phi_n(\vec{q}_n) \qquad (2.21) $$

Pour respecter le principe d'exclusion de Pauli, la fonction d'onde d'un système à n électrons doit être antisymétrique par rapport à l'échange de deux électrons. Pour cela, nous utilisons les déterminants de Slater.

## 2.4.2 Déterminants de Slater

La fonction d'onde électronique, définie dans le cadre de l'approximation orbitale comme un produit de spin-orbitales, ne permet pas de respecter le principe de Pauli. Ce principe est une condition à laquelle est soumise une fonction d'onde multiélectronique afin de satisfaire l'indiscernabilité des électrons :

Dans un système formé de n particules quantiques identiques, il n'existe pas de moyen de distinguer et de suivre l'un d'entre eux. Les deux particules peuvent échanger leurs positions, mais ceci laisse toutes les propriétés observables du système inchangées. La fonction d'onde électronique de ce système $\psi_e(1,2...i...j...n)$ doit satisfaire les postulats de la mécanique quantique en imposant des conditions sur la fonction d'onde quant à son comportement lors de l'échange des coordonnées des particules. Puisque les propriétés physiques ne peuvent pas dépendre de l'identité des deux particules, elles doivent rester invariantes à l'échange des coordonnées i et j. Entre autres, la densité de probabilité $|\psi_e|^2$ doit être inchangée en tout point selon :

$$ |\psi_e(1, 2\dots i\dots j\dots n)|^2 = |\psi_e(1, 2\dots j\dots i\dots n)|^2 \qquad (2.22) $$

Pour les électrons, la fonction d'onde est antisymétrique :

$$ \psi_e(...i...j...) = -\psi_e(...j..i...) \qquad (2.23) $$

Une fonction d'onde multiélectronique doit donc être antisymétrique par rapport à la permutation de deux électrons.

Le principe de Pauli a un effet énorme sur la chimie, en effet, il est relié aux propriétés magnétiques des électrons. La fonction d'onde électronique $\psi_e$ peut être exprimée sous la forme d'un produit antisymétrisé de spin-orbitales, déterminant de Slater :

$$ \psi_e = \frac{1}{\sqrt{n!}} \begin{vmatrix} \Phi_1(1) & \Phi_2(1) & \cdots & \Phi_n(1) \\ \Phi_1(2) & \Phi_2(2) & \cdots & \Phi_n(2) \\ \vdots & \vdots & \ddots & \vdots \\ \Phi_1(n) & \Phi_2(n) & \cdots & \Phi_n(n) \end{vmatrix} \qquad (2.24) $$

Les déterminants de Slater peuvent aussi être notés de manière abrégée :

$$ \psi_e = \hat{A} |\Phi_1 \Phi_2 \cdots \Phi_n| \qquad (2.25) $$

$$ \hat{A} = \frac{1}{N!} \sum_{p=0}^{N-1} (-1)^p \hat{P} = \frac{1}{N!} [\hat{1} - \sum_{ij} \hat{P}_{ij} + \sum_{ijk} \hat{P}_{ijk} - \dots] \qquad (2.26) $$
où $\hat{A}$ est l'opérateur d'antisymétrisation, $\hat{1}$ est l'opérateur identité, $P_{ij}$ génère toutes les possibilités de permutation de 2 électrons et $P_{ijk}$ génère toutes les possibilités de permutation de 3 électrons, etc.

## 2.5 Propriétés de l'énergie électronique

### 2.5.1 Principe variationnel

Le principe variationnel est à la base de tous les calculs modernes de structure mo- léculaire. Il se base sur le choix d'une fonction d'onde d'essai afin de calculer l'énergie électronique. Cette fonction d'onde approchée dépend d'un certain nombre de paramètres ajustables appelés paramètres variationnels. Ces coefficients sont variés de sorte à avoir une valeur minimale de l'énergie. Cette dernière est supérieure à l'énergie exacte (la valeur propre exacte de l'Hamiltonien électronique).

## 2.6 Méthode Hartree Fock

La méthode Hartree-Fock [117, 118, 119] permet d'obtenir une solution de l'équation de Schrödinger électronique pour des systèmes et des atomes polyélectroniques en utilisant une résolution variationnelle dans laquelle la fonction d'onde approchée est écrite sous la forme d'un déterminant de Slater. L'énergie HF est une première approximation de l'énergie exacte pour l'état fondamental du système considéré.

### 2.6.1 Energie Hartree-Fock

On cherche à déterminer l'énergie électronique, $E_{HF}$, qui peut s'écrire sous la forme :

$$E_{HF} = \langle \psi_e | \hat{H}_e | \psi_e \rangle \qquad (2.27)$$

où

$$\hat{H}_e = \sum_{i=1}^{n} \hat{h}(i) + \sum_{i=1}^{n} \sum_{j>i}^{n} \hat{g}(i, j) + \hat{V}_{NN} \qquad (2.28)$$

$\psi_e$ est un déterminant de Slater (équation 2.24).

A partir de l'expression (2.10), on peut définir un opérateur de coeur $\hat{h}(i)$ et un opérateur biélectronique $\hat{g}(i, j)$ :

$$\hat{h}(i) = -\frac{1}{2} \nabla_i^2 - \sum_{\alpha=1}^{N} \frac{Z_\alpha}{r_{i\alpha}} \qquad (2.29)$$

$$\hat{g}(i, j) = \frac{1}{r_{ij}} \qquad (2.30)$$

$$où \frac{1}{r_{i\alpha}} = \frac{1}{|r_i - R_{i\alpha}|} \text{ et } \frac{1}{r_{ij}} = \frac{1}{|r_i - r_{ij}|}$$
L'énergie électronique relative à l'opérateur de l'équation (2.28) peut être écrite en fonction d'intégrales à un et deux électrons :

$$E_{HF} = \sum_{i=1}^{n} h_i + \sum_{i=1}^{n} \sum_{j>i}^{n} (J_{ij} - K_{ij}) + V_{NN} \quad (2.31)$$

où

$$h_i = \langle \Phi_i(1) | \hat{h}(1) | \Phi_i(1) \rangle \quad (2.32)$$

$$J_{ij} = \langle \Phi_i(1) \Phi_j(2) | \hat{g}(1, 2) | \Phi_i(1) \Phi_j(2) \rangle \quad (2.33)$$

$$K_{ij} = \langle \Phi_i(1) \Phi_j(2) | \hat{g}(1, 2) | \Phi_j(1) \Phi_i(2) \rangle \quad (2.34)$$

L'intégrale de Coulomb $J_{ij}$ représente l'énergie de répulsion des électrons indépendemment de leurs spins. L'intégrale d'Echange $K_{ij}$ est la conséquence de la forme antisymétrique du déterminant de Slater. L'opérateur de répulsion nucléaire est indépendant des coordonnées des électrons et peut être intégré pour obtenir une constante telle que :

$$\langle \psi_e | \hat{V}_{NN} | \psi_e \rangle = V_{NN} \langle \psi_e | \psi_e \rangle = V_{NN} \quad (2.35)$$

L'expression de l'énergie 2.31 peut s'écrire aussi sous une forme plus symétrique :

$$E_{HF} = \sum_{i=1}^{n} h_i + \frac{1}{2} \sum_{i=1}^{n} \sum_{j=i}^{n} (J_{ij} - K_{ij}) + V_{NN} \quad (2.36)$$

Le facteur $\frac{1}{2}$ tient compte de la double somme sur tous les électrons.

$E_{HF}$ peut être exprimée en termes d'opérateur de Coulomb $\hat{J}(i)$ et d'Echange $\hat{K}(i)$ :

$$E_{HF} = \sum_{i}^{n} \langle \Phi_i | \hat{h}(i) | \Phi_i \rangle + \frac{1}{2} \sum_{i,j}^{n} [ \langle \Phi_j | \hat{J}(i) | \Phi_j \rangle - \langle \Phi_j | \hat{K}(i) | \Phi_j \rangle ] + V_{NN} \quad (2.37)$$

où

$$\hat{J}(i) |\Phi_j(2)\rangle = \langle \Phi_i(1) | \hat{g}(1, 2) | \Phi_i(1) \rangle |\Phi_j(2)\rangle \quad (2.38)$$

$$\hat{K}(i) |\Phi_j(2)\rangle = \langle \Phi_i(1) | \hat{g}(1, 2) | \Phi_j(1) \rangle |\Phi_i(2)\rangle \quad (2.39)$$

## 2.6.2 Méthode du "champ autocohérent", SCF

Afin de déterminer l'énergie Hartree-Fock, il faut minimiser $E_{HF}$ en faisant varier les spin-orbitales $\Phi_i$ selon le principe variationnel. Cette optimisation se fait en respectant les contraintes d'orthonormalisation ($\langle \Phi_i | \Phi_j \rangle = \delta_{ij}$) :

$$\delta E_{HF} - \sum_{i,j}^{n} \lambda_{i,j} (\langle \delta \Phi_i | \Phi_j \rangle - \langle \Phi_i | \delta \Phi_j \rangle) = 0 \quad (2.40)$$
où les $\lambda_{ij}$ sont des multiplicateurs de Lagrange [120].

A partir de l'équation (2.37), on déduit :

$$ \delta E_{HF} = \sum_{i}^{n} [\langle \delta\Phi_i | \hat{h}(i) | \Phi_i \rangle + \langle \Phi_i | \hat{h}(i) | \delta\Phi_i \rangle] + \frac{1}{2} \sum_{i,j}^{n} [\langle \delta\Phi_i | \hat{J}(j) - \hat{K}(j) | \Phi_i \rangle \\ + \langle \Phi_i | \hat{J}(j) - \hat{K}(j) | \delta\Phi_i \rangle + \langle \delta\Phi_i | \hat{J}(i) - \hat{K}(i) | \Phi_j \rangle + \langle \Phi_j | \hat{J}(i) - \hat{K}(i) | \delta\Phi_j \rangle] \quad (2.41) $$

Les troisième et cinquième termes sont identiques ainsi que les quatrième et sixième termes. Ces termes peuvent être collectés ensemble :

$$ \delta E_{HF} = \sum_{i}^{n} [\langle \delta\Phi_i | \hat{h}(i) | \Phi_i \rangle + \langle \Phi_i | \hat{h}(i) | \delta\Phi_i \rangle] + \sum_{i,j}^{n} [\langle \delta\Phi_i | \hat{J}(j) - \hat{K}(j) | \Phi_i \rangle + \langle \Phi_i | \hat{J}(j) - \hat{K}(j) | \delta\Phi_i \rangle] \quad (2.42) $$

$$ \delta E_{HF} = \sum_{i}^{n} [\langle \delta\Phi_i | \hat{F}(i) | \Phi_i \rangle + \langle \Phi_i | \hat{F}(i) | \delta\Phi_i \rangle] \quad (2.43) $$

On définit alors l'opérateur monoélectronique de Fock $\hat{F}_i$ :

$$ \hat{F}(i) = \hat{h}(i) + \sum_{j}^{n} (\hat{J}(j) - \hat{K}(j)) \quad (2.44) $$

L'opérateur de Fock est un opérateur monoélectronique effectif qui décrit l'énergie ci-nétique d'un électron, le potentiel d'attraction électron-noyau ainsi que le potentiel de répulsion des autres électrons à travers les opérateurs de Coulomb $\hat{J}(j)$ et d'Echange $\hat{K}(j)$. Il faut noter que $\hat{F}(i)$ est lié à $\delta E_{HF}$ et non $E_{HF}$ et que sa somme ne correspond pas à l'hamiltonien électronique $\hat{H}_e$.

Puisque les coefficients de Lagrange sont hermitiques, l'équation 2.40 devient :

$$ \sum_{i}^{n} \delta \langle \Phi_i | \hat{F}(i) | \Phi_i \rangle - \sum_{i,j}^{n} \lambda_{ij} \delta \langle \Phi_i | \Phi_j \rangle = 0 \quad (2.45) $$

On traite les variations $\langle \delta\Phi_i |$ et $|\delta\Phi_j \rangle$ comme indépendantes. On projette alors la relation ci-dessus sur $\langle \delta\Phi_i |$ :

$$ \sum_{i}^{n} \langle \delta\Phi_i | \hat{F}(i) | \Phi_i \rangle - \sum_{i,j}^{n} \lambda_{ij} \langle \delta\Phi_i | \Phi_j \rangle = 0 \quad (2.46) $$

On obtient ainsi les équations de HF :

$$ \hat{F}(i)\Phi_i = \sum_{j}^{n} \lambda_{ij}\Phi_j \quad (2.47) $$
Ces équations peuvent être simplifiées en utilisant les orbitales moléculaires canoniques $\Phi'_i$ associées aux valeurs propres $\epsilon'_i$ telles que :

$$ \hat{F}(i)\Phi'_{i} = \epsilon'_{i}\Phi'_{i} \qquad (2.48) $$

Les $\Phi'_i$ sont obtenues à partir de la diagonalisation de la matrice des multiplicateurs de Lagrange ($\lambda_{ij}=0$ pour $i \neq j$, $\lambda_{ii}=\epsilon_i$).

On reprend la notation $\Phi_i$ pour les orbitales moléculaires canoniques et $\epsilon_i$ pour les énergies correspondantes. Les spin-orbitales $\Phi_i$ correspondent à des combinaisons linéaires de M orbitales atomiques $\chi_\mu$ appelées fonctions de base :

$$ \Phi_i = \sum_{\mu=1}^{M} C_{i\mu}\chi_{\mu} \qquad (2.49) $$

On peut alors réécrire l'équation 2.47 selon l'expression suivante :

$$ \hat{F}_i \sum_{\mu=1}^{M} C_{i\mu}\chi_{\mu} = \epsilon_i \sum_{\mu=1}^{M} C_{i\mu}\chi_{\mu} \qquad (2.50) $$

En tenant compte de ce développement sur la base des fonctions $\chi_\mu$, le système d'équations de Fock résultant de l'application des variations linéaires comporte M équations. L'équation 2.50 peut s'écrire sous une forme matricielle. On obtient alors les **équations de Roothaan** [121] :

$$ FC = \epsilon SC \qquad (2.51) $$

C comprend les coefficients $C_{i\mu}$ qui sont les paramètres variationnels devant être optimisés.

$\epsilon$ est une matrice dont les éléments diagonaux sont les énergies $\epsilon_i$ des spin-orbitales moléculaires $\Phi_i$.

La matrice $F$ est représentative de l'opérateur de Fock dans la base d'orbitales atomiques $\chi_\mu$ et contient les intégrales de Fock :

$$ F_{\mu\lambda} = \langle \chi_{\mu} | \hat{F} | \chi_{\lambda} \rangle \qquad (2.52) $$

$S$ est la matrice des recouvrements des orbitales atomiques ayant pour éléments :

$$ S_{\mu\lambda} = \langle \chi_{\mu} | \chi_{\lambda} \rangle \qquad (2.53) $$

La résolution de l'équation 2.51 se fait d'une façon itérative. Il faut :

1/ partir d'un ensemble initial d'orbitales moléculaires (combinaisons linéaires des fonctions de base),

2/ construire l'opérateur de Fock correspondant,

3/ résoudre les équations 2.48 pour obtenir un nouvel ensemble d'orbitales,

4/ construire un nouvel opérateur de Fock,
5/ répéter la procédure jusqu'à convergence.

L'expression de $E_{HF}$ s'écrit en fonction des intégrales de Coulomb $J_{ij}$ (2.33) et d'Echange $K_{ij}$ (2.34) :

$$E_{HF} = \sum_{i}^{n} \epsilon_{i} - \frac{1}{2} \sum_{ij} [J_{ij} - K_{ij}] + V_{NN} \quad (2.54)$$

$$\epsilon_i = \langle \Phi_i | \hat{F}(i) | \Phi_i \rangle \quad (2.55)$$

Cette énergie Hartree-Fock s'approche de la valeur exacte de l'énergie électronique par valeurs supérieures, jusqu'à une limite fixée par la taille de la base.

Il faut noter que la résolution des équations Hartree-Fock se fait de deux façons selon que le système considéré est :

- à couche fermée (les électrons sont appariés) : c'est le formalisme RHF (restricted Hartree-Fock) [121]. Les spin-orbitales sont décrites à l'aide d'un ensemble de fonctions de base {$\chi_{\lambda}$} ($\lambda$=1,M). Un système matriciel (équation 2.51) connu sous le nom d'équation de Roothaan est alors obtenu.

- à couche ouverte : c'est le formalisme UHF (unrestricted Hartree-Fock) [122]. Les spin-orbitales sont construites à partir de deux ensembles {$\chi_{\lambda}^{\alpha}$} et {$\chi_{\lambda}^{\beta}$} correspondant respectivement aux deux états de spin $\alpha$ et $\beta$. On obtient alors deux équations de Roothaan couplées :

$$F^{\alpha} C^{\alpha} = \epsilon^{\alpha} S C^{\alpha} \quad (2.56)$$

$$F^{\beta} C^{\beta} = \epsilon^{\beta} S C^{\beta} \quad (2.57)$$

### 2.6.3 Les limites du modèle Hartree Fock

Le modèle Hartree-Fock conduit à la résolution rigoureuse du problème correspondant à un système à plusieurs électrons dans lequel les interactions $\frac{1}{r_{ij}}$ sont remplacées par les interactions électrostatiques moyennes données par les opérateurs de Coulomb $\hat{J}(i)$ et d'Echange $\hat{K}(i)$.

L'utilisation de la méthode HF avec une base de fonctions suffisamment importante tient compte d'environ 99% de l'énergie électronique. Néanmoins, les 1% restants sont souvent cruciaux pour décrire les phénomènes chimiques, comme par exemple la formation de liaisons. La différence entre l'énergie HF calculée et l'énergie exacte est appelée énergie de corrélation. Elle a pour expression :

$$E_{corr} = E_{exact} - E_{HF} \quad (2.58)$$

Pour approcher l'énergie exacte et la fonction d'onde exacte, il est donc nécessaire d'estimer la corrélation électronique qui se compose de deux contributions :
- La corrélation courte-portée, ou corrélation dynamique, est définie comme l'interaction à faible distance des électrons. Elle est prise en compte par les méthodes de type interaction de configurations (CI) [123]. Pour les systèmes où ce type de corrélation est important, les systèmes se trouvant dans leurs états fondamentaux près de la géométrie d'équilibre, la méthode Hartree-Fock va généralement surestimer les longueurs de liaisons.

- La corrélation statique est importante pour les molécules où l'état fondamental est décrit par plusieurs déterminants de Slater possédant de grands coefficients et pour des molécules excitées ou proches de la dissociation. Toutes les méthodes qui s'appuient sur un état obtenu après un calcul Hartree-Fock ne sont pas valables pour traiter la corrélation longue-portée. Par exemple, les forces de type van der Waals sont des forces de corrélation longue-portée, ayant un comportement asymptotique en $R^{-6}$.

Le développement des méthodes post Hartree-Fock a donc pour objectif la prise en compte d'une grande partie de l'énergie de corrélation.

## 2.7 Méthodes d'interaction de configurations (CI)

Afin de s'approcher de l'énergie exacte et de tenir compte de la corrélation électronique, il est indispensable de décrire la fonction d'onde électronique sous forme d'une combinaison linéaire de plusieurs configurations. La fonction d'onde HF, $\psi_0$, sert de point de départ ou de configuration de référence pour la description de la nouvelle fonction d'onde multiélectronique :

$$ \psi_{CI} = C_0\psi_0 + \sum_{I=1} C_I\psi_I \quad (2.59) $$

où $C_I$ sont les coefficients à déterminer par le principe variationnel et $\psi_I$ sont des déterminants construits à partir de l'excitation des spin-orbitales du déterminant $\psi_0$. En effet, la résolution des équations de Roothaan (éq. 2.51) pour un système à n électrons et M fonctions de base, conduit à n spin-orbitales moléculaires occupées et M-n spin-orbitales virtuelles. Les déterminants $\psi_I$ sont donc générés par l'occupation de ces spin-orbitales virtuelles suite à l'excitation des orbitales du déterminant Hartree-Fock.

En tenant compte de toutes les configurations qu'il est possible de construire à partir d'une base de spin-orbitales donnée (comprenant les spin-orbitales occupées et non occupées), le calcul de la corrélation électronique est total. Un tel calcul est appelé "Full Configuration Interaction". En utilisant une base infinie, la méthode d'interaction de configurations permet de résoudre l'équation de Schrödinger électronique d'une façon exacte. La génération d'un grand nombre de déterminants rend très difficile ce calcul.

**Principe variationnel :**

Les méthodes d'interaction de configuration reposent sur le principe variationnel. La fonction d'onde $\psi_{CI}$ est écrite sous la forme d'une combinaison linéaire de déterminants dont
## 2.7. Méthodes d'interaction de configurations (CI)

les coefficients $C_I$ sont ajustés de façon à obtenir une énergie minimale ($\frac{\partial E}{\partial C_I} = 0$).

$$ \psi_{CI} = C_0\psi_0 + C_s\psi_s + C_D\psi_D + C_T\psi_T + \dots = \sum_{I=0} C_I\psi_I \quad (2.60) $$

$\psi_0$ représente le déterminant de Slater qui sert de point de départ pour la description de la nouvelle fonction d'onde, $\psi_s$ le déterminant qui prend en compte toutes les simples excitations, $\psi_D$ le déterminant qui prend en compte toutes les doubles excitations et $\psi_T$ le déterminant qui prend en compte toutes les triples excitations.

L'énergie minimale est obtenue sous la contrainte de normalisation de la fonction d'onde $\Psi_{CI}$ via les multiplicateurs de Lagrange $\lambda$. La fonction de Lagrange $L$ s'exprime comme suit :

$$ L = \langle \psi_{CI} | \hat{H}_e | \psi_{CI} \rangle - \lambda [\langle \psi_{CI} | \psi_{CI} \rangle - 1] \quad (2.61) $$

avec

$$ \langle \psi_{CI} | \hat{H}_e | \psi_{CI} \rangle = \sum_{I=0} \sum_{J=0} C_I C_J \langle \psi_I | \hat{H}_e | \psi_J \rangle \quad (2.62) $$

$$ \langle \psi_{CI} | \psi_{CI} \rangle = \sum_{I=0} \sum_{J=0} C_I C_J \langle \psi_I | \psi_J \rangle = \sum_{I=0} C_I^2 \langle \psi_I | \psi_I \rangle = \sum_{I=0} C_I^2 \quad (2.63) $$

La dérivée de la fonction de Lagrange (équation 2.61) par rapport aux coefficients $C_i$ s'annule pour assurer la minimisation de l'énergie CI sous contrainte :

$$ \frac{\partial L}{\partial C_I} = 2 \sum_J C_J \langle \psi_I | \hat{H}_e | \psi_J \rangle - 2\lambda C_J = 0 \quad (2.64) $$

En introduisant la notation

$$ H_{IJ} = \langle \psi_I | \hat{H}_e | \psi_J \rangle $$

l'équation (2.64), pour chaque valeur de $I$, devient :

$$ \sum_J [H_{IJ} - \lambda \delta_{IJ}] C_J = 0 \quad (2.65) $$

Pour un ensemble de valeurs $I$, le problème variationnel est transformé en une résolution d'un ensemble d'équations séculaires que l'on peut écrire sous la forme matricielle suivante :

$$ \begin{pmatrix} H_{00} - \lambda & H_{01} & \cdots & H_{0J} & \cdots \\ H_{10} & H_{11} - \lambda & \cdots & H_{1J} & \cdots \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ H_{J0} & \cdots & \cdots & H_{JJ} - \lambda & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \end{pmatrix} \begin{pmatrix} C_0 \\ C_1 \\ \vdots \\ C_J \\ \vdots \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ \vdots \end{pmatrix} \quad (2.66) $$

Cette équation peut s'écrire aussi sous la forme compacte suivante :

$$ (\hat{H}_e - \lambda \hat{I}) C = 0 \quad (2.67) $$
La diagonalisation de la matrice CI (équation 2.66) permet d'obtenir un ensemble de valeurs propres correspondant aux énergies successives des différents états électroniques depuis l'état fondamental jusqu'aux états excités. Or ce calcul est extrêmement prohibitif. Pour le simplifier, on impose aux configurations $\psi_I$ de respecter des conditions de symétrie particulières [124]. Le principe est de construire ces configurations de telle sorte qu'elles possèdent une symétrie d'espace et de spin bien définie, autrement dit qu'elles puissent être représentées par le produit d'une fonction d'espace et de spin. On les appelle dans ce cas « Configurational State Functions » (CSFs). Si un déterminant de Slater n'a pas une symétrie de spin bien définie (pas fonction propre de $\hat{S}^2$), on construit $\psi_I$ comme une combinaison linéaire appropriée de déterminants. La même chose est réalisée pour la symétrie d'espace.

Plusieurs éléments de la matrice d'interaction de configuration peuvent s'annuler en prenant en compte plusieurs caractéristiques :

- L'opérateur hamiltonien électronique $\hat{H}_e$ (équation 2.9) ne couple pas les variables d'espace et de spin. Les éléments de matrice $\hat{H}_{IJ} = \langle \psi_I | \hat{H}_e | \psi_J \rangle$ associés à des configurations aux états de spin $\sigma_I$ et $\sigma_J$ orthogonaux sont nuls.

- L'hamiltonien électronique $\hat{H}_e$ (équation 2.9) appartient toujours à une représentation totalement symétrique. Si deux déterminants appartiennent à des représentations différentes, les éléments de matrice d'interaction de configuration sont alors nuls. La matrice (équation 2.66) peut ainsi être factorisée par blocs de symétries bien distincts, ce qui permet de séparer les calculs d'états électroniques de symétrie d'espace et de spin différents.

- A l'intérieur de chacun de ces blocs un certain nombre d'éléments de matrice sont nuls. En effet, si $\Phi_i$ et $\Phi_a$ sont les spin-orbitales (canoniques) issues des équations HF, elles sont fonctions propres de l'opérateur de Fock $\hat{F}$. Si on développe l'élément de matrice $\langle \psi_0 | \hat{H}_e | \psi_i^a \rangle$, on obtient :

$$ \langle \psi_0 | \hat{H}_e | \psi_i^a \rangle = \langle \Phi_i | \hat{h} | \Phi_a \rangle + \sum_j^n \left( \langle \Phi_i \Phi_j | \hat{g} | \Phi_a \Phi_j \rangle - \langle \Phi_i \Phi_j | \hat{g} | \Phi_j \Phi_a \rangle \right) \quad (2.68) $$

où $\psi_i^a$ est la configuration monoexcitée par rapport à $\psi_0$. En identifiant avec l'expression de l'opérateur de Fock (équation 2.44), on en déduit :

$$ \langle \psi_0 | \hat{H}_e | \psi_i^a \rangle = \langle \Phi_i | \hat{F} | \Phi_a \rangle = 0 \quad (2.69) $$

Les déterminants monoexcités $\psi_i^a$ n'interagissent donc pas avec la fonction de référence (théorème de Brillouin [125]). De plus, l'opérateur hamiltonien électronique $\hat{H}_e$ (équation 2.9) est constitué d'opérateurs mono et bi-électroniques. Or, on sait que de tels opérateurs ne peuvent coupler deux déterminants qui diffèrent par plus de deux spinorbitales (éléments de matrice nuls). La fonction de référence $\psi_0$ n'interagit donc qu'avec les configurations doublement excitées, $\psi_{ij}^{ab}$ :

$$ \langle \psi_0 | \hat{H}_e | \psi_{ij}^{ab} \rangle = \langle \Phi_i \Phi_j | \hat{g} | \Phi_a \Phi_b \rangle - \langle \Phi_i \Phi_j | \hat{g} | \Phi_b \Phi_a \rangle \quad (2.70) $$
La matrice CI ainsi construite est diagonale par blocs. Cependant, bien qu'elle possède de nombreux éléments de matrice nuls, le calcul Full CI implique un nombre infini de fonctions d'onde. Il est quasi impossible d'effectuer de tels calculs car le nombre de déterminants devient rapidement trop grand. Généralement, les méthodes utilisées sont dites CI tronquées, car on se limite à un nombre relativement faible d'excitations. Comme la prise en compte des mono-excitations n'apporte que peu d'amélioration par rapport au calcul Hartree-Fock, il faut au moins aller jusqu'à l'ordre deux des excitations afin d'observer les effets de la corrélation électronique. On parle alors de méthode CISD (Simples et Doubles excitations).

### 2.7.1 Méthode Multi Configuration Self-Consistent-Field (MCSCF)

La méthode Multi Configuration Self-Consistent-Field [126, 127, 128] est une méthode particulière d'interaction de configuration où les coefficients $C_I$ (éq. 2.59) de la combinaison des déterminants et les coefficients $C_{\mu i}$ (éq. 2.20) de la combinaison linéaire des orbitales atomiques $\chi_{\mu}$ sont optimisés. La minimisation de l'énergie par principe variationnel repose sur deux conditions :

$$ \frac{\delta E}{\delta C_I} = 0 \qquad (2.71) $$

$$ \frac{\delta E}{\delta C_{\mu i}} = 0 \qquad (2.72) $$

En partant d'une base de spin-orbitales {$\Phi_i$} bien adaptée au système à étudier, la Méthode Multi-Configuration Self-Consistent-Field permet de représenter la fonction d'onde avec un nombre limité de configuration CSF (Configurations State Functions). Cette combinaison linéaire forme ainsi la fonction d'onde de référence.

$$ \psi_{MCSCF} = \sum_{I} C_{I} \psi_{I} \qquad (2.73) $$

La forme des orbitales n'est plus fixée par un calcul préliminaire HF mais elle est adaptée à chaque étape du calcul MCSCF.

De façon générale, les méthodes MCSCF sont principalement utilisées pour générer une correction qualitative de la fonction d'onde électronique. Ceci permet la prise en compte de la partie statique de la corrélation électronique. L'augmentation du nombre des configurations permet alors d'inclure de plus en plus la corrélation dynamique jusqu'aux limites d'un calcul « full CI ». Dans ce cas, le traitement de la corrélation électronique est donc exact.

Le problème majeur de la méthode MCSCF est le choix du nombre limité de configurations nécessaires pour décrire les propriétés moléculaires du système considéré. Une autre méthode CASSCF est utilisée pour remédier à ce problème.
## 2.7.2 Méthode Complete Active Space SCF (CASSCF)

Une approche très utilisée pour surmonter les problèmes rencontrés dans la méthode MCSCF est la méthode « Complete Active Space Self-consistent Field » (CASSCF) [128, 129] où la sélection des configurations se fait selon une classification des orbitales moléculaires en sous-espaces :

- L'espace des orbitales inactives : les orbitales de plus basse énergie qui sont doublement occupées ou vides.

- L'espace des orbitales actives : les orbitales occupées de plus haute énergie et les orbitales inoccupées les plus basses en énergie obtenues à partir d'un calcul RHF [130]. Ces orbitales constituent l'espace actif.

- L'espace des orbitales moléculaires externes : les orbitales virtuelles qui ne sont jamais occupées et donc ne figurent dans aucune configuration.

La notation commune est [n,m] CASSCF où les n électrons sont distribués dans les m orbitales. Le développement full CI dans l'espace actif limite considérablement le nombre d'orbitales m et le nombre d'électrons n qui peuvent être traités par les méthodes CASSCF. Le choix des orbitales dites « importantes » devient alors crucial.

## 2.7.3 Méthode Multi Reference CI (MRCI)

La méthode MRCI (« Multi Reference CI ») [131] consiste à inclure dans la fonction d'onde toutes les configurations issues d'excitations électroniques simples et doubles à partir d'une fonction de référence, en général, une fonction issue d'un calcul de type MCSCF. Puisque la fonction de référence contient des configurations simplement et doublement excitées, la fonction d'onde MRCI est constituée de configurations triplement et quadruplement excitées dans l'espace interne (ensemble des orbitales inactives et actives) et simplement et doublement excitées dans l'espace externe (ensemble des orbitales virtuelles). On récupère ainsi une grande partie de la corrélation dynamique.

La fonction d'onde MRCI a pour expression :

$$ \psi_{MRCI} = \sum_{I} C_{I} \psi_{I} + \sum_{S} \sum_{a} C_{S}^{a} \psi_{S}^{a} + \sum_{D} \sum_{a,b} C_{D}^{ab} \psi_{D}^{ab} \quad (2.74) $$

a et b constituent les orbitales de l'espace externe, S et D sont les excitations simples et doubles des configurations de référence. $\psi_I$, $\psi_S^a$ et $\psi_D^{ab}$ représentent respectivement les configurations internes et les excitations simples et doubles d'une CSF de référence vers l'espace externe.

La méthode MRCI permet de décrire environ 95% de l'énergie de corrélation. Malheureusement, cette méthode n'est applicable que pour les petits systèmes moléculaires car elle est très coûteuse en temps de calcul et en capacité mémoire d'ordinateurs.
## 2.8 Méthode des Clusters Couplés (CC)

### 2.8.1 Méthode des Clusters Couplés CCSD(T)

La méthode des clusters couplés est l'une des plus utilisées à l'heure actuelle pour les calculs de précision, car elle permet de récupérer une fraction importante de la corrélation électronique [132, 133]. Cette méthode introduit dans la fonction d'onde de référence des corrections d'un certain type jusqu'à l'ordre infini.

La fonction d'onde s'écrit sous la forme :

$$ \psi^{cc} = e^{\hat{T}}\psi_0 \qquad (2.75) $$

avec

$$ e^{\hat{T}} = 1 + \hat{T} + \frac{1}{2}\hat{T}^2 + \frac{1}{6}\hat{T}^3 + \dots = \sum_{k \ge 0}^{\infty} \frac{\hat{T}^k}{k!} \qquad (2.76) $$

où $\hat{T}$ est l'opérateur cluster défini par :

$$ \hat{T} = \hat{T}_1 + \hat{T}_2 + \hat{T}_3 + \dots + \hat{T}_N \qquad (2.77) $$

L'opérateur $\hat{T}_i$ agit sur $\psi_0$ pour donner des excitations simples, doubles, triples, etc. On peut les définir comme suit :

$$ \hat{T}_1 \psi_0 = \sum_{i}^{\text{occ}} \sum_{a}^{\text{vir}} t_i^a \psi_i^a \qquad (2.78) $$

$$ \hat{T}_2 \psi_0 = \sum_{i<j}^{\text{occ}} \sum_{a<b}^{\text{vir}} t_{ij}^{ab} \psi_{ij}^{ab} \qquad (2.79) $$

Les coefficients t des développements sont nommés amplitudes.

A partir des équations (2.76) et (2.77), $e^{\hat{T}}$ peut s'écrire sous la forme suivante :

$$ e^{\hat{T}} = 1 + \hat{T}_1 + (\hat{T}_2 + \frac{1}{2}\hat{T}_1^2) + (\hat{T}_3 + \hat{T}_2\hat{T}_1 + \frac{1}{6}\hat{T}_1^3) + (\hat{T}_4 + \hat{T}_3\hat{T}_1 + \frac{1}{2}\hat{T}_2^2) + \frac{1}{2}\hat{T}_2\hat{T}_1^2) + \frac{1}{24}\hat{T}_1^4) + \dots \qquad (2.80) $$

Le premier terme génère la fonction de référence, le second tous les états simplement excités, la première parenthèse tous les états doublement excités, la seconde parenthèse les états triplement excités, etc.

L'équation de Schrödinger électronique s'écrit sous la forme :

$$ \hat{H}_e e^{\hat{T}} \psi_0 = E_{cc} e^{\hat{T}} \psi_0 \qquad (2.81) $$

Multiplions par $\psi_0^*$ à gauche et intégrons :

$$ \langle \psi_0 | \hat{H}_e e^{\hat{T}} | \psi_0 \rangle = E_{cc} \langle \psi_0 | e^{\hat{T}} \psi_0 \rangle = E_{cc} \langle \psi_0 | (1 + \hat{T}_1 + \hat{T}_2 + \dots) \psi_0 \rangle \qquad (2.82) $$

$$ \langle \psi_0 | \hat{H}_e e^{\hat{T}} | \psi_0 \rangle = E_{cc} \qquad (2.83) $$
Puisque l'Hamiltonien électronique (équation 2.28) ne contient que des opérateurs mono et biélectroniques, l'équation (2.83) peut être réécrite comme suit :

$$
\begin{align}
\langle \psi_0 | \hat{H}_e (1 + \hat{T}_1 + \hat{T}_2 + \frac{1}{2} \hat{T}_1^2) | \psi_0 \rangle &= E_{cc} \nonumber \\
\langle \psi_0 | \hat{H}_e | \psi_0 \rangle + \langle \psi_0 | \hat{H} | \hat{T}_1 \psi_0 \rangle + \langle \psi_0 | \hat{H}_e | \hat{T}_2 \psi_0 \rangle + \frac{1}{2} \langle \psi_0 | \hat{H}_e | \hat{T}_1^2 \psi_0 \rangle &= E_{cc} \nonumber \\
E_0 + \sum_{i}^{occ} \sum_{a}^{virt} t_i^a \langle \psi_0 | \hat{H}_e | \psi_i^a \rangle + \sum_{i<j}^{occ} \sum_{a<b}^{virt} (t_{ij}^{ab} + t_i^a t_j^b - t_i^b t_j^a) \langle \psi_0 | \hat{H}_e | \psi_{ij}^{ab} \rangle &= E_{cc} \tag{2.84}
\end{align}
$$

Lorsque les orbitales HF sont utilisées pour la construction du déterminant de Slater, les éléments de matrice du premier terme $\langle \psi_0 | \hat{H}_e | \psi_i^a \rangle$ s'annulent selon le théorème de Brillouin [125]. Le deuxième terme est une somme sur les intégrales biélectroniques.

$$
E_0 + \sum_{i<j}^{occ} \sum_{a<b}^{vir} (t_{ij}^{ab} + t_i^a t_j^b - t_i^b t_j^a) (\langle \Phi_i \Phi_j | \Phi_a \Phi_b \rangle - \langle \Phi_i \Phi_j | \Phi_b \Phi_a \rangle) = E_{cc} \quad (2.85)
$$

L'énergie Coupled Cluster est déterminée par les amplitudes d'ordre un $t_i^a$ et d'ordre deux $t_{ij}^{ab}$ et par des intégrales sur les orbitales moléculaires à deux électrons. Cependant, la prise en compte de toutes les excitations n'est possible que pour les systèmes de petite taille. Pour traiter des systèmes plus conséquents, il faut «tronquer» les équations CC à un certain niveau d'excitation. Ceci mène à une approximation de $E_{CC}$.

Les méthodes des clusters couplés tronqués sont classées selon le nombre d'excitations permises au niveau de l'opérateur $\hat{T}$ comme suit :

- CCS pour les excitations simples

- CCD pour les excitations doubles

- CCT pour les excitations triples

- CCQ pour les excitations quadruples

La «troncation» la plus répandue consiste à s'arrêter aux doubles excitations ($\hat{T} = \hat{T}_1 + \hat{T}_2$) qui correspondent aux méthodes CCSD (Coupled Cluster restricted to Single and Double excitations). La troncation d'ordre supérieur avec $\hat{T} = \hat{T}_1 + \hat{T}_2 + \hat{T}_3$ se nomme CCSDT. On complète parfois le calcul CCSD par l'inclusion des triples excitations de façon perturbative. On parle alors de la méthode CCSD(T) qui est intermédiaire entre CCSD et CCSDT et qui a l'avantage d'être peu coûteuse en temps de calcul et applicable à des systèmes à plusieurs électrons en utilisant des bases relativement grandes.

La seule limitation de la méthode des Clusters Couplés provient de la fonction d'onde de référence qui est monoconfigurationnelle et qui ne permet donc d'appliquer ces méthodes qu'aux états fondamentaux de molécules, ou bien à l'état électronique le plus bas de chaque symétrie et multiplicité de spin.
## 2.8.2 Méthode des Clusters Couplés explicitement corrélés : CCSD(T)-F12

La convergence de l'énergie de corrélation évolue très lentement en fonction de la taille de la base. La raison de cette défaillance est liée à la difficulté de décrire la singularité qui apparaît lorsque la distance interélectronique, $r_{12}$, tend vers zéro. A une coalescence d'électrons, la fonction d'onde se comporte comme [134] :

$$ \lim_{r_{12} \to 0} \left( \frac{\partial \psi}{\partial r_{12}} \right) = \frac{1}{2} \psi(r_{12} = 0) \qquad (2.86) $$

Afin de résoudre ce problème et d'accélérer la convergence vers la limite des bases, les méthodes explicitement corrélées sont utilisées. Elles prennent en compte explicitement les distances électron-électron. La fonction d'onde utilisée dans l'approche CCSD-F12 avec excitations simples et doubles a pour expression [135] :

$$ \psi_{CCSD} = e^{\hat{T}_1 + \hat{T}_2} \psi_0 \qquad (2.87) $$

où $\psi_0$ est la fonction d'onde de référence HF et les opérateurs d'excitation simples et doubles sont [135] :

$$ \hat{T}_1 = t_a^i \hat{E}_{ai} \qquad (2.88) $$

$$ \hat{T}_2 = T_{ab}^{ij} \hat{E}_{ai} \hat{E}_{bj} + \mathcal{T}_{\alpha\beta}^{ij} \hat{E}_{\alpha i} \hat{E}_{\beta j} \qquad (2.89) $$

Les indices $i$ et $j$ correspondent aux orbitales occupées. $a$ et $b$ correspondent aux orbitales virtuelles $\alpha$ et $\beta$ à l'orbitale de la base complète. $\hat{E}_{ai}$ et $\hat{E}_{\alpha i}$ $\hat{E}_{bj}$ sont les opérateurs d'excitation à un et deux électrons, respectivement. $t_a^i$ et $T_{ab}^{ij}$ sont les amplitudes simples et doubles. La théorie F12 fait introduire des amplitudes additionnelles, $\mathcal{T}_{\alpha\beta}^{ij}$, approximées par [135, 136] :

$$ \mathcal{T}_{\alpha\beta}^{ij} = \mathcal{T}_{mn}^{ij} \mathcal{F}_{\alpha\beta}^{mn} \qquad (2.90) $$

$$ \mathcal{F}_{\alpha\beta}^{mn} = \langle mn | \hat{Q}_{12} F_{12} | \alpha\beta \rangle \qquad (2.91) $$

$$ \hat{Q}_{12} = 1 - |rs\rangle\langle rs| - |mx\rangle\langle mx| - |xm\rangle\langle jm| \qquad (2.92) $$

$F_{12} \equiv F(r_{12})$ est un facteur de corrélation adapté à courte portée et les $\mathcal{T}_{kl}^{ij}$ sont les amplitudes utilisées dans le traitement F12. Les indices $r$ et $s$ désignent la base complète des orbitales moléculaires. $x$ est la base auxiliaire complémentaire.

Le projecteur $\hat{Q}_{12}$ assure de façon explicite l'orthogonalité des configurations corrélées à la fonction de référence Hartree-Fock ainsi que les doubles excitations [135] :

$$ |\psi_{ij}^{mn}\rangle = \mathcal{F}_{\alpha\beta}^{mn} \hat{E}_{\alpha i} \hat{E}_{\beta j} |\psi\rangle \qquad (2.93) $$

Les termes explicitement corrélés ont amélioré la fonction d'onde décrivant des électrons à faibles distances interatomiques. Le produit des orbitales spatiales $\phi_i(\vec{r}_1) \phi_j(\vec{r}_2)$ est remplacé par une fonction de corrélation à courte portée telle que :

$$ |u_{ij}(\vec{r}_1, \vec{r}_2)\rangle = \mathcal{T}_{mn}^{ij} \hat{Q}_{12} \hat{F}_{12}(r_{12}) |\phi_m(\vec{r}_1) \phi_n(\vec{r}_2)\rangle \qquad (2.94) $$
Selon Knizia, Adler et Werner [136], les termes pour lesquels $mn = ij$ et $mn = ji$ sont les plus importants. Si dans l'expression de la fonction d'onde de référence le terme $|u_{ij}\rangle$ est introduit, la probabilité de trouver deux électrons (i,j) proches l'un de l'autre s'annule. Le facteur de corrélation à courte portée $\hat{Q}_{12}\hat{F}_{12}|\phi_m\phi_n\rangle$ permet de décrire correctement la fonction d'onde pour une distance interatomique qui tend vers zéro [136]. Le facteur est pris comme une simple fonction de Slater qui peut être approximée par une combinaison linéaire de gaussiennes [135] :

$$F(r_{12}) = -\frac{1}{\gamma} \exp(-\gamma r_{12}) \approx \sum_{i} c_i \exp(-\alpha_i r_{12}^2) \quad (2.95)$$

Pour la méthode CCSD(T)-F12a, le calcul de l'énergie de corrélation est approximé selon Adler, Knizia et Werner [135] et elle correspond à la somme de l'énergie de corrélation CCSD(T) et la correction MP2-F12 (E(MP2-F12)-E(MP2)).

## 2.9 Théorie de la perturbation de Møller-Plesset

Une manière plus approchée de traiter le problème de la corrélation électronique est d'utiliser la théorie des perturbations qui consiste à appliquer une perturbation $\hat{W}$ à un hamiltonien de référence $\hat{H}_0$ dont on connaît les valeurs propres et les vecteurs propres. L'hamiltonien du système s'écrit alors :

$$\hat{H}_e^p = \hat{H}_0 + \lambda \hat{W} \quad (2.96)$$

$$\hat{H}_0 \psi_i = E_i \psi_i, \quad i = 1, 2, \dots, \infty \quad (2.97)$$

où l'indice $p$ indique que l'hamiltonien est perturbé et $\lambda$ est un paramètre déterminant l'amplitude de la perturbation. Les solutions de l'opérateur hamiltonien non perturbé forment une base complète orthonormée.

Les méthodes de perturbation [137] peuvent être utilisées pour corriger les approches qui utilisent l'approximation de particules indépendantes. L'équation de Schrödinger électronique perturbée est donnée par l'équation suivante :

$$\hat{H}_e^p \psi = W \psi \quad (2.98)$$

Dans le cas où $\lambda = 0$, $\hat{H}_e^p = \hat{H}_0$, $\psi = \psi_0$ et $W = E_0$. Lorsque la perturbation augmente, l'énergie et la fonction d'onde sont développées en puissance de $\lambda$ :

$$W = \lambda^0 W_0 + \lambda^1 W_1 + \lambda^2 W_2 + \lambda^3 W_3 + \dots \quad (2.99)$$

$$\psi = \lambda^0 \psi_0 + \lambda^1 \psi_1 + \lambda^2 \psi_2 + \lambda^3 \psi_3 + \dots \quad (2.100)$$
Tous les termes de corrections sont orthogonaux avec la fonction d'onde de référence :

$$
\begin{align*}
\langle \psi | \psi_0 \rangle &= 1 \\
\langle \psi_0 + \lambda^1 \psi_1 + \lambda^2 \psi_2 + \lambda^3 \psi_3 \dots | \psi_0 \rangle &= 1 \\
\langle \psi_0 | \psi_0 \rangle + \lambda^1 \langle \psi_1 | \psi_0 \rangle + \lambda^2 \langle \psi_2 | \psi_0 \rangle + \lambda^3 \langle \psi_3 | \psi_0 \rangle + \dots &= 1 \\
\langle \psi_{i \neq 0} | \psi_0 \rangle &= 0 \tag{2.101}
\end{align*}
$$

L'équation de Schrödinger électronique (2.98) devient :

$$
\begin{align}
(\hat{H}_0 + \lambda \hat{W})(\lambda^0 \psi_0 + \lambda^1 \psi_1 + \lambda^2 \psi_2 + \lambda^3 \psi_3 + \dots) &= \nonumber \\
(\lambda^0 W_0 + \lambda^1 W_1 + \lambda^2 W_2 + \lambda^3 W_3 + \dots)(\lambda^0 \psi_0 + \lambda^1 \psi_1 + \lambda^2 \psi_2 + \lambda^3 \psi_3 + \dots) &= \tag{2.102}
\end{align}
$$

Dans le cadre de la théorie de Moller Plesset (MP), l'opérateur non perturbé correspond
à la somme des opérateurs de Fock $\hat{F}(i)$ :

$$
\begin{align}
\hat{H}_0 &= \sum_{i=1}^{n} \hat{F}(i) = \sum_{i=1}^{n} [\hat{h}(i) + \sum_{j=1}^{n} (\hat{J}(j) - \hat{K}(j))] \nonumber \\
&= \sum_{i=1}^{n} \hat{h}(i) + \sum_{i=1}^{n} \sum_{j=1}^{n} \langle \hat{g}(i, j) \rangle = \sum_{i=1}^{n} \hat{h}(i) + 2 \langle \hat{V}_{ee} \rangle \tag{2.103}
\end{align}
$$

Néanmoins, la répulsion électronique moyenne $\langle \hat{V}_{ee} \rangle$ est comptée doublement par la somme sur les opérateurs de Fock. Ainsi, l'opérateur de perturbation, appelé potentiel de fluc- tuation, s'écrit :

$$
\hat{W} = \hat{H}_e^p - \hat{H}_0 = \sum_{i=1}^{n} \sum_{j>i}^{n} \hat{g}(i, j) - \sum_{i=1}^{n} \sum_{j=1}^{n} \langle \hat{g}(i, j) \rangle = \hat{V}_{ee} - 2 \langle \hat{V}_{ee} \rangle \quad (2.104)
$$

Compte tenu des expressions des équations (2.99) et (2.100), nous pouvons attribuer à la
méthode MP une certaine hiérarchie : MP0, MP1, MP2, etc, permettant d'introduire de
mieux en mieux la corrélation électronique :

-MP0 : la fonction d'onde d'ordre zéro n'est autre que le déterminant de HF et l'énergie d'ordre zéro est la somme des énergies mono-électroniques $\epsilon_i$ (équation 2.55) :

$$
W_0 = \langle \psi_0 | \hat{H}_0 | \psi_0 \rangle = \langle \psi_0 | \sum_{i=1}^{n} \hat{F}(i) | \psi_0 \rangle = \sum_{i=1}^{n} \epsilon_i \quad (2.105)
$$

A ce niveau de perturbation, la répulsion inter-électronique est ignorée.

-MP1 : l'énergie calculée est obtenue en ajoutant au terme d'énergie MP0 les intégrales de Coulomb et d'Exchange :

$$
W_1 = \langle \psi_0 | \hat{W} | \psi_0 \rangle = \langle \hat{V}_{ee} \rangle - 2 \langle \hat{V}_{ee} \rangle = - \langle \hat{V}_{ee} \rangle \quad (2.106)
$$

L'énergie du premier ordre correspond exactement à l'énergie HF.
–MP2 : L'énergie de corrélation électronique n'est prise en compte qu'à partir des termes de correction de second ordre. L'énergie calculée au niveau MP2 est la somme de l'énergie MP1 (E(HF)) plus un terme de correction E₂ qui représente l'énergie issue de la répulsion électronique. Les expressions des corrections de second ordre de l'énergie incluent des éléments de matrice de perturbation entre l'état de référence HF et les états excités [43]. Puisque la perturbation est décrite par des opérateurs à deux électrons, tous les éléments de matrice relatifs à des états excités triple, quadruple, etc, sont nuls. L'énergie de correction de second ordre qui représente la première contribution dans l'expression de l'énergie de corrélation inclut la somme des déterminants doublement excités :

$$W_2 = \sum_{j>i}^{occ} \sum_{b>a}^{vir} \frac{\langle \psi_0 | \hat{W} | \psi_{ij}^{ab} \rangle \langle \psi_{ij}^{ab} | \hat{W} | \psi_0 \rangle}{E_0 - E_{ij}^{ab}} \quad (2.107)$$

où i et j dénomment deux électrons de l'orbitale occupée, a et b sont deux électrons de l'orbitale virtuelle. La formule explicite de la correction MP2 est donnée par :

$$E_2 = \sum_{j>i}^{occ} \sum_{b>a}^{vir} \frac{\langle \langle \Phi_i \Phi_j | \Phi_a \Phi_b \rangle \rangle \langle \Phi_i \Phi_j | \Phi_b \Phi_a \rangle}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b} \quad (2.108)$$

L'avantage principal des méthodes perturbatives est la cohérence en taille. Cette propriété signifie que l'énergie varie linéairement avec la taille du système. Ceci est nécessaire pour un calcul correct des énergies d'interaction. Cependant, la théorie des perturbations présente certains défauts. L'inconvénient majeur provient du choix de considérer la fonction d'onde d'ordre zéro comme approximation raisonnable pour la fonction d'onde réelle. En effet, les termes de correction seront plus grands si la fonction d'onde HF est de qualité médiocre et le développement perturbatif ne converge pas systématiquement. De plus, ce n'est pas une méthode variationnelle.

## 2.10 La théorie de la fonctionnelle de la densité

La théorie de la fonctionnelle de la densité (DFT) [138, 139] est une alternative aux théories citées précédemment qui se fondent sur une fonction d'onde multiélectronique. Pour calculer l'énergie et déterminer les propriétés de l'état fondamental d'un système à n électrons, il n'est plus nécessaire de connaître la fonction d'onde $\Psi(r_1, r_2...r_n)$ mais plutôt la densité électronique $\rho(r)$. La fonction d'onde d'un système à n électrons dépend des 3n coordonnées d'espace (n nombre total de particules du système considéré) mais aussi de n coordonnées de spin, entraînant une dépendance à 4n variables. La densité électronique dépend uniquement des 3 coordonnées d'espace et de la coordonnée de spin du système considéré. Ainsi, alors que la complexité de la fonction d'onde augmente avec le nombre de variables, la densité électronique, elle, reste indépendante de la taille du système. Contrairement aux méthodes ab initio qui utilisent un hamiltonien exact et font des approximations sur la fonctions d'onde, les méthodes DFT font des approximations sur l'hamiltonien et calculent exactement la densité électronique associée. Les premières bases 64
de la DFT ont été données en 1964 par Hohenberg et Kohn, [138] qui ont démontré la possibilité de calculer les propriétés d'un système à N électrons en interaction dans un potentiel externe à l'aide de la densité électronique. Ensuite, Kohn et Sham [139] ont rendu possible les applications de la DFT en chimie quantique.

## 2.10.1 Densité électronique

La densité électronique $\rho(\vec{r})$ est définie comme étant la distribution moyenne dans l'espace d'un électron dans le champ des autres électrons dans un élément de volume centré sur la position $\vec{r}$ tel que :

$$ \rho(\vec{r}) = n \int |\Psi(\vec{r}_1, \vec{r}_2, \dots, \vec{r}_n)|^2 d\vec{r}_2 \dots d\vec{r}_n \quad (2.109) $$

$\vec{r}$ inclut implicitement les coordonnées de spin ; l'intégrale sous-entend la sommation sur les états de spin [139].

La densité électronique possède deux propriétés importantes :

$$ \lim_{\vec{r} \to \infty} \rho(\vec{r}) = 0 \\ \int \rho(\vec{r}) d\vec{r} = n $$

## 2.10.2 Théorème de Hohenberg et Kohn

En mécanique quantique, si l'on connaît le nombre d'électrons d'un système de particules et le potentiel externe $E_{Ne}$ (toute interaction entre les électrons et les noyaux, est connue) auquel il est soumis, la fonction d'onde de l'état fondamental et la densité électronique associée sont déterminées. La réciproque fut démontrée en 1964. C'est le premier théorème de Hohenberg et Kohn : [138]

« Pour tout système de particules en interaction dans un potentiel externe, ce potentiel est uniquement déterminé, à une constante additive près, par la densité $\rho(\vec{r})$ du système dans son état fondamental ».

Ce premier théorème énonce que toute observable d'un état fondamental stationnaire non-dégénéré peut être calculée, exactement en théorie, par le biais de la densité électronique de cet état fondamental.

Le second théorème de Hohenberg et Kohn [138] traite le principe variationnel pour les méthodes DFT : « L'énergie de l'état fondamental est le minimum de la fonctionnelle $E[\rho]$ »

tel que :

$$ E[\rho] = \min_{\rho(\vec{r})} E[\rho(\vec{r})] \quad (2.110) $$
Ces théorèmes ne sont valables que pour l'état fondamental non dégénéré. La fonctionnelle énergie est décomposée en termes d'énergie cinétique T et de potentiels d'interaction noyaux-électrons $E_{Ne}$ et de répulsion noyau-noyau $E_{ee}$, toutes variables de $\rho$ :

$$E[\rho] = T[\rho] + E_{ee}[\rho] + E_{Ne}[\rho] \qquad (2.111)$$

Les deux premiers termes sont rassemblés dans une fonctionnelle universelle $F[\rho]$ indépendante du potentiel externe.

Le terme décrivant l'attraction entre les électrons et les noyaux s'écrit :

$$E_{Ne}[\rho] = - \sum_{\alpha}^{N} \int \frac{Z_{\alpha} \rho(\vec{r})}{|\vec{R}_{\alpha} - \vec{r}|} d\vec{r} \qquad (2.112)$$

L'énergie $E_{ee}[\rho]$ est divisée en une partie de Coulomb $J[\rho]$ et une partie d'Echange $K[\rho]$ par analogie à la théorie de HF (équation 2.36). La fonctionnelle $J[\rho]$ est donnée par l'expression classique :

$$J[\rho] = \frac{1}{2} \iint \frac{\rho(\vec{r})\rho(\vec{r}')}{|\vec{r} - \vec{r}'|} d\vec{r}d\vec{r}' \qquad (2.113)$$

Pour déterminer les énergies cinétique $T[\rho]$ et d'échange $K[\rho]$, une première approximation est de considérer le nuage électronique comme un gaz uniforme d'électrons indépendants : [130]

$$T_{TF}[\rho] = C_F \int \rho^{5/3}(\vec{r}) d\vec{r} \qquad (2.114)$$

$$K_D[\rho] = -C_x \int \rho^{4/3}(\vec{r}) d\vec{r} \qquad (2.115)$$

où

$$C_F = \left(\frac{3}{10}\right) \times (3\pi^2)^{2/3} \qquad (2.116)$$

$$C_x = \left(\frac{3}{4}\right) \left(\frac{3}{\pi}\right) \qquad (2.117)$$

La fonctionnelle d'énergie nommée $E_{TF}[\rho]$ a pour expression :

$$E_{TF}[\rho] = T_{TF}[\rho] + E_{Ne}[\rho] + J[\rho] \qquad (2.118)$$

Ce modèle est connu sous le nom du théorème de Thomas-Fermi (TF). L'ajout du terme $K_D[\rho]$ constitue le modèle de Thomas-Fermi-Dirac.

Le théorème de Hohenberg et Kohn a été appliqué en chimie quantique par Kohn et Sham (1965). L'objectif était de trouver une fonctionnelle d'échange-corrélation permettant de connecter la densité avec l'énergie et de traiter la corrélation électronique.
## 2.10.3 Théorème de Kohn et Sham

Kohn et Sham (1965) [139] proposèrent d'introduire la notion d'orbitale dans le formalisme DFT dont l'idée est de diviser la fonctionnelle de l'énergie cinétique en deux parties :

- un terme majoritaire qui peut être calculé exactement représente l'énergie cinétique calculée selon le modèle de particules indépendantes. Le modèle KS est étroitement lié à la méthode HF. Les mêmes formules pour l'énergie cinétique, le potentiel électrons-noyaux et la répulsion de Coulomb sont utilisés.

- l'autre partie est un terme de correction.

En se basant sur des méthodes de fonction d'onde, Kohn et Sham proposèrent d'introduire la notion d'orbitale dans le formalisme DFT. Dans ce nouveau cadre, l'énergie cinétique dépend explicitement des orbitales et s'écrit :

$$ \hat{T}[\rho] = \sum_{i=1}^{n} \langle \psi_e | - \frac{1}{2} \nabla_i^2 | \psi_e \rangle \quad (2.119) $$

$\psi_e$ est la fonction d'onde électronique de l'état fondamental.

Kohn et Sham [139] adaptèrent ensuite un modèle Hartree-Fock à particules indépendantes à la méthode DFT.

Supposons qu'il existe une fonction d'onde à un seul déterminant de Slater $\psi_s$ reproduisant la densité électronique de la fonction d'onde électronique exacte dans le cas où les électrons associés aux orbitales $\Phi_i$ n'interagissent pas entre eux. Dans ce nouveau cadre, l'expression de la fonctionnelle d'énergie cinétique des électrons sans interaction est donnée par :

$$ \hat{T}_S[\rho] = \sum_{i=1}^{n} \langle \Phi_i | - \frac{1}{2} \nabla_i^2 | \Phi_i \rangle \quad (2.120) $$

S indique que l'énergie cinétique est calculée à partir d'un seul déterminant de Slater.

Cependant, l'interaction électron-électron est introduite par le terme d'échange-corrélation, $E_{xc}$. La fonctionnelle de l'énergie de Kohn et Sham s'écrit alors comme suit :

$$ E_{DFT}[\rho] = T_S[\rho] + J[\rho] + E_{Ne}[\rho] + E_{xc}[\rho] \quad (2.121) $$

En égalisant $E_{DFT}[\rho]$ à l'énergie exacte $E[\rho]$ (équation 2.111), $E_{xc}[\rho]$ s'exprime comme :

$$ E_{xc}[\rho] = (T[\rho] - T_S[\rho]) + (E_{ee}[\rho] - J[\rho]) \quad (2.122) $$

La première parenthèse correspond à la différence entre les énergies cinétiques des électrons indépendants et des électrons en interaction. Cette contribution n'existe pas dans le terme de corrélation des méthodes de fonctions d'onde. $E_{ee}$ est l'énergie totale liée aux
répulsions électroniques et $J[\rho]$ l'énergie coulombienne classique de répulsion électron-électron.

Les composantes qui constituent l'essentiel de la fonctionnelle énergie totale, c'est-à-dire $T_S[\rho]$, $J[\rho]$ et $E_{Ne}[\rho]$ peuvent être calculées exactement à partir de la densité électronique qui est obtenue facilement dans le modèle orbitalaire :

$$ \rho(\vec{r}) = \sum_{i}^{n} |\Phi_i|^2 \qquad (2.123) $$

L'application du principe variationnel aux fonctions $\phi_i$ (partie spatiale des spin-orbitales) conduit aux équations de Kohn Sham qui sont similaires aux équations de Hartree-Fock :

$$ \hat{h}_{KS} \phi_i = \epsilon_i \phi_i \qquad (2.124) $$

$$ \left[-\frac{1}{2}\nabla_i^2 + \hat{V}_{eff}\right]\phi_i = \epsilon_i \phi_i \qquad (2.125) $$

$$ \left[-\frac{1}{2}\nabla_i^2 + \hat{V}_{Ne}(r) + \int \frac{\rho(r)}{r-r'} dr + \hat{V}_{xc}(r)\right]\phi_i = \epsilon_i \phi_i \qquad (2.126) $$

où $\hat{h}_{KS}$ est un opérateur monoélectronique dont les $\phi_i$ sont les vecteurs propres. $\hat{V}_{eff}$ est appelé opérateur de potentiel effectif de Kohn-Sham.

Une approche SCF est appliquée afin de résoudre les équations KS (2.126) à savoir :

- Choix d'un jeu d'essai d'orbitales,

- Calcul de la densité d'essai,

- Construction du potentiel effectif de Kohn-Sham,

- Diagonalisation de la matrice $\langle \phi_i | h_{KS} | \phi_j \rangle$.

Les orbitales obtenues constituent le jeu d'essai pour l'itération suivante. Cette procédure est répétée jusqu'à convergence de l'énergie.

Malgré le fait que la DFT soit une théorie exacte, la fonctionnelle permettant d'obtenir $E_{xc}$ reste inconnue. Il est donc nécessaire de l'approximer. Plusieurs types d'approximations existent :

- L'approche LDA (Local Density Approximation) qui consiste à considérer la densité comme celle d'un gaz d'électrons homogène, c'est-à-dire que la densité varie lentement en fonction de la position.

- L'approche LSDA (Local Spin Density Approximation) qui est un modèle LDA prenant en compte le spin $(\alpha, \beta)$,

- Les méthodes GGA (Generalised Gradient Approximation) où l'énergie d'échange-corrélation dépend non seulement de la densité électronique mais aussi de ses dérivées,

- Les méthodes méta-GGA où la fonctionnelle dépend de la densité électronique, de son gradient et du Laplacien de la densité électronique,

- Les fonctionnelles hybrides où la fonctionnelle d'échange contient différents pourcentages du terme d'échange obtenu par les méthodes HF.
Ces différentes approches sont détaillées dans l'annexe A.

## 2.11 Bases atomiques

Les orbitales moléculaires selon la méthode LCAO sont définies comme étant une combinaison linéaire d'orbitales atomiques modélisées par des fonctions de base. Dans les premiers calculs de chimie quantique, les orbitales de type Slater (STO) ont été utilisées pour décrire la distribution spatiale des électrons :

$$f_{STO} = Nr^n x^{l_x} x^{l_y} x^{l_z} \exp(-\xi r) \qquad (2.127)$$

N est le facteur de normalisation, n est le nombre quantique principal de la couche électronique, $\xi$ est une constante liée à la charge effective du noyau, souvent estimée à l'aide des règles empiriques proposées par Slater, $l_x$, $l_y$ et $l_z$ sont des entiers naturels dont la somme définit le type d'orbitale utilisée.

Pour surmonter les difficultés numériques des fonctions STO, des fonctions de type gaussienne (GTO) sont utilisées, de forme générale :

$$f_{GTO} = N x^{l_x} x^{l_y} x^{l_z} \exp(-\alpha r^2) \qquad (2.128)$$

$\alpha$ est un paramètre représentant l'extension radiale. La dépendance en $r^2$ de la partie radiale des fonctions GTOs, rend ces dernières moins appropriées, par rapport aux fonctions STO, pour représenter le comportement des électrons lorsqu'ils sont très proches du noyau et lorsqu'ils y sont très éloignés. Cependant, les orbitales gaussiennes sont mieux adaptées aux calculs numériques.

Plusieurs gaussiennes sont nécessaires pour obtenir le comportement d'une orbitale de Slater ; une fonction atomique s'exprime alors sous la forme d'une combinaison linéaire de fonctions gaussiennes élémentaires :

$$X(CGTO) = \sum_{i}^{k} c_g f_{PGTO} \qquad (2.129)$$

Les coefficients $c_g$ sont des coefficients de contraction fixés par des calculs préalables sur l'état fondamental des atomes isolés. $k$ orbitales primitives (Primitive Gaussian Type Orbitals) sont contractées en une orbitale CGTO (Contracted Gaussian Type Orbitals). Cette contraction est utile pour les orbitales décrivant les électrons de coeur alors que pour les orbitales de valence qui sont fortement déformées lors de la formation des liaisons chimiques, moins de fonctions contractées sont utilisées afin de garder une certaine flexibilité. Dans ce cas, chaque orbitale atomique est représentée par plusieurs fonctions contractées. On appelle ce type de base *Valence Double*, *Triple*, *Quadruple* et *Quintuple* zeta (VDZ, VTZ,...) lorsqu'à chaque OA sont associées deux, trois ...fonctions contractées.
2.11.1 Bases split valence

Les bases « split valence » consistent à décrire les orbitales de coeur à l'aide d'une com-
binaison de fonctions gaussiennes et les orbitales de valence à partir de bases multi-zeta.
Parmi les bases les plus connues introduites par Pople et al. [140], on peut citer des bases
comme 6-31G, 6-311G...,etc. La nomenclature pour ces bases indique la façon dont elles
sont construites. Par exemple, pour un atome de carbone dans la base 6-31G, le nombre
6 indique que l'orbitale 1s est représenté par 6 fonctions gaussiennes. les orbitales de
valence 2s et 2p sont décrites par 3 gaussiennes et une fonction gaussienne, respectivement.

2.11.2 Fonctions de polarisation

Les bases « split valence » ont permis de changer la taille des orbitales mais pas
leurs formes. Lors d'une réaction chimique, les atomes se rapprochent de plus en plus et
les OAs de valence se déforment. Afin d'augmenter la flexibilité de la base utilisée, des
orbitales d'ordre supérieur sont ajoutées aux orbitales considérés. Par exemple, pour la
base 6-31G(d,p), une fonction de type *d* est ajoutée à chaque atome lourd (comme le
carbone et l'oxygène) et une orbitale de type *p* est ajoutée à chaque atome d'hydrogène.
Les bases polarisées développées par Dunning [141] trouvent un grand succès puisqu'elles
permettent une convergence de l'énergie électronique vers la limite de base complète.
Le nom de ces bases commence avec cc-p, pour correlation consistent polarized. Pour
les orbitales de valence, ces bases sont *Double*, *Triple*, *Quadruple* et *Quintuple* zeta.
Ces bases sont alors notées cc-pVDZ, cc-pVTZ...etc. Elles incluent successivement des
fonctions de polarisation de plus en plus grandes (d, f, g...).

2.11.3 Fonctions diffuses

Pour les molécules comportant des liaisons ioniques ou contenant des paires libres
d'électrons, il est indispensable de tenir compte de la densité électronique significative
loin du centre des noyaux. Pour cela, des fonctions diffuses sont ajoutées à la base. Ce sont
des fonctions non contractées et elles sont dotées d'un faible exposant leur permettant
de représenter l'espace diffus situé à grande distance du noyau. Une fonction diffuse est
représentée par le signe+. Pour la base 6-31++G(d,p), les fonctions s et p supplémentaires
sont ajoutées pour les atomes lourds alors que pour la base 6-31++G(d,p), des fonctions
diffuses de type s sont ajoutées à l'atome d'hydrogène. Les bases « correlation consistent
polarized » (cc-pVXZ) peuvent être augmentées par des fonctions diffuses servant à décrire
les interactions à plus longue portée, il faut alors ajouter le préfixe aug- au nom de la base.
### 2.11.4 Bases explicitement corrélées

Les méthodes explicitement corrélées (section 2.8.2) ont été conçues pour accélérer la convergence de l'énergie de corrélation en fonction de la taille de base utilisée. Un ensemble de bases a été optimisé afin de les utiliser avec ces méthodes [142, 143]. Appelées cc-pVXZ-F12 (X=D,T,Q), elles sont similaires en taille et en composition aux bases aug-cc-pVXZ et aug-cc-pV(X+D)Z. Les bases explicitement corrélées présentent l'avantage d'améliorer la convergence vers la limite CBS (complete basis set) [142, 143]. Leurs performances ont été évaluées en utilisant la méthode MP2-F12 pour calculer les énergies pour des molécules contenant des atomes appartenant à la première et à la deuxième colonne du tableau périodique. Par exemple, la plus petite base cc-pVDZ-F12 récupère près de 99% de l'énergie de corrélation. La performance de ces nouvelles bases explicitement corrélées a été évaluée par le calcul des propriétés énergétiques de 15 réactions chimiques [142]. Par rapport à la limite CBS, l'erreur maximale trouvée avec cc-pVTZ-F12 est de -0.2 kcal/mol. [142].

### 2.11.5 Bases complètes : CBS

Un des principaux objectifs de la chimie théorique moderne est le calcul précis des propriétés réactionnelles. Ceci a été obtenu pour une grande variété de réactions par l'utilisation d'une base complète CBS développée par Petersson et ses collaborateurs [144, 145, 146]. Cette procédure élimine toute erreur due à la troncation de la base. Les modèles CBS utilisent la convergence asymptotique obtenue en extrapolant des calculs utilisant une base finie. Par exemple, les bases de Dunning ont un comportement de type :

$$E_X = E_{CBS} + B\exp(-CX) \quad (2.130)$$

où $E_{CBS}$ désigne l'énergie extrapolée vers la base infinie et X le nombre cardinal de la base (n=2 pour cc-pVDZ, n=3 pour cc-pVTZ...). L'extrapolation à base infinie est obtenue si l'énergie pour trois n différents est connue.

### 2.11.6 Erreur de superposition de base (Basis set superposition error BSSE)

Lorsque deux espèces, A et B, interagissent ou se rapprochent, les fonctions de base servant à les décrire se superposent : chacune acquiert un espace supplémentaire en prenant les orbitales atomiques de l'autre. Des erreurs de superposition de bases sont donc susceptibles de se produire si on ne prend pas en compte l'espace complet de A et de B. Les énergies de faibles intensités notamment les énergies de Van der Waals sont alors surestimées. Ce phénomène est plus marqué si on utilise des petites bases. Un moyen d'évaluer la correction à apporter à cette erreur est l'utilisation de la méthode de Boys et Bernardi [147], dite méthode de « contrepoids » (Counterpoise Correction). Prenons l'exemple de deux systèmes A et B ayant deux bases a et b, respectivement.
L'énergie d'interaction est donnée par l'expression suivante :

$$ \Delta E(AB) = E_{AB}^{ab} - E_{A}^{a} - E_{B}^{b} $$

où l'exposant indique la base utilisée et l'indice le système considéré. Ainsi, $E_{AB}^{ab}$ représente l'énergie du complexe AB évaluée avec la base ab et calculée pour la géométrie du complexe. Les entités isolées sont décrites par leurs propres bases : $E_{A}^{a}$ est l'énergie de A avec la base a et $E_{B}^{b}$ est l'énergie de B avec la base b.

La correction de contrepoids est ajoutée à l'énergie d'interaction afin de corriger cette erreur :

$$ \Delta E^{CP}(AB) = \Delta E(AB) + \Delta E(CP) $$

La correction est définie par :

$$ \Delta E(CP) = E_{AB}^{ab} - E_{A}^{ab} - E_{B}^{ab} $$

où $E_{A}^{ab}$ est l'énergie de A et $E_{B}^{ab}$ est l'énergie de B, calculées avec la base d'orbitales atomiques ab.

## 2.12 Conclusion

Au cours de ce chapitre, nous avons résumé les méthodes ab initio de la chimie théorique, basées sur la notion de la fonction d'onde, de la plus simple à la plus précise qui permet une meilleure estimation de la corrélation électronique. Nous avons par ailleurs introduit les méthodes de la fonctionnelle de la densité basées sur la densité électronique. Dans notre travail, pour la réaction d'abstraction d'atome d'hydrogène à partir du tétraméthylsilane, nous avons optimisé les structures électroniques des fragments et nous avons calculé les fréquences harmoniques et le chemin d'énergie minimum avec la méthode perturbative, MP2 et avec une base assez étendue, la cc-pVTZ. Afin de prédire des propriétés énergétiques de façon plus précise, la méthode CCSD(T) a été utilisée pour un calcul « single point energy ». Afin d'étudier la dynamique de la réaction, une surface d'énergie potentielle au niveau CCSD(T)/cc-pVTZ//MP2-cc-pVTZ a été construite.

Pour l'étude des interactions des atomes d'hydrogènes avec les nano-clusters de silicates, nous avons effectué des calculs basés sur la DFT. Nous avons choisi la fonctionnelle hybride (MPWB1K) puisqu'elle permet de bien décrire les liaisons de faible interactions [148]. De plus, elle a été aussi utilisée dans plusieurs études en chimie interstellaire [32, 31, 29, 149] (voir Annexe A, paragraphe (A.4.2) pour plus de détails).
Chapitre 3

# Dynamique moléculaire

## 3.1 Introduction

Notre objectif est de décrire aussi précisément que possible les processus élémentaires tels qu’ils se déroulent au niveau microscopique lors d’une réaction chimique et de relier ces résultats aux observations que l’on peut faire au niveau macroscopique (grandeurs thermodynamiques et propriétés cinétiques) et également aux expériences de laboratoire (faisceaux moléculaires croisés et laser) qui se sont développées ces dernières années. Les expériences tel que la fluorescence induite par laser (LIF) et l’ionisation multiphotonique résonante couplée à une analyse par spectrométrie de masse (REMPI) permettent l’étude de l’influence sur les taux de réaction des différents états quantiques des réactants et de comprendre les mécanismes conduisant à la distribution des états quantiques des produits. Certaines expériences permettent aussi de connaître la distribution angulaire des produits (sections différentielles). Il est évident que des théories de plus en plus sophistiquées s’avèrent nécessaires pour interpréter ces expériences.

Du point de vue microscopique et dans le cadre de l'approximation de Born-Oppenheimer, de nombreuses théories ont été développées depuis les travaux pionniers de Hirschfelder et al. [150] et Hulburt et Hirschfelder [151] sur la collision colinéaire H+H₂. Les méthodes (classiques ou quantiques) diffèrent selon que la réaction est sans barrière, avec barrière ou si elle présente des intermédiaires multiples. Les réactions sans barrière sont très rapides et les constantes de vitesse varient de 10⁻⁹ à 10⁻¹¹ cm³ s⁻¹ molec⁻¹. La variation des constantes de vitesse en fonction de la température ne suit pas la loi d'Arrénius et parfois les constantes diminuent lorsque la température augmente. Pour ces réactions, qui sont souvent très exothermiques, la théorie classique de capture peut souvent être appliquée [152].

Les réactions pour lesquelles la surface de potentiel présente plusieurs maxima et minima sont un véritable challenge tant du point de vue du calcul de la PES et du fit de celle-ci que des calculs de dynamique. Même des réactions triatomiques sont alors très
difficiles à étudier. On peut citer la réaction : $O(^1D) + H_2 \rightarrow OH + H$ qui a été l'objet de nombreuses études expérimentales [153, 154] et théoriques [155, 156]. L'étude des mécanismes de réaction montre qu'à basse énergie de collision, la réaction se produit via une insertion de O entre les deux atomes H alors qu'à haute énergie le mécanisme est plutôt un mécanisme d'abstraction [157]. Le même comportement est observé pour la réaction $Si + O_2 \rightarrow SiO + O$ [158]. L'approche QCT (Quasi Classical Theory) peut être une approche intéressante qui donne de bons résultats, cependant elle ne permet pas de tenir compte des effets quantiques importants lorsque des atomes légers (H par exemple) sont concernés.

Ces dernières années, de nombreux travaux théoriques ont été consacrés aux réactions avec barrière, en particulier pour les triatomiques. Les sections différentielles, obtenues par approche quantique, des réactions "benchmark" $H+H_2$ et $H+D_2$ sont en remarquable accord avec les résultats expérimentaux [159]. La réaction $F + \text{ortho}-D_2$ a été l'objet d'un très grand intérêt par suite de l'existence de résonances de Feshbach qui produisent des pics très prononcés dans la section différentielle. L'accord théorie-expérience est exceptionnel [160]. Il faut noter que cet accord repose pour l'étude théorique sur la prise en compte (i) de trois états moléculaires couplés dont les PESs ont été calculées avec une très grande précision, (ii) des couplages non-adiabatiques entre ces surfaces, (iii) et sur étude quantique de la dynamique incluant tous les degrés de liberté.

L'étude des réactions impliquant de plus grosses molécules nécessite un certain nombre d'approximations, tant pour l'obtention de la PES que pour les calculs de dynamique. Ces systèmes présentent plusieurs difficultés : tout d'abord, il faut construire une surface d'énergie potentielle précise à partir de calculs *ab initio* de structure électronique et en obtenir un fit. Ensuite, il faut traiter tous les degrés de liberté des réactants ou produits afin de réaliser les calculs de dynamique. Il y a pour N atomes, 3N-6 degrés de liberté ce qui rend très rapidement les calculs irréalisables. C'est pourquoi de nouvelles approches de dynamique ont été développées. On peut citer par exemple la méthode "multiconfigurational time-dependent Hartree" (MCTDH) qui permet de traiter quantiquement la dynamique pour de nombreux degrés de liberté couplés [161]. Récemment, une méthode a été développée qui demande la construction d'une surface de potentiel à dimensionnalité réduite avec un nombre relativement limité de points *ab initio*, suivie par des calculs de dynamique quantique qui traitent explicitement les liaisons qui se brisent et qui se forment [162]. C'est cette méthode 2D à dimensionnalité réduite que nous avons utilisée.

Une première partie de ce chapitre (3.2) est dédiée à la description de la PES d'une collision réactive ainsi que des points caractéristiques de cette surface. Notamment, on examinera le chemin réactionnel intrinsèque et son rôle dans le calcul des constantes de vitesse dans le cadre de la théorie de l'état de transition. Les probabilités semi-classiques déterminées à partir de ces données sont décrites dans le paragraphe (3.3).

La surface d'énergie potentielle peut être utilisée dans un traitement quantique de la dynamique nucléaire à travers un modèle à dimensionnalité réduite qui s'applique à
des réactions chimiques ayant une approche colinéaire. Dans ce modèle, la méthode de dynamique quantique indépendante du temps permet de traiter explicitement les liaisons chimiques qui se forment et se brisent entre les atomes dans le cadre de la théorie de la diffusion. Dans ce chapitre, nous décrivons ce modèle déjà développé pour étudier les problèmes de collisions réactives. Nous commençons par définir le système de coordonnées et l'hamiltonien utilisés. Nous présentons ensuite la méthode de propagation de la R-matrice empolyée pour résoudre les équations couplées permettant d'extraire la matrice S de diffusion et d'en déduire des probabilités d'état-à-état.

## 3.2 Surface d'énergie potentielle (PES)

Dans le cadre de l'approximation de Born-Oppenheimer, la détermination des valeurs propres de l'équation de Schrödinger électronique permet d'obtenir une surface d'énergie potentielle (PES). La construction de cette surface est indispensable à la résolution de l'équation de Schrödinger nucléaire qui permet la caractérisation des mouvements des noyaux lors d'une collision. A une configuration géométrique donnée de l'ensemble des noyaux correspond un point représentatif sur la PES ayant une énergie notée V. Parmi ces structures, il y a des points stationnaires correspondant à un minimum de l'énergie potentielle tel que les composantes du vecteur gradient sont nulles :

$$F_i(x_{eq}) = \left( \frac{\partial V(x)}{\partial q_i} \right)_{x_{eq}} = 0 \qquad (3.1)$$

où $q_i$ est le déplacement interne à partir d'une position initiale $x_0$ tel que $q_i=x_i-x_0$. $x_{eq}$ représente les coordonnées du minimum.

La PES permet également la détermination de tous les chemins reliant les réactants (Q) aux produits (P). Ces points stationnaires sont représentés sur la figure 3.1 par deux vallées (correspondant aux structures d'équilibre). Ce sont des points à partir desquels tout mouvement dans une direction quelconque de la PES, associé à une déformation structurale donnée de la molécule, conduit nécessairement à une valeur supérieure de l'énergie. En partant de la première vallée vers la seconde, certains chemins peuvent franchir un "col". Ces chemins sont identifiés comme "chemins réactionnels". Le point selle (col) correspond à un maximum local d'énergie dans une direction particulière de la PES et à un minimum dans toutes les autres directions. Un tel point correspond à une géométrie (ou état) de transition de la molécule, reliant deux structures à l'équilibre. Ces chemins réactionnels permettent de représenter le système moléculaire en évolution et dépendent donc de la répartition de l'énergie initiale des réactifs sur les degrés de liberté. Afin de déterminer la nature de ces points stationnaires, une analyse des modes normaux de vibration est indispensable.

### 3.2.1 Modes normaux de vibration

Une optimisation géométrique ne révèle pas la nature du point stationnaire. Il faut, pour identifier celle-ci, calculer les fréquences caractéristiques de la géométrie optimisée.
<figure><img src="image_7.png" /><figcaption>FIGURE 3.1 – Surface d'énergie potentielle du système A + BC → AB + C.</figcaption></figure>

Ce calcul est rendu possible grâce aux programmes de structures électroniques tels que Gaussian [163], Molpro [164], Gamess [165], Nwchem [166], etc.

Le calcul des fréquences permet de représenter les mouvements des atomes dans le système considéré, selon ses degrés de liberté. Pour une molécule composée de N atomes, le nombre de degrés de liberté est égal à 3N. Parmi ces degrés de liberté, trois correspondent à une translation de toute la molécule et trois sont les rotations de la molécule (2 pour une molécule linéaire). Les 3N-6(5) degrés de liberté restant correspondent aux déformations structurales de la molécule (longueurs et angles de liaison), que l'on nomme vibrations internes. Ces mouvements internes se décomposent en modes normaux de vibration et correspondent à des fréquences propres de vibration. Ces modes de vibration sont déterminés à partir de la diagonalisation, pour chaque point stationnaire, de la matrice Hessienne, H. Cette matrice est formée des dérivées partielles secondes du potentiel par rapport au coordonnées cartésiennes pondérées par les masses atomiques :

$$H_{ij} = \left( \frac{\partial^2 V}{Q_i Q_j} \right)_{x_{eq}} \qquad (3.2)$$

où $Q_i = \sqrt{m_i} \times q_i$.

Les vecteurs propres de H sont les modes normaux de vibration et les valeurs propres correspondantes, $\epsilon_i$, sont reliées aux fréquences propres :

$$\epsilon_i = \omega_i^2 \qquad (3.3)$$

Les informations extraites de la diagonalisation de la hessienne permettent de distinguer différents points caractéristiques sur une PES :

- Si les valeurs propres de la matrice hessienne sont toutes positives, alors les fréquences
sont toutes réelles et le point stationnaire est un minimum local.

- Si l'une d'entre elles est négative, c'est-à-dire correspond à une fréquence imaginaire, le point stationnaire correspond alors à un état de transition (TS).

- Si plus d'une d'entre elles est négative, c'est un point selle d'ordre supérieur.

### L'énergie du point zéro (ZPE)

Dans l'approximation de l'oscillateur harmonique, l'énergie quantifiée de vibration peut être écrite selon l'expression suivante :

$$V_n = (n + \frac{1}{2})\hbar\omega \quad (3.4)$$

où $\omega$ la fréquence vibrationnelle et n le nombre quantique de vibration. L'énergie fondamentale de vibration correspond alors à $\hbar\omega/2$. Cette énergie la plus basse est appelée énergie du point zéro (ZPE). L'origine de la ZPE se trouve dans le principe d'incertitude de Heisenberg. Si la position d'une particule est connue, sa quantité de mouvement ne peut être déterminée précisément et donc son énergie cinétique ne peut pas être nulle.

Lors de l'optimisation de la géométrie d'un point stationnaire de la surface d'énergie potentielle, on obtient une énergie potentielle qui correspond à un minimum de la courbe d'énergie potentielle à laquelle on ajoute la ZPE de tous les modes vibrationnels, à leurs positions d'équilibre respectifs, afin de tenir compte de cette fluctuation quantique.

### 3.2.2 Chemin d'énergie minimum (MEP)

Après avoir optimisé les géométries des points stationnaires, le calcul des fréquences permet de déterminer leur nature et notamment le point selle. Cependant, il convient de vérifier que l'état de transition obtenu est bel et bien celui que l'on recherchait à l'origine et non pas celui d'un autre chemin réactionnel. Pour cela, un calcul d'énergie point par point est effectué en définissant des variations de faible amplitude à partir de la géométrie du TS selon ce mode de déformation. Un chemin d'énergie minimum (MEP) est alors localisé en suivant le chemin de la plus grande pente caractérisé par une valeur maximale du gradient. Qualitativement, cette action permet, à partir de la géométrie de l'état de transition, de faire varier de part et d'autre, les coordonnées du système afin de trouver le chemin d'énergie minimum (MEP) qui connecte les deux vallées correspondant aux réactants et aux produits. On appellera $s$ la coordonnée sur ce chemin de réaction. L'approche numérique appelée « coordonnée intrinsèque de réaction » (IRC) [167] permet à partir du TS de déterminer la MEP, où les coordonnées cartésiennes sont "mass-weighted" i.e. pondérées par la racine carrée des masses atomiques.

La détermination de la MEP est indispensable pour le calcul des taux de réaction avec des méthodes semi-classiques qui tiennent en partie compte des effets quantiques (effet tunnel).
<figure><img src="image_8.png" /><figcaption>FIGURE 3.2 – Barrière de potentiel à une dimension.</figcaption></figure>

## 3.3 Probabilités semi-classiques

Les atomes légers tels que l'hydrogène peuvent traverser les barrières de potentiel par effet tunnel, ce qui ne peut s'effectuer classiquement. Cet effet quantique est important aux basses températures. C'est l'effet dominant pour la description des probabilités lors de réactions dominées par l'existence d'un état de transition (TST, voir chapitre 4).

On se propose d'étudier le comportement d'une particule de masse *m* et d'énergie *E* quand elle traverse une barrière carrée de potentiel (voir figure 3.2) telle que :

$$ \begin{cases} V = V(x) & \text{pour } x_1 \le x \le x_2 \rightarrow \text{zone 1} \\ V = 0 & x < x_1 \text{ et } x > x_2 \rightarrow \text{zone 2} \end{cases} $$

Ce problème peut se traiter exactement en imposant la continuité de la fonction d'onde et de sa dérivée à chaque discontinuité du potentiel. On ne donnera ici que les résultats principaux (voir [168] et [169] pour plus de détails).
L'équation de Schrödinger indépendante du temps s'écrit :

$$ \left[ \frac{d^2}{dx^2} + \frac{2m}{\hbar^2} (E - V(x)) \right] \Psi(x) = 0 \qquad (3.5) $$

Soit :

$$ k = \sqrt{\frac{2mE}{\hbar^2}} \qquad \mathcal{K} = \sqrt{\frac{2m(E-V)}{\hbar^2}} \qquad (3.6) $$

- Pour la zone 2, la fonction d'onde s'écrit :

$$ x < x_1 \rightarrow \Psi(x) = A e^{ikx} + B e^{-ikx} \qquad (3.7) $$

$$ x > x_2 \rightarrow \Psi(x) = F e^{ikx} + G e^{-ikx} \qquad (3.8) $$

- Pour la zone 1, on cherche $\Psi(x)$ sous la forme :

$$ \Psi(x) = e^{(i/\hbar)\sigma(x)} \qquad (3.9) $$
L'équation (3.5) devient :

$$ \frac{1}{2m} \left(\frac{d\sigma}{dx}\right)^2 - \frac{i\hbar}{2m} \frac{d^2}{dx^2}\sigma = E - V(x) \qquad (3.10) $$

$$ \frac{1}{2m}(\sigma')^2 - \frac{i\hbar}{2m}\sigma'' = E - V(x) \qquad (3.11) $$

Supposons que :

$$ \sigma''/(\sigma')^2 \ll 1 \qquad (3.12) $$

l'équation (3.11) devient :

$$ \frac{1}{2m}(\sigma')^2 = E - V(x) \qquad (3.13) $$

Si $E > V(\mathbf{x})$, on a alors :

$$ \sigma(x) = \pm \int \sqrt{2m(E-V)} dx \qquad (3.14) $$

La fonction d'onde est de la forme :

$$ \Psi(x) = C e^{i\mathcal{K}x} + D e^{-i\mathcal{K}x} \qquad (3.15) $$

Si $E < V(\mathbf{x})$, la région $x_1 \le x \le x_2$ est classiquement interdite, le passage se produit par effet quantique. L'équation (3.15) est encore valable à condition de remplacer :

$$ \mathcal{K} = i\kappa $$

avec

$$ \kappa = \sqrt{\frac{2m(V-E)}{\hbar^2}} $$

En introduisant les relations de continuité de la fonction d'onde et de sa dérivée en $x_1$ et $x_2$, on obtient des relations entre les coefficients $A, B, C, D, F$ et $G$. Ceci permet de déduire l'amplitude de la fonction d'onde après la traversée de la barrière en fonction de l'amplitude initiale [168]. On obtient finalement la probabilité semi-classique $P(E)$, à partir du carré de la fonction d'onde transmise :

$$ P(E) = \exp[-2\theta(E)] \qquad (3.16) $$

$$ \text{avec } \theta(E) = \sqrt{\frac{2m(V-E)}{\hbar^2}}(x_2 - x_1) \qquad (3.17) $$

Ce résultat se généralise facilement au cas d'une barrière de potentiel de forme quelconque en :

$$ \theta(E) = \int \sqrt{\frac{2m(V(x) - E)}{\hbar^2}} dx \qquad (3.18) $$

On remarque ainsi que les paramètres physiques importants pour le coefficient de transmission sont l'énergie et la hauteur de la barrière, mais également sa forme, notamment
sa largeur. On retrouve ce résultat important dans la théorie de l'état de transition.

**Remarque** : L'approximation (équation 3.12) peut se réécrire sous la forme :

$$ \frac{1}{2\pi} \frac{d\lambda}{dx} \ll 1 \qquad (3.19) $$

où $\lambda = \frac{2\pi}{K}$ est la longueur d'onde de de Broglie de la particule dans la zone $x_1 \le x \le x_2$. Cette approximation impose que la longueur d'onde de de Broglie ne varie que faiblement sur des distances de l'ordre de la longueur d'onde elle-même. C'est l'approximation « quasi classique » ou « semi-classique ».

## 3.4 Méthodes de dynamique quantique : Modèle à dimensionnalité réduite (2D)

L'objectif de la méthode de dynamique quantique est d'étudier les mécanismes réactionnels à l'échelle moléculaire. Les collisions du type $H_a + H_b YZ \rightarrow H_a H_b + YZ$, ayant un chemin réactionnel colinéaire, peuvent être étudiées en appliquant un modèle à dimensionnalité réduite à deux degrés de liberté (2D).

Nous allons d'abord appliquer ce modèle à une réaction triatomique en spécifiant les coordonnées hypersphériques utilisées pour décrire le processus réactionnel (section 3.4.1). L'extension de ce modèle à des réactions plus complexes est présentée dans le paragraphe qui suit (3.4.2). La construction de la PES est détaillée dans le paragraphe (3.4.3). La résolution des équations couplées en utilisant le formalisme de la R-matrice [170] est donnée dans les paragraphes (3.4.4) et (3.4.5). L'expression de la matrice de diffusion et le calcul des sections efficaces sont décrits dans les paragraphes (3.4.6) et (3.4.7). Quelques notions fondamentales de la théorie élémentaire de la diffusion sont introduites dans l'annexe B.

### 3.4.1 Cas d'une réaction triatomique : A + BC → AB + C

Le problème de la collision réactive pour un système triatomique est relativement simple à formuler. La résolution se fait en utilisant une approche appelée « coupled channel (CC) »[171]. L'approche générale à adopter pour résoudre le problème de diffusion quantique pour ce type de collisions se base sur :

- Un développement de la fonction d'onde totale sur une base appropriée,

- Une propagation des équations couplées,

- Une extraction de la matrice de diffusion S en appliquant les conditions limites aux régions asymptotiques.

Dans le cas du traitement exact du processus collisionnel, le nombre d'équations couplées à propager devient important lorsque tous les états vibrationnels et rotationnels énergétiquement accessibles doivent être considérés. Ce nombre peut varier de quelques
centaines à plusieurs milliers [171]. La majorité de ces états sont les états rotationnels (en particulier la dégénérescence des états projetés, $J_z$) associés aux états vibrationnels.

En raison des limites de calcul, des approximations ont été développées. La première, appelée approximation des états couplés « coupled states approximation (CS) » [171], suppose que les différents états projetés peuvent être découplés. Un autre modèle, nommé « reduced dimensionality exact quantum (RDEQ) theory » [171] se basant sur l'approche CS, introduit d'autres approximations pour traiter les mouvements de rotation et de vibration pour les réactions triatomiques de type A + BC → AB + C. Une des approximations adoptées est celle qui suppose que la réaction suit un chemin colinéaire : c'est le modèle « rotating linear model » (RLM) [171, 172]. Le mode de pliage est supposé lent par rapport aux autres degrés de liberté internes et par conséquent, il est traité adiabaticquement dans les calculs de diffusion. Dans ces conditions, il n'y a plus que deux degrés de liberté à traiter. En supposant de plus que le système moléculaire est un rotateur rigide linéaire (RLM), une expression analytique de l'opérateur moment cinétique total est obtenue. L'expression de l'hamiltonien du modèle RLM est donnée en coordonnées hypersphériques (définies à la fin de ce paragraphe) par :

$$ \hat{H} = -\frac{\hbar^2}{2\mu} \frac{\partial^2}{\partial \rho^2} - \frac{\hbar^2}{2\mu\rho^2} \frac{\partial^2}{\partial \delta^2} + \frac{3\hbar^2}{8\mu\rho^2} + \frac{J(J+1)}{2\mu\rho^2} + V(\rho, \delta) \quad (3.20) $$

Le terme $J(J+1)$ correspond au mouvement relatif quantifié pour un moment cinétique total $\hat{J}$. $V(\rho, \delta)$ est la surface d'énergie potentielle de la réaction colinéaire A + BC. Cette surface est bidimensionnelle ($V(\rho, \delta) = V_{2D}(\rho, \delta)$) puisqu'elle décrit explicitement le mouvement des deux degrés de liberté internes. Le système moléculaire est représenté donc par :

$$ A + BC(s') \rightarrow AB(s) + C $$

où $s'$ et $s$ représentent les vibrations internes initiales et finales des liaisons BC et AB, respectivement.

Cependant, le modèle RLM est généralement insuffisant pour prédire les propriétés dynamiques car il néglige la contribution des degrés de liberté restants [171].

Pour faire face aux limites du modèle RLM, Walker et Hayes [172] ont développé le modèle *Bending Corrected Rotating Linear Model (BCRLM)*. Dans cette approche, les liaisons qui se forment et se cassent sont traitées explicitement dans les calculs de diffusion, alors que le degré de liberté restant (mode de pliage) est pris en compte au moyen de l'énergie du point zéro :

$$ V(\rho, \delta) = V_{2D}(\rho, \delta) + ZPE_{pliage}(\rho, \delta) $$

Les calculs de diffusion quantiques rencontrent plusieurs difficultés en particulier l'augmentation de la taille des matrices à propager et à diagonaliser. L'adoption du modèle à dimensionnalité réduite a permis de simplifier l'étude de ce problème. De plus, il faut choisir le même système de coordonnées permettant de décrire tous les canaux des réactants et des produits à la fois. Ceci est possible grâce aux coordonnées hypersphériques [173].
<figure><img src="image_9.png" /><figcaption>FIGURE 3.3 – Système de coordonnées de Jacobi pour la réaction colinéaire $A + BC \rightarrow AB + C$. « × » dénote l’emplacement des centres de masse de AB et BC.</figcaption></figure>

## Système de coordonnées de Jacobi et coordonnées hypersphériques

Chaque point de la PES correspond à une configuration nucléaire donnée. Les coordonnées hypersphériques sont utilisées pour décrire chaque point de cette surface. Ce système de coordonnées présente l'avantage de tenir compte de tous les réarrangements moléculaires des réactants et des produits. Ces coordonnées peuvent être déduites des coordonnées de Jacobi qui sont utiles pour décrire les mouvements de grandes amplitudes rencontrés au cours de réactions chimiques (par exemple, ceux d'isomérisation) [174].

Pour la réaction $A + BC \rightarrow AB + C$, les coordonnées de Jacobi des réactants correspondent à la distance $r_{BC}$ entre B et C et à la distance $R_A$ entre A et le centre de masse de BC. Rappelons que cette définition n'est valable que pour une réaction colinéaire. La figure 3.3 présente l'ensemble de ces coordonnées. Pour les grandes distances atome-diatomique, $R_A$ décrit le mouvement translationnel de A par rapport à BC alors que les faibles distances $r_{BC}$ décrivent le mouvement vibrationnel de BC. Les coordonnées de Jacobi des produits, $R_C$ et $r_{AB}$, sont définies de manière similaire, et décrivent la situation où l'atome C et la molécule AB sont éloignés.

Afin de transformer les coordonnées de Jacobi en une représentation qui permet une description uniforme des canaux asymptotiques des réactants et des produits, [175] elles
doivent tout d'abord être pondérées selon : [172]

$$R_r = \left[ \frac{M_A}{\mu} \right]^{(1/2)} R_A \qquad (3.21)$$

$$r_r = \left[ \frac{M_{BC}}{\mu} \right]^{(1/2)} r_{BC} \qquad (3.22)$$

$$R_p = \left[ \frac{M_C}{\mu} \right]^{(1/2)} R_C \qquad (3.23)$$

$$r_p = \left[ \frac{M_{AB}}{\mu} \right]^{(1/2)} r_{AB} \qquad (3.24)$$

(3.25)

Les indices r et p se rapportent aux reactants et aux produits, respectivement. La masse réduite $\mu$ est définie par :

$$\mu = \left[ \frac{m_A m_B m_C}{m_{ABC}} \right]^{(1/2)} \qquad (3.26)$$

Les facteurs de masse sont définis comme suit : [175]

$$M_A = \frac{m_A(m_B + m_C)}{m_A + m_B + m_C} \qquad (3.27)$$

$$M_{BC} = \frac{m_B m_C}{m_{BC}} \qquad (3.28)$$

$$M_C = \frac{m_C m_{AB}}{m_{ABC}} \qquad (3.29)$$

$$M_{AB} = \frac{m_A m_B}{m_{AB}} \qquad (3.30)$$

où $m_A$, $m_B$ et $m_C$ sont les masses atomiques.

La pondération des coordonnées de Jacobi permet de réduire la complexité de traitement de ces coordonnées en décrivant chaque type de mouvement avec la même masse réduite, $\mu$. Ces coordonnées sont liées à travers la matrice de rotation suivante : [176]

$$\begin{pmatrix} R_r \\ r_r \end{pmatrix} = \begin{pmatrix} \cos\theta & \sin\theta \\ \sin\theta & -\cos\theta \end{pmatrix} \begin{pmatrix} R_p \\ r_p \end{pmatrix} \qquad (3.31)$$

$\theta$ est le " skewing angle" défini comme suit :

$$\theta = \tan^{-1} \left[ \frac{m_B m_{ABC}}{m_A m_C} \right] \qquad (3.32)$$

La figure 3.3 montre que les coordonnées de Jacobi ne sont appliquées que pour des régions limitées de l'espace. En effet, elles permettent de décrire d'une manière indépendante les canaux asymptotiques des réactifs et des produits respectivement en fonction des coordonnées ($R_r$, $r_r$) et ($R_p$, $r_p$). Pour pouvoir décrire simultanément tout réarrangement
<figure><img src="image_10.png" /></figure>

FIGURE 3.4 – Les coordonnées hypersphériques ($\rho, \delta$) définies à partir des coordonnées de Jacobi des réactants ($R_r, r_r$) et des produits ($R_p, r_p$) pour la réaction $A + BC \rightarrow AB + C$. Deux « contours plots » (en bleu) de la surface d'énergie potentielle aux faibles valeurs de $\rho$ (région entourant l'état de transition) et aux grandes valeurs de $\rho$ (puits de potentiel des réactants $A + BC$ et des produits $AB + C$).

moléculaire dans les équations de dynamique réactionnelle, les coordonnées hypersphériques sont utilisées. Ainsi, les systèmes de coordonnées des réactifs et des produits sont liées dans la région d'interaction.

Les coordonnées hypersphériques (figure 3.4) peuvent être définies en fonction des coordonnées de Jacobi pondérées. Le rayon hypersphérique $\rho$, qui représente la coordonnée de diffusion, est indépendant de l'indice du réarrangement i (où i correspond à r ou p).

$$ \rho = (R_i^2 + r_i^2)^{1/2} \qquad (3.33) $$

alors que l'hyperangle $\delta_i$ (qui permet de décrire la vibration interne) en dépend

$$ \delta_i = \tan^{-1}[\frac{r_i}{R_i}] \qquad (3.34) $$

Les hyperangles $\delta_r$ et $\delta_p$ sont liés au skewing angle $\theta$ via cette relation : [176]

$$ \delta_r + \delta_p = \delta_{max} = \theta \qquad (3.35) $$

Il est donc possible de remplacer $\delta_r$ et $\delta_p$ par une seule coordonnée angulaire, l'hyperangle $\delta$, qui varie de 0 à $\delta_{max}$. Ainsi, un seul ensemble de coordonnées de Jacobi ($R_r, r_r$), par exemple, suffit à définir les coordonnées hypersphériques ($\rho, \delta$). La transition des coordonnées de Jacobi pondérées des réactants à ceux des produits se fait à travers une rotation d'angle $\delta_{max}$.
## 3.4.2 Cas d'un système polyatomique quelconque

L'extension du modèle à dimensionnalité réduite aux collisions polyatomiques quelconques nécessite l'application d'autres approximations dans lesquelles le nombre de degrés de liberté à traiter explicitement dans le calcul de dynamique réactionnelle (mode réactif) est varié. Les modes spectateurs sont introduits dans le potentiel effectif via leurs énergies du point zéro, ZPEs.

La complexité mathématique pour un système polyatomique possédant 3N-6 degrés de liberté augmente avec N et un traitement quantique exact devient alors quasi-impossible. En conséquence, on utilise une extension de l'approche BCRLM à tous les modes non pris en compte explicitement dans le modèle 2D. Cette approximation revient à réduire ce système complexe à un processus de collision pseudo-triatomique approprié de type $H_a+H_b$YZ, Y est l'atome à partir duquel l'atome d'hydrogène ($H_b$) va être éliminé. Dans cette approche, les liens chimiques qui se forment ($H_a-H_b$) et se cassent ($H_b-Y$) sont traités explicitement dans le calcul de dynamique alors que les autres modes sont traités comme spectateurs. Ainsi, la ZPE des 3N-8 modes spectateurs est incluse dans le potentiel effectif $V(\rho, \delta)$ :

$$V(\rho, \delta) = V_{2D}(\rho, \delta) + ZPE_{spectateurs}(\rho, \delta) \quad (3.36)$$

L'expression (3.20) du hamiltonien est alors utilisée en considérant que le potentiel $V(\rho, \delta)$ correspond au potentiel effectif (équation 3.36).

## Système de coordonnées de Jacobi et coordonnées hypersphériques

L'extension des coordonnées hypersphériques (paragraphe précédent) à la réaction $H_a + H_b$YZ → $H_a$H$_b$ + YZ (figure 3.5) se fait par l'adaptation des définitions des coordonnées de Jacobi. Les coordonnées de Jacobi décrivant les produits correspondent à la distance r entre $H_a$ et $H_b$ et à la distance R entre le centre de masse de $H_a$H$_b$ et Y (voir figure 3.5). Dans notre approche, le centre de masse de YZ, la molécule à partir de laquelle l'hydrogène est éliminé, est assimilé à Y. Cette approximation a été adoptée pour l'étude de plusieurs réactions d'abstraction à partir de CH$_3$OH, C$_2$H$_6$, C$_3$H$_8$ [177, 178, 179, 180, 181], CH$_3$NH$_2$ [179] et n-C$_4$H$_{10}$ [182]. Cette approximation, qui permet de réduire la complexité du calcul du centre de masse, a conduit dans tous les systèmes étudiés à des taux de réaction en bon accord avec l'expérience [182].

Les facteurs de masse utilisés pour définir les coordonnées hypersphériques ainsi que le « skewing angle » sont :

$$M_1 = \frac{(m_Y + m_Z) \times (m_{H_a} + m_{H_b})}{m_{H_a} + m_{H_b} + m_Y + m_Z} \quad (3.37)$$

$$M_2 = \frac{m_{H_a} \times m_{H_b}}{m_{H_a} + m_{H_b}} \quad (3.38)$$

$$M_3 = \frac{m_Y \times m_Z}{m_Y + m_Z} \quad (3.39)$$

$$\mu = [M_1 M_2 M_3]^{(1/2)} \quad (3.40)$$
<figure><img src="image_11.png" /><figcaption>FIGURE 3.5 – Système de coordonnées colinéaires pour la réaction H_a + H_b YZ → H_a H_b + YZ.</figcaption></figure>

où m_H, m_Y et m_Z sont les masses atomiques.

Les coordonnées hypersphériques sont obtenues à partir de la transformation suivante [175] :

$$ \left(\frac{M_1}{\mu}\right)R^2 = [\rho\cos(\delta)]^2 \qquad (3.41) $$

$$ \left(\frac{M_2}{\mu}\right)r^2 = [\rho\sin(\delta)]^2 \qquad (3.42) $$

Notons que cette approche quantique à dimensionnalité réduite 2D ne s'applique qu'aux réactions colinéaires. Les énergies électroniques des points de la PES et les énergies du point zéro des modes spectateurs sont obtenues par des calculs *ab initio*. La représentation de potentiel effectif par une fonctionnelle analytique nécessite un fit adapté aux caractéristiques de la surface (voir paragraphe suivant).

### 3.4.3 Construction de la PES à dimensionnalité réduite (RD)

L'approche à dimensionnalité réduite [183, 184, 185] adoptée dans notre travail pose sur l'optimisation de la géométrie du système réactif, le calcul de l'énergie avec des méthodes *ab initio* performantes et un calcul de fréquence pour chaque point ρ, δ de la surface. La PES est ensuite obtenue à partir d'un ajustement (fit) à une fonctionnelle analytique construite à l'aide de fonctions de Morse, assurant ainsi une réduction du nombre des points *ab initio* à calculer. Les calculs *ab initio* sont effectués de façon à caractériser les régions les plus importantes de la surface, sans pour autant augmenter de manière significative le coût du calcul. Néanmoins, ce fit doit être précis afin que le traitement dynamique quantique soit efficace.

La méthode que nous avons utilisée pour construire la surface d'énergie potentielle (2D) est décrite ci-dessous :
Optimisation de la géométrie

La surface d'énergie potentielle (PES) est calculée à partir d'une grille de points (R, r) en
utilisant des programmes de chimie quantique tels que Gaussian09 [163] et Molpro [164].
L'optimisation directe des géométries en coordonnées hypersphériques n'est pas possible
avec ces codes. A chaque point de la grille (R,r), la géométrie est optimisée en figeant
les modes actifs et en permettant aux autres degrés de liberté d'être totalement relaxés.
L'optimisation est donc faite sous contrainte.

Dans notre étude d'abstraction d'un atome d'hydrogène à partir du tétraméthysilane, l'optimisation est faite au niveau MP2 (voir paragraphe 2.9 pour plus de détails sur la méthode de calcul), avec une base assez étendue, cc-pVTZ par exemple. Il faut s'assurer que tous les domaines la PES sont bien représentés par les points *ab initio*.

Calcul de l'énergie : $V_{2D}$

Afin d'atteindre un degré de précision plus important, le calcul du potentiel est effectué
pour chaque point de la grille (R,r) ( ou « single point energy ») avec des méthodes *ab
initio* plus sophistiquées.

Les méthodes des clusters couplés (voir 2.8 pour plus de détails sur la méthode de calcul) avec des bases étendues sont les plus utilisées pour les calculs de précision puisqu'elles prennent en compte un grand pourcentage de la corrélation électronique. La méthode CCSD(T) est le plus souvent utilisée pour ce calcul des énergies des points *ab initio* servant à construire la surface bidimensionnelle. Dans notre travail, nous avons effectué des calculs au niveau CCSD(T)/cc-pVTZ.

Des travaux récents [181, 186] ont mis en évidence la performance des méthodes
explicitement corrélées telles que CCSD(T)-F12a (voir paragraphe 2.8.2 pour plus de
détails sur la méthode de calcul) et qui permettent de calculer des propriétés énergétiques
très précises avec des bases de petites tailles, réduisant ainsi l'effort numérique. Nous
avons testé ce type de base dans l'étude de la réaction d'abstraction d'atomes H à partir
du tétraméthylsilane (voir chapitre 6, paragraphe 6.2).

Calcul de l'énergie de point zéro des modes spectateurs : $ZPE_{spectateurs}$

Le calcul des fréquences doit être effectué avec la même méthode et la même base que celles utilisées pour l'optimisation de la géométrie (MP2/cc-pVTZ). Pour chaque point de la grille (R, r), on extrait de la matrice hessienne la contribution des modes réactifs afin de ne conserver que la contribution des modes spectateurs à l'énergie du point zéro. La méthode de projection utilisée pour obtenir l'énergie de point zero $ZPE_{spectateurs}$ est décrite ci-après.

Méthode de projection rectiligne des modes réactifs

L'analyse vibrationnelle ne se fait que pour les points stationnaires ayant un gradient d'énergie nul. Cependant, les géométries optimisées sous contrainte pour construire la surface de dimensionnalité réduite possèdent un gradient non nul par rapport aux deux modes actifs traités explicitement dans les calculs de diffusion. Leurs mouvements vibra-
tionnels peuvent être mélangés avec les autres modes.

Pour pouvoir déterminer l'énergie du point zéro des modes spectateurs et ainsi calculer $V(\rho, \delta)$ (équation 3.36), les modes actifs doivent être découplés des premiers avant la diagonalisation de la hessienne $H$ (équation (3.2)).

Ce formalisme a été utilisé dans des travaux antérieurs pour les modes translationnels et rotationnels [187] et pour des modes réactifs en appliquant une contrainte via un opérateur de projection [188]. Les coordonnées utilisées dans notre travail sont les coordonnées rectilignes.

Soit un sous-espace $\mathbf{C}$ contenant $d$ modes actifs et un sous-espace complémentaire $\mathbf{C}1$ contenant $F-d$ modes spectateurs tel que $F=3N-6$. Ces deux sous-espaces appartiennent à un même espace nommé $\mathbf{R}$. La détermination de la ZPE des modes spectateurs nécessite la projection de la hessienne de dimension $3N-6$ de l'espace $\mathbf{R}$ sur le sous-espace $\mathbf{C}1$ de dimension $F-d$.

Tout vecteur $\mathbf{u}$ appartenant à $\mathbf{R}$ aura une projection dans $\mathbf{C}$ et dans $\mathbf{C}1$ à travers un opérateur de projection $\mathbf{P}^C$ et $\mathbf{P}^{C1}$ tel que :

$$ \mathbf{u} = \mathbf{u}^C + \mathbf{u}^{C1} \qquad (3.43) $$

$$ \mathbf{u}^C = \mathbf{P}^C \mathbf{u} \qquad (3.44) $$

$$ \mathbf{u}^{C1} = \mathbf{P}^{C1} \mathbf{u} \qquad (3.45) $$

Les opérateurs de projection possèdent les propriétés suivantes :

$$ (\mathbf{P}^C)^2 = \mathbf{P}^C \qquad (3.46) $$

$$ \mathbf{P}^C + \mathbf{P}^{C1} = 1 \qquad (3.47) $$

$$ \mathbf{P}^C \mathbf{P}^{C1} = 0 \qquad (3.48) $$

On aura alors

$$ \mathbf{P}^{C1} = 1 - \mathbf{P}^C \qquad (3.49) $$

Ainsi, dans le sous espace $\mathbf{C}1$ de dimension $F-d$, la hessienne $H^P$ qui doit être diagonalisée est donnée par : [188]

$$ H^P = (1 - \mathbf{P}^C)H(1 - \mathbf{P}^C) \qquad (3.50) $$

$\mathbf{P}^C$ est l'opérateur de projection relatif aux élongations des liaisons traitées explicitement dans les calculs de diffusion (les mouvements le long de $H_a-H_b$ et de $H_b$-YZ). Il est construit dans un système de coordonnées cartésiennes pondérées par les masses atomiques du système moléculaire étudié : [188]

$$ \mathbf{P}^C = e_{(H_a H_b)} e_{(H_a H_b)}^\top + e_{(H_b - YZ)} e_{(H_b - YZ)}^\top \qquad (3.51) $$

Les vecteurs unitaires orthogonaux $e_{(H_a H_b)}^\top$ et $e_{(H_b - YZ)}^\top$ sont les transposés des vecteurs $e_{(H_a H_b)}$ et $e_{(H_b - YZ)}$. Ces derniers sont construits individuellement à partir de l'orthogonálisation par la procédure standard de Gram-Schmidt [189] des vecteurs unitaires $e_{(H_a H_b)}^0$
et $e^0_{(H_b-YZ)}$ le long de $H_a-H_b$ et $H_b-YZ$. Ils possèdent chacun seulement 6 éléments non nuls relatifs aux coordonnées cartésiennes de $H_a$, $H_b$ et de $H_b$ et $YZ$, respectivement.

Le vecteur $e^0_{(H_a H_b)}$ a pour expression : [188]

$$e^0_{(H_a H_b)} = \begin{pmatrix} 0 \\ \vdots \\ R^{(H_a H_b)}_{H_a, x} \\ R^{(H_a H_b)}_{H_a, y} \\ R^{(H_a H_b)}_{H_a, z} \\ 0 \\ \vdots \\ R^{(H_a H_b)}_{H_b, x} \\ R^{(H_a H_b)}_{H_b, y} \\ R^{(H_a H_b)}_{H_b, z} \end{pmatrix} \quad (3.52)$$

Les six composantes non nulles sont définies comme suit : [188]

$$R_{H_{a,i}} = \frac{N}{d_{(H_a H_b)}} \left[ \left( \frac{m_{H_b}}{m_{H_a}} \right)^{1/2} (x_{(H_a,i)} - x_{(H_b,i)}) \right] \quad (3.53)$$

$$R_{H_{b,i}} = \frac{N}{d_{(H_a H_b)}} \left[ \left( \frac{m_{H_a}}{m_{H_b}} \right)^{1/2} (x_{(H_b,i)} - x_{(H_a,i)}) \right] \quad (3.54)$$

où N est une constante de normalisation, $d_{(H_a H_b)}$ est la distance interatomique entre $H_a$ et $H_b$ et $x_{(H_a,i)}$ et $x_{(H_b,i)}$ sont les coordonnées cartésiennes pondérées telles que $i=x,y,z$.

La hessienne $H^P$ est ainsi déterminée à partir des équations (3.50) et (3.51). La diagonalisation de cette matrice donne seulement 3N-8 valeurs propres relatives aux modes spectateurs. Les modes translationnels (3 modes), rotationnels (3 modes) et réactifs (2 modes) n'interviennent pas.

Récemment, une autre méthode de projection a été mise en oeuvre : la méthode de projection curviligne. Elle nécessite la transformation du gradient et de la hessienne en coordonnées curvilignes. Ce formalisme a l'avantage d'éliminer les phénomènes non physiques liés aux fréquences imaginaires obtenues pour les points autres que ceux de la MEP [180]. Cependant, cette méthode n'a pas été utilisée (puisqu'elle nécessite un développement mathématique) avancé dans l'étude de la réaction d'abstraction d'un atome d'hydrogène à partir du tétraméthylsilane (voir chapitre 6).

## Méthode d'ajustement

Dans notre approche, la forme analytique de la surface de potentiel effectif est obtenue à partir d'un ajustement à une double fonctionnelle de Morse dépendant de 29 paramètres
(C{1,2..29}) et ayant pour expression [162] :

$$
V(\rho, \delta) = \alpha[(1 - \exp(\mathcal{A} \times \delta - \mathcal{B}))^2 - C14] \\
+ \beta[(1 - \exp(\mathcal{D} \times \delta - \mathcal{E}))^2 - C28] + C29
$$

$\alpha$ et $\beta$ sont les fonctions amplitudes de la variable $\rho$ :

$$
\alpha = (C1 + (C2 \times \rho^{C3}) \times \exp(-C4 \times \rho))
$$

$$
\beta = (C15 + (C16 \times \rho^{C17}) \times \exp(-C18 \times \rho))
$$

$\mathcal{A}, \mathcal{B}, \mathcal{D}$ et $\mathcal{E}$ sont des polynômes du second degré de la forme :

$$
\mathcal{A} = (C5 + C6 \times \rho + C7 \times \rho^2)
$$

$$
\mathcal{D} = -(C19 + C20 \times \rho + C21 \times \rho^2)
$$

$$
\mathcal{B} = (C8 + C9 \times \rho + C10 \times \rho^2) - \log(C11 + C12 \times \rho + C13 \times \rho^2)
$$

$$
\mathcal{E} = -(C22 + C23 \times \rho + C24 \times \rho^2) - \log(C25 + C26 \times \rho + C27 \times \rho^2)
$$

Le but est d'obtenir par méthode itérative un ensemble de paramètres qui minimisent
la fonctionnelle non linéaire :

$$
\sum_{i=1}^{m} [V_{\eta}(\rho_{i}, \delta_{i}) - V_{\eta}^{ab \ initio}(\rho_{i}, \delta_{i})]^{2} \qquad (3.55)
$$

où *m* est le nombre de points de la grille tel que *m* ≥ 29 (nombre de paramètres {C} à optimiser) et $V_{\eta}^{ab \ initio}$ est le potentiel *ab* *initio* calculé pour chaque point de la grille en utilisant une méthode précise (dans notre cas c'est la CCSD(T)).

Le résultat de l’algorithme d’ajustement dépend de l’ensemble des paramètres initiaux.
De plus, l’ajustement doit être précis dans les régions physiquement importantes de la
surface de potentiel, i.e. la région proche de l’état de transition ainsi que les fonds des
puits.

**3.4.4 Equations couplées**

Dans le modèle 2D que nous avons utilisé, la fonction d'onde du système est obtenue en résolvant l'équation de Schrödinger suivante :

$$
\hat{H}\Psi = E\Psi \tag{3.56}
$$

$E$ représente l'énergie totale et $\hat{H}$ est l'hamiltonien donné par l'équation (3.20) où $V(\rho, \delta)$ est le potentiel effectif (équation (3.36)). Pour résoudre l'équation (3.56), l'approche suivie consiste à diviser l'espace de coordonnées en secteurs $i$ concentriques de largeur constante. La résolution est ainsi faite entre $\rho_{min}$ et $\rho_{max}$ pour $N_{\gamma}$ secteurs, $N_{\gamma}$ étant 90