List of Tables

x LIST OF TABLES 5.2 Ground state MgHe interaction energies calculated at the CCSD(T) level of theory in combination of the aug-cc-pV5Z basis set for He and the indicated basis set for Mg, except in cases 1,2 where the basis indicated is also used for He. In all cases, only the Mg 1s electrons have been kept frozen in calculations, exception made for case 1 where only the Mg 3s electrons have been correlated. The complete set of bond functions (33211) has been used except in case a where no bond functions are used and case b where only the 332 bond functions have been used. All energies are in 10 -6 a.u. and are counter-poise corrected. . . . . . . . . . . . . . . . . . . . . . . . . . .

Equilibrium distances

, r e , and the corresponding potential energy values, ǫ, using different sets of bond functions. All the the electrons of Mg have been correlated in the calculations, except the 1s electrons. E 0 is the ground state energy and r 0 the ground state distance expectation value. . . . . . .

5.4

Ab initio counter-poise corrected interaction energies : calculations have been carried out with an aug-cc-pV5Z (aug-cc-pVQZ) for He (Mg) basis set in combination of the complete set of bond function (33211). All energies are in ×10 -6 a.u. ∆V (r) is the relative error between V CCSDT f it The last two columns give the binding energies for a given N. All energies and standard deviations are expressed in cm -1 . The helium pair potential of Jeziorska [35] was used in combination with the CCSDT potential for Mg-helium interaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.3 DMC mean energies for pure and corresponding Mg doped helium clusters.

The last two columns give the binding energies for a given n. All energies and standard deviations are expressed in cm -1 . The data are obtained using the helium pair potential of Ref. [34] in combination with the CCSDT potential for the MgHe pair interaction. . . . . . . . . . . . . . . . . . . . . 9.4 DMC mean energies for pure and corresponding Ca doped helium clusters.

The last two columns give the binding energies for a given n. All energies and standard deviations are expressed in cm -1 and the time step, δt in atomic units. The helium pair potential of Ref. [34] have been used in the DMC calculations. The Ca-helium pair interaction was described by the CaHe CCSD(T) potential. . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.5 Comparison between the DMC binding energies (E n ) and the extrapolated energies (E f it ) from the fitting process according to Eq. 9.17. δE n and ∆E n are respectively DMC error bars and absolute errors of the fitting function with respect to the DMC points. Energy units are cm -1 . . . . . . . . . . .

A. [START_REF] Gspann | Physics of Electronic and Atomic Collisions[END_REF] Hartree-Fock energies for the MgHe complex and its components. The binding energies -E Bind. -are counter-poise corrected with respect to the (33211) bond functions. The aug-cc-pV5Z and aug-cc-pCVQZ basis sets were respectively used for He and Mg in combination with the bond functions (33211). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 CCSD(T) energies for the MgHe complex and its components. The binding energies are counter-poise corrected with respect to the (33211) bond functions. The Mg basis was aug-cc-pCVQZ+BFs and the He one was aug-cc-pV5Z+BFs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 CCSDT energies for the MgHe complex and its components. The binding energies are counter-poise corrected with respect to the bond functions (BFs=33211). The Mg basis was aug-cc-pCVQZ+BFs and that of He was aug-cc-pV5Z+BFs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 A. 4 Hartree-Fock energies for the MgHe complex and its components. Atomic basis sets are aug-cc-pV5Z+BFs and aug-cc-pVQZ+BFs for respectively He and Mg. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 A.5 CCSD(T) energies for the MgHe complex and its components. The Mg basis was aug-cc-pCVQZ and the He one was aug-cc-pV5Z. The binding energies are counter-poise corrected. . . . . . . . . . . . . . . . . . . . . . . 177 A.6 CCSD(T) energies for the MgHe complex and its components. The binding energies are counter-poise corrected. The Mg basis was aug-cc-pCVQZ and the He one was aug-cc-pV5Z. The complete set of bond functions (BFs) has also been used in combination of the atomic basis sets. . . . . . . . . . 178 A.7 MP2 energies for the MgHe complex and its components. The binding energies are counter-poise corrected with respect to the bond functions (BFs=33211). The Mg basis was aug-cc-pCVQZ+BFs and the He one was aug-cc-pV5Z+BFs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 A.8 MP4 energies for the MgHe complex and its components. Binding energies are counter-poise corrected with respect to the bond functions (33211). The Mg basis was aug-cc-pCVQZ+BFs and that of He was aug-cc-pV5Z+BFs. 180 A.9 CCSD(T) energies for the MgHe complex. Cp means that energies are counter-poise corrected. All energies are expressed in micro-hartree. The Mg basis was aug-cc-pCVQZ and the He one was aug-cc-pV5Z (see Fig. 5.5).181

List of Figures

4.1 Potential energy curves of the CaHe 1 Σ + state obtained at the CCSD(T) level of theory : circles picture the curves obtained by the use of cc-pVDZ and aug-cc-pVDZ respectively for calcium and helium. The second curves represents the use of cc-pVTZ for calcium and aug-cc-pVTZ for helium at the same level of theory. For both curves, the complete set of mid-bond functions (33211) has been used. . . . . . . . . . . . . . . . . . . . . . . . .

4.2

Influence of the bond functions for CaHe X 1 Σ + state at the CCSD(T) level. Basis sets are aug-cc-pV5Z for He and aug-cc-pCVQZ for Ca : curve with circles are obtained without bond functions, whole squares take into account the addition of the 322 set of mid-bond functions and diamonds result from the use of the complete set of mid-bond functions (33211) . . . 4.3 X 1 Σ + CaHe potential energy curves derived from different ab initio approaches in combination with aug-cc-pV5Z for He and aug-cc-pCVQZ for Ca plus bond functions (33211): curve with triangles are the CCSD(T) results, squares are for MP2, stars for MP4. Dots are results from Ref. [68] (calculated at the CCSD(T) level with the wtMCP pseudo-potential for Ca) and smooth curve represents the unpublished data from W. Meyer. . .

5.1

The core-valence correlation influence studied at the CCSD(T) level of theory in combination with the following basis sets : VQZ for both He and Mg, and aug-cc-pCVQZ ( aug-cc-pV5Z) for Mg (He). All curves are calculated taking into account the complete set of bond functions (33211). Notice that the CCSD(T) level with aug-cc-pCVQZ ( aug-cc-pCV5Z) for Mg (He) + BF is abbreviated as b1. . . . . . . . . . . . . . . . . . . . . . .

5.2

Influence of bond functions considering the CCSD(T),b1 method. . . . . . c) solid lines draw the results of the cor1 method (Eq. 5.6), dots those from the cor2 method (Eq. 5.7) and dashed lines those from the cor3 method (Eq. 5.8). 55 5.7 CBS extrapolation for the correlation energy of the He (a), Mg (b), MgHe (c) using for each system the three fitting methods (cor1, cor2 and Eq. 5.5). Ab initio points are performed at the method II level of theory. Fig. (d) shows the fitting HF energies : Curve with dots does not take into account the double ζ energies. In Fig. (a,b andc) solid lines draw the results of the cor1 method (Eq. 5.6), dots those from the cor2 method (Eq. 5.7) and dashed lines those from (Eq. 5.8). . . . . . . . . . . . . . . . . . . . . . . . 57 5.8 CBS extrapolation for the correlation energy of the He (a), Mg (b), MgHe (c) using for each system the two fitting methods (cor1, cor2). Ab initio points are performed at the method III level of theory. In all figures, curves with dots and solid lines do not take into account the double ζ energies. In Fig. (a,b andc) solid and dashed lines draw the results of the cor1 method (Eq. 5.6). The other curves are from the cor2 method (Eq. 5.7). . . . . . . 58 7.1 An expanded view of the OCS IR spectrum in pure 4 He droplets with N4 He = 6000 atoms (A) and in pure 3 He droplets with N3 He = 12000 atoms (B). The depletion is plotted versus the change in wave number with respect to the origin of the spectrum in (A) at ν 0 = 61.64 cm -1 . From Ref. [16]. . 86 7.2 A series of OCS IR spectra in 3 He similar to Fig. 7.1 (B) but with increasing average numbers N4 He of added 4 He atoms [ N4 He =0 (A), 7 (B), 25 (C), 35 (D), 60 (E), and 100 (F)]. The change in wave number in the abscissa is with respect to ν 0 = 2061.71 cm -1 . From Ref. [16]. . . . . . . . . . . . . .

8.13

Radial helium densities (ρ) of the He 12 cluster : curves with triangles and diamonds show the densities respectively obtained with a non optimized wave function and an optimized one. Horizontal solid line is for the bulk helium density. ρ is in Å -3 and R com in Å . . . . . . . . . . . . . . . . . . . 9.1 Comparison of the 3 1 P 0 1 ← 3 1 S 0 transition of Mg atoms picked up by He nanodroplets (thin line) and solvated in bulk liquid helium (thick line). This figure is reported in Ref. [28]. . . . . . . . . . . . . . . . . . . . . . . 112 9.2 Binding energies of the Mg@He N complexes when using helium pair potential of [35]. Solid line is the fitting function from Eq. 9.17. . . . . . . . . . 116 9.3 Absolute DMC energies obtained by using the Jeziorska pair potential for the helium atoms [35]. Circles are the energies of pure helium clusters (He n ) and triangles are that of the MgHe n clusters. . . . . . . . . . . . . . 117 9.4 Absolute DMC energies of the MgHe n clusters obtained by using the Aziz pair potential for the helium atoms [34]. . . . . . . . . . . . . . . . . . . . 118 9.5 Binding energies of the Mg@He N complexes when using helium pair potential of [127]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 9.6 Absolute DMC energies (a) of the CaHe N clusters (crosses) obtained by using the Aziz pair potential for the helium atoms [34] and our CaHe CCSD(T) potential. Circles depict the helium energies. of the Mg atom in the Mg@He n complexes according to Ref. [31]. . . . . . 127 9.8 Radial helium density (ρ in Å -3 ) in the MgHe 220 droplet (a) and density profiles for Mg@ 4 He N drops for N = 300, 500, 1000, 2000, 3000, and 5000 (from [32])(b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 9.9 From [30]: ion yield curves for Ca, Sr, Ba, and Xe atoms picked up by a helium nanodroplet beam (on the left). The shape difference reflects the fact that the metal atoms are located at the droplet surface and the xenon atom is solvated inside the droplet. The same conclusion is true in the case of the Mg atom (the right of the figure). . . . . . . . . . . . . . . . . . . . 130 9.10 Radial probability densities of Ca in helium clusters whose size are respectively, from the left to the right of the figure, 15, 32, 53, 83 and 133. These distributions are obtained when using helium pair potential from Ref. [34].

Résumé

Ces dernières années, les agrégats d'hélium superfluides ont fait l'objet de nombreuses études aussi bien expérimentales que théoriques. Le fruit de ces études a permis le développement de méthodes spectroscopiques innovantes (HENDI) utilisant les nanogouttes d'hélium comme l'ultime matrice, exploitant ainsi la très faible température de ce milieu particulier et sa faible intéraction avec les dopants pour une meilleure résolution spectrale. Cependant, un nombre important de questions subsistent quant aux agrégats d'hélium dopés, particulièrement, ceux dopés par les alcalino-terreux. En effet, la position d'une impureté au sein de la gouttelette d'hélium est loin d'être un problème trivial pour certaines espèces telles les alcalino-terreux. Ceci est particulièrement vrai dans le cas où l'impureté est l'atome de magnésium. Des preuves expérimentales d'un état solvaté du magnésium sont annoncées dans la littérature tandis que de récentes expériences laissent penser à une position plutôt surfacique du magnésium dans les agrégats d'hélium. Du point de vue théorique, la même ambiguïté persiste quant à la position de Mg dans la nanogouttes d'helium. Dans le but de contribuer a une meilleure compréhension des clusters d'hélium dopés par les métaux alcalino-terreux (Mg et Ca), nous avons, au cours de ce travail, du déterminer avec précision les energies d'intéractions des états fondamentaux des systèmes van der Waals CaHe et MgHe. Pour ce faire, des méthodes ab initio telles les approaches des clusters couplés (CC) mais aussi perturbationnelles (MP2 et MP4) ont été appliquées à ces deux systmes avec succès. Les meilleurs potentiels d'intéraction ont été utilisés par la suite comme potentiels d'intéractions de paire dans l'approche Monte Carlo à diffusion (DMC) en combinaison de deux types de potentiel d'intéraction pour l'hélium. Aussi bien pour CaHe n que pour MgHe n , des simulations DMC ont été produites depuis n = 1 jusqu' à n = 220, le résultat principal en est une position surfacique de l'impureté quelque soit l'alcalino-terreux considéré. Dans le cas particulièrement délicat des clusters d'hélium dopés par le magnésium, des calculs de DMC avec des contraintes géométriques montrent que le potentiel radial effectif de Mg dans He 20 et He 50 est plutôt plat. En fin, sont présentés également les résultats concernant la recombinaison dynamique de deux atomes de magnésium à l'intérieur d'un agregat d'environ 2000 atomes d'hélium utilisant une méthode basée sur un potentiel effectif pour l'intéraction He-He.

Abstract

During the last decades, superfluid helium clusters have been widely studied both experimentally and theoretically. As a result of the latter studies, a new spectroscopic domain has emerged (HENDI) where helium nanodroplets are used as ultimate matrices for accurate spectroscopic measurements, taking advantage of their very low temperature and their weak interaction with the impurity. However, many questions still are remaining about the helium nanodroplets, especially those doped with alkaline-earth atoms. In fact, the simple position of an impurity in this medium is far from being trivial for some doping species like the alkaline-earth atoms. This is particularly true when the impurity is the magnesium atom. Experimental evidence of a completely Mg solvated state is announced in the literature whereas very recent experiments advance the opposite situation for the Mg atom (near the surface). From the theoretical point of view, the position of the Mg atom in the helium droplet still remain ambiguous in the actual literature. In order to contribute to a better understanding of the alkaline-earth (Mg and Ca) doped helium clusters, we have determined, in this work, accurate interaction energies for both CaHe and MgHe van der Waals systems. For this aim, ab initio methods such as the coupled clusters (CC) as well as Møller-Plesset approaches (MP2 and MP4) have been successfully applied to both systems. The best interaction potentials have been then used as pair interactions for the diffusion Monte Carlo (DMC) approach in combination with two accurate helium pair interactions. For both CaHe n and MgHe n , DMC calculations have been carried out for n = 1 up to 220, the result was a surface location of the dopant whatever the latter is. In the particularly delicate case of Mg doped helium clusters, constrained DMC calculations have been performed for He 20 and He 50 . The results were a very flat energy profile in both cases. Finally, results concerning the dynamics of recombination of two Mg atoms based on an effective potential for helium inside an almost 2000 helium atom cluster are given.

Since the first experiments on doped 4 He droplets [START_REF] Gspann | Physics of Electronic and Atomic Collisions[END_REF][START_REF] Goyal | [END_REF][3][4][5][6][7][8], this ultra-cold environment has become an attracting and intriguing domain of physical and chemical research for all the scientific community : experimentalists [9] as well as theoreticians [10]. Why these exciting features of helium droplets ? We would simply say because this unusual medium presents unique properties. In fact the helium nanodroplets still are liquid at a temperature [11][12][13][14][15] of about 0.4 K. Furthermore, these droplets are superfluid and therefore present absence of internal friction which allows a quasi-free rotation of the molecular dopants [12,[16][17][18]. On the other hand, the extremely weak interaction energies with a doping impurity are at the same time a great advantage for spectroscopic experiments and a real challenge for theoreticians who aim at an accurate determination of those interactions. Among many other interesting properties of helium nanodroplets, we can cite the ability of the latter to rapidly dissipate the excess energy of excited impurities [19]. All these exceptional characteristics of helium droplets make them ideal matrices in order to carry out highly accurate spectroscopic experiments. The latter constitute a recently emerged domain: the HElium Nanodroplet Isolation Spectroscoy (known under the acronym HENDI). Finally, helium nanodroplets have given to scientists the unique opportunity for studying chemical reactions and very reactive species in this medium [20][21][22][23].

However, the highly quantum nature of helium nanodroplets in addition to its interaction weakness with dopants make the "simple" solvation process of an atomic dopant not really understood nor rationalized up to now.

Considering the alkaline atom cases, despite their extremely weak van der Waals interactions with helium, they support at least one bound state and theoretical calculations as well as experiments have demonstrated a surface location for all alkaline atoms [24][25][26][27].

On the other hand the location of the neutral alkaline-earth atoms is more ambiguous especially for the magnesium atom. Some experiments advance a completely solvated state of this impurity. In fact Reho et al. [28] measured the laser-induced fluorescence spectra of the 3 1 P 0 1 ← 3 1 S 0 transition of Mg solvated in helium nanodroplets and compared their result to the same kind of spectra of Mg in bulk liquid helium measured by Moriwaki et al. [29]. The relatively good overlap between the two spectra led Reho et al. to suppose a completely solvated state of Mg. Other experiments imply a surface location of the Mg atom. Ren and co-workers [30] have recently measured the electron energy dependence of the ionization yield of alkaline-earth and Xenon atoms picked-up by helium clusters and found a qualitative shape difference between the yield curves of species completely solvated in the droplet (Xe) and species located in the surface region (Ca, Sr and Ba). The measurements also suggest a surface location of Mg.

From the theoretical point of view the situation is not clearer. In fact, Mella et al. [31] have performed diffusion quantum Monte Carlo (DMC) calculations of the Mg@He n systems with n up to 50. The results of these calculations led to ambiguous mixed cluster structures since the position of the Mg impurity in the helium nanodroplets seems to depend on the cluster size. Secondly, the results published by Mella et al. appear to fail to reproduce a bound state for very small Mg doped clusters.

On the other hand, an alternative theoretical approach has been applied to the Mg doped helium clusters by Hernando et al. [32,33]. Using a density functional theory (DFT) treatment, the authors found a completely solvated state of Mg. However, we are not very confident in these results since the DFT is a semi-empirical method the results of which directly depend on the parametric functional used to describe the correlation part of the energy.

In order to bring some contributions to the understanding of the Mg doped helium nanodroplets, we have performed accurate diffusion Monte Carlo calculations exactly (within statistical errors) solving the many-body Schrödinger equation for MgHe n system, with n up to 220, using two types of helium pair potentials [34,35] in combination with a very accurate MgHe ground state pair potential which we calculated. The less ambiguous CaHe n has also been treated by the DMC method and some results are also reported.

The weakness of the van der Waals interactions, which are briefly introduced in Chapter 2, involving helium and a dopant implies the choice of appropriate and accurate ab initio methods to determine the dopant-helium interactions. For this aim, ab initio methods are briefly introduced in Chapter 3. The most accurate among the latter have been used in order to calculate the ground state potential energy interactions of CaHe and MgHe respectively in Chapter 4 and Chapter 5. While Chapter 6 is dedicated to the presentation of the variational quantum Monte Carlo (VMC) and diffusion quantum Monte Carlo (DMC) methods, Chapter 7 exposes some interesting characteristics of helium droplets. In Chapter 8, the computational details allowing accurate and robust DMC calculations are given. DMC results for both Ca@He n and Mg@He n are presented in Chapter 9. Finally, a classical dynamics simulation based on an effective potential for helium, which takes into account the large zero point energy and well estimates the helium density, is presented in Chapter 10.

Chapter 2

Weak interactions 2.1 Introduction

The weak interaction energies between atomic or molecular system have been named van der Waals interactions in recognition of Johannes van der Waals work about the equation of state of ideal gas :

P V = RT (2.1)
which he modified into the following equation of state to better take into account the imperfection aspect of real gases:

(P + a V 2 )( V -b) = RT (2.2)
In Eqs. 2.1 and 2.2 R is the gas constant, T , the temperature and V , the molar volume of the fluid. However, in equation 2.2 two new parameters appear in the equation of state : a and b. The first reflects the strength of molecular or atomic attraction and the second (b) is proportional to the volume of the particle. The work of van der Waals has encouraged scientists of the 20 th century (Debye, Keesom, London,...) to study the interactions between molecules and atoms.

Dispersion interactions

After the main investigations of London in this domain [36][37][38], four primary types of intermolecular interaction energies have been established : electrostatic, induction, dispersion and exchange. These interactions belong to two main categories : long range interactions and short range interactions. Long range effects occur by definition at large Chapter 2. Weak interactions interatomic (intermolecular) distances, where the wave functions of each components do not overlap whereas the short range interactions are due to the overlap of the wave functions.

Here, we are confronted with dispersion forces since all the pair interactions we need in this study (He 2 , CaHe and MgHe) only imply atomic species which obviously do not have permanent multipoles. Moreover, the closed shell nature of the systems excludes chemical bonding.

London provided the first expression of the long range attractive interactions (known as London forces) based on quantum mechanics [START_REF] Rigby | The Forces between Molecules[END_REF] :

U dis (r) = -3α 2 E I 4(4πǫ 0 ) 2 r 6 = -C 6 r 6 (2.3)
where α and E I are respectively the polarizability and the ionization energy of the system and ǫ 0 the permittivity of vacuum. On the other hand the short range interaction, due to the overlap of the electronic cloud of atoms, can be generally expressed as

U rep (r) = Ae -βr (2.4)
where A and β are real parameters depending on the nature of the system and r is the interatomic (or intermolecular) distance.

Combination of the long range interaction, ie Eq. 2.3, and the short range repulsion, ie Eq. 2.4, leads to the basic description of the system by means of the following expression:

U tot = U rep + U dis (2.5) = Ae -βr - C 6 r 6 (2.6) (2.7)
Since the dispersion energy contains higher order terms arising from the contributions of other instantaneous multipoles such as quadrupoles and octopoles, the total function with which we model (fit) the pair interaction energies of our systems is written as :

V pair (r) = Ae -βr - C 6 r 6 - C 8 r 8 - C 10 r 10 (2.8)
where C i are i th order dispersion coefficients.

Chapter 3 Electronic Structure Methodology

In this chapter we will describe the main electronic structure theories in order to understand how the time-independent Schrödinger equation is solved. Most of these methods have been used to determine the electronic energy curves of the MgHe and CaHe van der Waals systems. All the following methods, which allow to solve the nonrelativistic timeindependent Schrödinger equation, are based on the Born-Oppenheimer approximation.

It assumes a separation between the motion of the electrons and the nuclei, leading to two eigenvalue equations: the electronic and the nuclear Schrödinger equations. In this chapter we will show how we can extract the two main unknowns, the energy and the wave function, from those equations using different approaches.

Time-independent Schrödinger Equation

All ab initio methods which are presented in this chapter to perform quantum chemical calculations, aim at solving the stationary nonrelativistic Schrödinger equation [START_REF] Szabo | Modern Quantum Chemistry[END_REF]. This (time-independent) equation is given as

Ĥ|Ψ = E|Ψ (3.1)
In the previous equation the total energy, E, is calculated by applying the time-independent molecular Hamiltonian, Ĥ, on the wave function, |Ψ , describing the system's properties.

The molecular Hamiltonian can be described in terms of distinct parts as Ĥ = Te + TN + VNe + VNN + Vee (3.2)

In the expression above one can notice that the spin contributions due to the electronic and nuclear particles of the system are neglected. The different parts of the right hand side of equation ( 3.2) are explicitly expressed as :

Te = - h2 2m e i △ i (3.3) TN = - h2 2 A △ A M A (3.4) V N e = - e 2 4πǫ 0 A i Z A r iA (3.5) V N N = e 2 4πǫ 0 A B>A Z A Z B R AB (3.6) V ee = e 2 4πǫ 0 i j>i 1 r ij (3.7)
where h is the reduced Planck's constant, △ is the Laplacian operator, m e and M A are the mass of the electron and the nucleus, respectively. These masses are used for the definition of the two kinetic energy operators, Te and TN , respectively for the electrons and the nuclei. The term V N e represents the Coulomb potential between an electron i and a nucleus A. Finally, the last two terms, V ee and V N N , represent the overall Coulomb repulsion between pairs of electrons and nuclei, respectively. Z A describes the atomic number of nucleus A, r ij the distance between electron i and j, r iA the distance between electron i and nucleus A, R AB that between two nuclei A, B. Notice that e and ǫ 0 are respectively the elementary charge and the permittivity constant of vacuum. A purely potential operator V can be defined as

V = V N e + V N N + V ee (3.8)
In order to solve equation ( 3.1), some approximations are necessary, among them, the main approximation used is the Born-Oppenheimer approximation.

Born-Oppenheimer Approximation

Unfortunately, the Schrödinger equation, eq. 3.1, is not exactly solvable for most systems, except in the case of one-electron systems such as the hydrogen atom or the H + 2 molecule. One way to overcome this inability is a separation between the nuclear terms and the electronic terms [START_REF] Helgaker | Molecular Electronic-Structure Theory[END_REF]. The most common way to do so is applying the Born-Oppenheimer approximation (BOA) [START_REF] Born | [END_REF] to equation 3.1. Since the mass difference between the electron and the nucleus is very important, for example in the case of the helium atom, the mass ratio of the nucleus to the electron is more than 3600, electrons are assumed to instantaneously adapt to the nuclear motion. In other words, the geometry of the nuclei can be considered as fixed when studying the electrons. Applied to the time-independent Schrödinger equation, the BOA leads to two distinct equation : the electronic Schrödinger equation and the nuclear Schrödinger equation which are introduced in the following.

Electronic Schrödinger Equation

Under the Born-Oppenheimer approximation (BOA) a separation between the motion of the electrons and the nuclei is done. Ĥe , the electronic Hamiltonian in Eq. 3.2, is the part of the total Hamiltonian which accounts for the electrons. In this electronic operator, Ĥe , the kinetic energy term for the nuclei and the electrostatic nuclear interaction, are dropped :

Ĥe = Te + V N e + V ee (3.9) = - h2 2m e i △ i - 1 4πǫ 0 A i Z A e 2 r iA + i j>i e 2 r ij (3.10)
The purely electronic Schrödinger equation then becomes

Ĥe |Ψ e (r i , R A ) = E e (R A )|Ψ e (r i , R A ) (3.11)
In Eq. 3.11 |Ψ e (r i , R A ) is the electronic wave function and explicitly depends on the electron coordinates, r i , and parametrically on the nuclear coordinates, R A . The electronic energies, E e (R A ), also parametrically depend on the vector-coordinate R A .

In practice, the stationary electronic Schrödinger equation, Eq. 3.11, is solved and E e (R A ) determined. The potential energy operator for a given electronic state -that is to say a point in the potential energy surface -is determined as

Ve (R A ) = E e (R A ) + VNN (R A ) (3.12)
The stationary points evaluated at different nuclear coordinates, by operation of the Ve operator onto the wave function, will produce potential energy surfaces or potential energy curves in case of diatomic molecules.

Nuclear Schrödinger Equation

Another consequence of the Born-Oppenheimer approximation is the appearance of the following nuclear Schrödinger equation

ĤN |Ψ N (r i , R A ) = E N |Ψ N (r i , R A ). (3.13)
The average values of the electronic coordinates calculated over the electronic wave function can be applied to the motion of the nuclei. The nuclear Hamiltonian, ĤN , is written as the sum of the electronic operator V e (R A ) of Eq. 3.12 and the kinetic energy operator for the nuclei TN given in Eq. 3.2 :

ĤN = - h2 2M A A △ A + 1 4πǫ 0 A B>A Z A Z B R AB + E e (R A ) (3.14) = - h2 2M A A △ A + Ve (R A ) (3.15)
The nuclear energy, E N , in Eq. 3.13, now describes the total energy, E t , in the Born-Oppenheimer approximation and accounts for the electronic and vibrational energy of a given system. The BOA is suited for cases when electronic and nuclear Hamiltonians are well separated, in other words, when nuclear motion cannot cause a change in the electronic state. The potential energy curves described in this work concern the electronic ground states, which are well separated from the excited electronic states. This justifies the use of the BOA.

The Variational Principle

Given an arbitrary wavefunction for a system, the corresponding total energy is, by definition, the expectation value of the Hamiltonian operator. This concept is translated into the following equation :

E = Ψ| Ĥ|Ψ Ψ|Ψ (3.16)
The wave function can be expanded in terms of eigenfunctions of the Hamiltonian operator, Ĥ, as :

|Ψ = k c k |Φ k (3.17)
where the coefficients in the expansion are complex numbers defined as :

c k = Φ k |Ψ .
Introducing Eq. 3.17 in Eq. 3.16 leads to the following equations :

E = k,l c * k c l Φ k | Ĥ|Φ l k,l c * k c l Φ k |Φ l (3.18) = k,l c * k c l E l Φ k |Φ l k,l c * k c l Φ k |Φ l (3.19) = l |c l | 2 E l l |c l | 2 (3.20)
By definition of Eq. 3.17, the ground state energy is E 0 . It is therefore clear that :

E 0 ≤ E (3.21)
The variational theorem gives an approximation to the ground state wave function : given a parameterized wave function, its parameters are adjusted to obtain the lowest energy. The optimized function is then an approximation to the ground state wave function.

Solving the Electronic Schrödinger Equation

The solutions of the electronic Schrödinger equation (Eq. 3.11) for different geometries provide potential energy surfaces or curves. These curves will be used as pair potentials in DMC and dynamic calculations. In order to solve Eq. 3.11 it is important to appropriately represent the electronic wave functions. A simple and approximate way to represent a wave function is the use of Slater determinants.

Slater Determinant

According to the Pauli exclusion principle, a wave function must be antisymmetric under exchange of two fermions. This principle is satisfied by using a Slater determinant to antisymmetrize a wave function which can then describe the ground state of a molecule. The Slater determinant in its abbreviated form,

ψ = |ϕ 1 ϕ 2 • • • ϕ n | (3.22)
is written, in a less compact shape, as

ψ = 1 √ n! ϕ 1 (x 1 ) ϕ 2 (x 1 ) • • • ϕ n (x 1 ) ϕ 1 (x 2 ) ϕ 2 (x 2 ) • • • ϕ n (x 2 ) . . . . . . . . . . . . ϕ 1 (x n ) ϕ 2 (x n ) • • • ϕ n (x n ) (3.23)
is described by mono-particle functions, ϕ(x), also called spin orbitals. A spin orbital is composed as a product of a spatial orbital, φ(r), and a spin function α(σ) or β(σ)

ϕ(x i ) = ϕ(r, σ) = φ(r) α(σ) β(σ) (3.24)
The spin-orbitals are assumed to be orthonormal that is to say, they verify the following relation

ϕ i |ϕ j = δ ij (3.25)
where δ ij is the Kronecker delta which equals 0 if i = j and 1 otherwise. The variational principle (see section 3.2) states that the best wave function Ψ 0 (x) is the one that gives the lowest energy

E 0 = Ψ 0 (x)| Ĥe |Ψ 0 (x) (3.26)
Now, after having applied the BOA to the Schrödinger equation, in a first step the energy minimization, within the model of an antisymmetrized product, according to the variational principle, leads to the Hartree-Fock equations.

Hartree-Fock Equations

The Hartree-Fock (HF) method [43][44][45] is an approximation which determines the ground state energy and wave function for an N-electron system. It approximates the exact wave function by a single antisymmetrized product (a determinant) that is optimized by solving the Hartree-Fock equations iteratively, in a process known as the self-consistent field (SCF) method.

Minimizing the energy in Eq. 3.26 is made by varying the space part of the spin-orbitals. The Hartree-Fock equation f (

x i )|ϕ(x i ) = ε i |ϕ(x i ) (3.27)
is an eigenvalue equation whose solution gives the minimum energy of the system, E 0 , and also the optimal spin-orbitals, ϕ(x i ). The single particle Fock operator f (x i ) = ĥ(x i ) + vHF (x i ) (3.28) with ĥ(

x i ) = - h2 2m e △ i + A Z A e 2 4πǫ 0 r iA vHF (x i ) = N j Ĵj (x i ) + Kj (x i )
is defined by the single electron Hamiltonian, written as ĥ(x i ), and the effective Hartree-Fock potential operator vHF (x i ). The latter is defined by two components : the first term is the classical Coulomb operator

Ĵj (x 1 )|ϕ i (x 1 ) = ϕ * j (x 2 ) e 2 4πǫ 0 r12 ϕ j (x 2 )dx 2 |ϕ i (x 1 ) . (3.29)
This operator represents the interaction of one electron, in this case electron i, with the remaining N -1 electrons. The second term is known as the exchange operator and is of purely quantum nature :

Kj (x 1 )|ϕ i (x 2 ) = ϕ * j (x 2 ) e 2 4πǫ 0 r12 ϕ i (x 2 )dx 2 |ϕ j (x 1 ) . (3.30)
The two-electron potential operator 1 r12 describes the electrostatic potential energy between electron 1 and electron 2. Application of K causes the position and spin of electron i to be exchanged with that of electron j.

Restricted Closed-Shell Hartree-Fock

When solving the Hartree-Fock (HF) equation, Eq. 3.27, one needs to evaluate the shape of the spin orbitals. For a closed-shell molecule, the pair of spin orbitals have the same spatial function :

ϕ 2i (x) = ψ i (r)α(σ) i = 1, 2, ...K ψ i (r)β(σ) (3.31)
Introducing this equation into the HF Eq. 3.27 results in two distinct HF equations,

f (x 1 )|ψ i (r 1 )α(σ 1 ) = ε i |ψ i (r 1 )α(σ 1 ) (3.32) f(x 1 )|ψ i (r 1 )β(σ 1 ) = ε i |ψ i (r 1 )β(σ 1 ) (3.33)
given in terms of their spin functions, α(σ i ) and β(σ i ). In order to get an expression exclusively in terms of the spatial orbitals, it is necessary to remove the spin function from the Fock operator. To do so, the spin orbital is replaced by its spatial orbital and spin functions. In a closed-shell system, the contributions from the α terms and those from β are equal. It is then sufficient to multiply Eq. 3.32 from the left by α * (σ i ) and integrate over its spin

α * (σ 1 ) f(x 1 )α(σ 1 )dσ 1 |ψ i (r 1 ) = ε k |ψ k (r 1 ) (3.34)
In Eq. 3.34, f is the closed-shell Fock operator and is written in the following form f (r 1 ) = ĥ(r 1 ) +

N/2 i=1 2 Ĵi (r 1 ) -Ki (r 1 ) (3.35)
where N/2 is the number of α or β electronic particles. Introducing the permutation operator, the Fock operator can be written in terms of spatial orbitals as f (r 1 ) = ĥ(r 1 ) +

N/2 i=1 ψ * (r 2 ) e 2 4πǫ 0 r12 (2 -P12 )ψ i (r 2 )dr 2 (3.36)
This involves that the Coulomb and exchange operators are now defined with respect to the spatial orbitals and are respectively written as

Ĵi (r 1 )|ψ k (r 1 ) = ψ * i (r 2 ) e 2 4πǫ 0 r12 ψ i (r 2 ) dr 2 |ψ k (r 1 ) (3.37) and Ki (r 1 )|ψ k (r 1 ) = ψ * j (r 2 ) e 2 4πǫ 0 r12 ψ k (r 2 ) dr 2 |ψ j (r 1 ) . (3.38)
Since the Fock operator has been defined in terms of the spatial orbitals, the Roothaan-Hall equations provide an algebraic method that allows to solve the Hartree-Fock equations. These equations are introduced in the following section.

Roothaan-Hall Equations

In the previous section it has been demonstrated that the spatial orbitals were obtained by integrating out the spin functions. Now, The Hartree-Fock Eq. 3.27 can be rewritten in terms of these spatial orbitals following the equation below

f (r 1 )|ψ k (r 1 ) = ε k |ψ k (r 1 ) (3.39)
Unfortunately, in case of molecular systems, Eq. 3.39 can not be solved analytically, it is hence converted into a set of algebraic equations by introducing a set of known functions, as introduced by Roothaan and Hall [46,47]. The advantage of the Roothaan-Hall equation formalism is that the molecular orbitals can be expanded into a linear combination of known one-electron functions :

|ψ k = K µ=1 C µk |φ µ (3.40)
In the previous expression, the basis functions are written with Greek indices and the molecular orbitals with Latin ones. Inserting Eq. 3.40 into Eq. 3.39, multiplying from the left with φ * ν and integrating gives

K µ C µk Fνµ φ * ν (r 1 ) f (r 1 )φ ν (r 1 )dr 1 = ε k K µ C µk Sνµ φ * ν (r 1 )φ ν (r 1 )dr 1 (3.41)
Eq. 3.41 can be written in a more compact shape as

FC = SCε (3.42)
where F is the Fock matrix, S the overlap matrix and the matrix C which contains the coefficients for the molecular orbitals φ i in the column i. The molecular orbital energies are given by the diagonalized matrix ε. All these matrices are hermitian which means that for real orbitals, they are real and symmetric.

The Fock matrix F depends on the coefficients, C. This makes the Roothaan-Hall equations nonlinear and in order to solve a nonlinear system an iterative approach is used : the self-consistent field procedure which is discussed below (Section 3.3.5).

Self-consistent Field

The Self-consistent Field procedure (SCF) is an iterative method which is used to solve the Hartree-Fock equations [43]. It has also been applied to the Roothaan-Hall equations. The SCF method first solves the nonlinear equations with a guessed set of molecular orbital coefficients, obtained for example by diagonalizing the one-electron part of the Hamiltonian. Since the Fock operator depends on these coefficients, after one iteration new coefficients are calculated with new energies. The procedure is repeated until selfconsistency has been reached and the Hartree-Fock energy is then given as

E HF = Ψ 0 (x)| Ĥe |Ψ 0 (x) (3.43)
The Hartree-Fock energies depend on the basis set quality and the Hartree-Fock limit is reached (saturation of the one electron basis) when increasing the basis set quality does not change the energy.

Configuration Interaction (CI)

The correlation energy, E corr , is defined with respect to the Hartree-Fock energy and the total non-relativistic energy (≡ E e from eq. 3.11) as

E corr = E T ot. -E HF (3.44)
This correlation energy is a negative quantity since the Hartree-Fock energy defines an upper bound to the total exact energy. To account for the correlation energy, many methods exist. One way to do so is to expand the exact electronic wave function into a linear combination of Slater determinants (CI). A single determinant ground-state wave function for the N electrons, following the HF scheme, is composed with a set of 2K spin orbitals

|Ψ ξ 0 = |ϕ 1 ϕ 2 ...ϕ i ϕ j ...ϕ N (3.45)
The configuration interaction method is based on the fact that the single Slater determinant is not sufficient to correctly describe the wave function. In fact, many other determinants can be formed from different possible excitations of the electrons. A single excited determinant with respect to the HF determinant is written as (motion of and electron from an occupied spin orbital to a virtual spin orbital)

|Ψ ξ i = |ϕ 1 ϕ 2 ...ϕ ξ ϕ j ...ϕ N (3.46)
This determinant describes the fact that a single electron is relocated from its occupied spin orbital i to one of the virtual spin orbitals ξ. By analogy a doubly excited determinant can be formed by relocating two electrons from their initial positions, i and j, to two virtual spin orbitals, ξ or κ

|Ψ ξκ ij = |ϕ 1 ϕ 2 ...ϕ ξ ϕ κ ...ϕ N . (3.47)
The excitations of higher order continue until all N electrons have been moved to various virtual spin orbitals. The full-CI wave function is a sum of the excited determinants and the HF determinant :

|Φ 0 = |Ψ 0 + i,ξ c ξ i |Ψ ξ i + i>j;ξ>κ c ξκ ij |Ψ ξκ ij + i>j>k;ξ>κ>η c ξκη ijk |Ψ ξκη ijk + ...
A linear combination of the initial configurations is done so that the configurations are eigenfunctions of the spin and angular momentum operators. The configurations obtained by this linear combination are termed CSF (configuration state functions). In practice, for large systems a full CI calculation is computationally infeasible. A common way to overcome this inability is to truncate the CI method to only singles and doubles excitations (CISD). But a problem appears with the truncation: the truncated method is no longer size consistent nor size extensive. In fact, a method is size extensive when the calculated energy for N non-interacting atoms equals the sum of N times the energy of the single atoms. Whereas the size consistency notion refers to the fact that at the dissociation of a molecule, the energy of the molecule should be equal to the sum of energies of each part.

Perturbational theories

Until now, we have dealt with variational methods which solve the many-body timeindependent Schrödinger equation. All these methods are governed by the variational theorem (Section 3.2). Another alternative approach to perform accurate quantum calculations is using perturbational methods, where the remaining correlation contribution to the wave function beyond a HF-SCF or MCSCF wave function is assumed to be small.

Møller-Plesset Perturbation Theory

The Møller-Plesset Perturbation Theory (MPPT) [48,49] uses the perturbation theory to calculate the coefficients of the CI expansion. Basically, the electronic correlation effects are assumed as a perturbation, V , in front of the all-electron Fock operator, F , (described in Section 3.3.2). This assumption leads to decomposition of the Hamiltonian as :

Ĥ = F + V (3.48)
The V operator in Eq. 3.48 is called the fluctuation operator. The starting point is the Hartree-Fock wave function (the unperturbed state). The Rayleigh-Schrödinger perturbation theory is applied and gives the perturbative corrections to the wavefunction Ψ (1) , Ψ (2) , ..., which are constructed from the single, double, etc. excitations as shown in the CI expansion of Eq. 3.4. According to the perturbation theory, the corrections to the energy to first order (E (1) ), second order (E (2) ), ... and the corresponding contributions to the coefficients (c ξ i , c ξκ ij , etc.) in the CI expansion are determined. Perturbation theory to first order in the energy just gives the Hartree-Fock energy.

Concerning the second order energy correction, E (2) , in the basis of the occupied (i, j, ...) and unoccupied (ξ, κ, ...) spin orbitals, it is evaluated to be [START_REF] Helgaker | Molecular Electronic-Structure Theory[END_REF] :

E (2) = 1 4 i,j;α,β ij||αβ 2 ε i + ε j -ε α -ε β (3.49)
In the previous equation ( 3.49), the quantity, E (2) , is termed the MP2 correlation energy and ε i , ε j , ... are the Hartree-Fock orbital energies. The quantity ij||αβ is defined as the difference between two-electron integrals over spin orbitals :

ij||αβ = ij|αβ -ij|βα . (3.50)
Møller-Plesset perturbation theory to second order (in the energy), MP2, is a computationally inexpensive correlation method commonly used, especially for large systems.

A higher order perturbation theory (in energy) exist : MP4. The latter needs the second order correction to the wavefunction, Ψ (2) . MP4 has, therefore, contributions from single, double, triple and quadruple excitations. However, accounting for triple excitations is more difficult (and time consuming) than accounting for the quadruples. For this reason the triple excitations are often neglected, giving MP4(SDQ) theory. The MP4 theory which explicitly accounts for triple excitations is denoted MP4(SDTQ).

When higher accuracy is required, it appears more interesting to use coupled cluster methods. The latter make the subject of the next section.

Coupled Cluster Theory

Electron correlation calculations can be performed by means of another approach which is the coupled cluster (CC) theory [50,51]. Via the cluster operator, T , a coupled cluster wavefunction is formulated in these terms :

Ψ = e T Ψ 0 (3.51)
where Ψ 0 is the HF determinant. The cluster operator , T , is constructed from one body, two body, ..., cluster terms, T1 , T2 , etc. These cluster terms represent single excitation, double excitation, etc. operators :

T = T1 + T2 + T3 + ... (3.52)
The right hand side terms of Eq. 3.52 are explicitly written as : The operators {â + a } are creation operators that is to say they generate an electron in spin orbital a whereas {â i } are annihilation operators which remove an electron from orbital i. Once applied together, they represent the excitation of an electron from orbital i to orbital a. The numerical coefficients, {t a i }, {t ab ij }, {t abc ijk }, etc., are known as the cluster amplitudes. Asymptotically expanding the e T operator leads to the following equations : In a similar manner, as in the CI case, the coupled cluster is in practice truncated after double excitations :

T1 = i,
e T = 1 + T + 1 2 T 2 + 1 3! T 3 + ... ( 3
Ψ = e T1 + T2 Ψ 0 (3.63) = {1 + T1 + T2 + 1 2 T 2 1 + T1 T2 + 1 2 T 2 1 T2 + 1 6 T 3 1 + 1 2 T 2 1 + 1 24 T 4 1 + ...}Ψ 0 (3.64)
Since the wave function is not linear in the cluster amplitudes, {t a i },{t ab ij }, etc., it cannot be calculated using standard eigenvalue methods, it is rather obtained iteratively by solving the Schrödinger equation in the space of the configurations used (reference state, single and double excitations). Iteratively solved are the following equations :

Ψ 0 | Ĥ|e T Ψ 0 = E (3.65) Φ a i | Ĥ|e T Ψ 0 = t a i E (3.66) Φ ab ij | Ĥ|e T Ψ 0 = t ab ij E (3.67)
In Eqs. 3.65, 3.66 and 3.67, the eigenvalue, E, represents the coupled cluster energy. Solving the latter three equations leads to the CCSD method (coupled cluster with singles and doubles excitations). To increase the accuracy of the CCSD method, the triple excitations are accounted for by a perturbative treatment, using the coupled cluster wave function as the reference state [52] :

∆E triples = i<j<k;a<b<c (1 + T1 + T2 + T3 )Ψ 0 | Ĥ|Φ abc ijk (1 + T2 )Ψ 0 | Ĥ|Φ abc ijk (ε i + ε j + ε k ) -(ε a + ε b + ε c ) (3.68)
where ε i , ε j , ... are the HF orbital energies.

The method obtained (CCSD with perturbative triples) is called CCSD(T). It has been used in this work for generating accurate energies for CaHe and MgHe molecules.

Truncating the cluster expansion after the T3 term, becomes an interesting alternative since modern computers have large memory and high frequency CPU. When explicitly including the connected components of the triple excitations, the method is denoted as the CCSDT method. Although the computational cost is nowadays too high to allow this method to be commonly used, we have managed to obtain a complete potential energy curve for the MgHe ground state complex.

As the coupled cluster theory is based on a single reference determinant (HF determinant), the accuracy of the results strongly depend on the fact that the reference determinant is dominant in the coupled cluster expansion for all the molecular geometries.

Basis Sets

In this work the basis sets used are the augmented correlation consistent (cc) basis sets, aug-cc-pVXZ [53,54]. They have been constructed so that, as the cardinal number of the basis set, X, increases, the description of electron correlation becomes systematically better and predictable. The minimal basis set is cc-pVDZ. The latter is improved by adding atomic functions which are chosen so as to maximise their contribution to the electron correlation. In this way a systematic improvement in the description of the correlation energy is observed when the cardinal, X, increases. The main advantage of such systematic improvements of electron correlation description is that the energies from a sequence of correlation consistent calculations can be fitted to smooth monotonic functions. It is a way to extrapolate the finite-basis energies to the complete basis set limit (see Section 3.9).

The aug-cc-pCVXZ basis sets [55] have also been used for the alkaline earth atoms (Ca and Mg); these basis sets are based on the corresponding aug-cc-pVXZ. However, the aug-cc-pCVXZ basis sets have additional tight functions added which allow them a better description of the correlation of core electrons and between core and valence electrons. These additional correlations are of extreme importance especially in the case of alkaline-earth bound with the helium atom since the interaction is very weak.

Basis Set Superposition Error

Using finite basis sets in calculations of potential energy surfaces involves the presence of basis set superposition errors. This phenomenon is due to the fact that for a given bound molecule AB, the atom A can be stabilised by the close presence of the basis functions of atom B and vice versa. The system is hence not only bound by any true interaction between A and B but also by this additional superposition effect. A possible and approximate correction of this effect is obtained via the counterpoise method of Boys and Bernardi [56]. This method involves the calculation of the energy of each atom or fragment both with its basis functions, E A , E B , and with the basis functions of the entire system E A(B) , E (A)B . This counterpoise correction is respectively given for A and B by:

∆E CP A = E A(B) -E A (3.69) ∆E CP B = E (A)B -E B (3.70)
The total counter poise correction to the interaction energy is the sum of the counterpoise corrections, ∆E CP A + ∆E CP B . The counter poise corrected interaction energy is finally written as:

E corrected AB = E A + E B -E AB + ∆E CP A + ∆E CP B (3.71)
Another way to accurately determine interaction energies can be achieved by the complete basis set extrapolation limit which is the subject of the next section.

Complete Basis Set Extrapolation

In order to extrapolate energies to the complete basis set, complete basis set extrapolation methods represent, in principle, the highest level of theory available for this aim. These methods employ an accurate ab initio method in combination with the correlation consistent basis sets. This method gives reliable results for strongly bound systems. It seems, however, to not be satisfying for weakly bound systems (refer to Chapter 5).

Binding energy definition

In the case of a diatomic molecule, the binding energy is defined as the difference between the molecular electronic energy and those of its components. This definition, of course, can be generalised to any molecule which contains more than two atomic components. For the AB molecule, we can formulate1 its binding energy as

E bind = E AB -(E A + E B ) (3.72)
In the previous relation, E bind is the binding potential energy between atom A and B, E A and E B their respective energies. Each of these energies can be decomposed into two contributions: a Hartree-Fock part and a correlation one.

Hartree-Fock energy

The Hartree-Fock part of the energy, E HF , is expected to evolve, when increasing the number of ζ, x as

E HF (x) = E ∞ HF + αe -βx (3.73)
where E ∞ HF is the Hartree-Fock energy for an infinite basis, α and β are the fitting parameters. The Dunning-Feller [53,57,58] (exponential form (Eq. 3.73)) has extensively been shown to better extrapolate the HF part of the binding energy [59] than any power law does.

The correlation energy

A large number of extrapolation schemes have been proposed in the literature for this purpose over the last years. One of the most commonly used is the " 1 3 " scheme of Helgaker et al. [60]:

E corr (x) = A 1 + B 1 x 3 (3.74)
Another extrapolation scheme not so far from the previous one has been proposed by Gdanitz [61] and is expressed in the following way

E corr (x) = A 2 + B 2 (x -1 4 ) 3 (3.75)
Two possible other extrapolation schemes have been used in this work. The first one is (Eq. 3.76) similar to the extrapolation scheme of Eq. 3.74, except that the exponent is now an adjustable parameter (C 3 ). The other scheme assumes both Hartee-Fock energy and correlation energy can be extrapolated using the same functional form (Eq. 3.77).

E corr (x) = A 3 + B 3 x C 3 (3.76) E corr (x) = A 4 + B 4 e -C 4 x (3.77)
In Eqs. 3.73 to 3.77, x is the cardinal number of the basis set, E ∞ HF , α, β, A i , B i and C i are fitted parameters. The extrapolation schemes should use at least energy points obtained at up to the cc-pVQZ basis sets. Morever, basis sets containing diffuse functions are necessary for weakly bound molecules (alkaline-earth atoms bound with He in our case). Although, CBS extrapolation usually gives accurate results for classical system ,ie, strongly bound molecules via a real chemical bond, one can notice in Section 5.5 that it is not really the case for very weakly bound van der Waals system.

Chapter 4

The CaHe X 1 Σ + state

Introduction

In this chapter, we deal with different electronic structure calculations for the CaHe ground state potential energy curves. These calculations employ large sets of basis functions centered on He atom and Ca atom with addition of diffuse mid-bond functions to better describe the electronic binding interactions. Appropriate ab initio methods, such as MP2, MP4, and CCSD(T), are used to determine the van der Waals potential energy curves for the CaHe molecule. The contributions of core-valence correlation effects and mid-bond functions effects have been studied. Also, the basis set superposition error (BSSE) has been taken into account using the counter poise approximation of Boys and Bernardi [56]. Hence, all energies reported in this chapter are counter-poise corrected.

Computational details

The binding interaction between the closed shell Ca and He (1s 2 ) atoms belongs to the range of very weakly bound van der Waals systems. Because of this very small interaction the CaHe ground state system requires an accurate and appropriate ab initio treatment. All calculations have been carried out in the C 2v symmetry point group by means of the Molpro suite of programs [START_REF] Werner | Molpro 2006.10 is a package of ab initio programs[END_REF]. We first performed electronic structure calculations at the coupled cluster level including explicitly single and double excitations and a perturbative treatment of triple excitations (denoted as CCSD(T)) [52] in combination with several kinds of Dunning type basis sets (aug)-cc-pVXZ where X=D, T or Q. Then we calculate the interaction energy of our system using the second and fourth-order Møller -Plesset perturbation method, respectively symbolized by MP2 and MP4 [48]. All ab initio elec-tronic potential energy points fall into the distance range from 6.0 a 0 (Bohr) to 35.0a 0 with variable intervals so that the PES is well represented around the minimum. The basis sets used in this work were (aug-)cc-pCVXZ for the calcium atom (X=D,T,Q) [START_REF] Koput | [END_REF] and aug-cc-pVXZ for the helium atom (X=D,T,Q,5) [64]. However, X should reach a minimum value in order to obtain smooth potential energy curve (see Section 4.2.1). After this we treat the CaHe molecule with the CCSD(T) method in combination with large and diffuse basis set : aug-cc-pV5Z for the helium atom and the well adapted aug-cc-pCVQZ for the calcium atom. Furthermore interatomic bond functions [65] have been used in order to better describe the binding interaction. These interatomic functions have been placed at the geometrical center of the system, however an other alternative approach can be used : placing the bond function at the center of mass of the molecule. This latter method was not used in this work since following the arguments in Ref. [66], no real difference is found between these two choices. In the present work, all electrons -except the (1s 2 ) of the calcium atom-were correlated. For all electronic correlation calculations (CCSD(T), MP2, MP4), the starting point (the reference determinant) was the restricted Hartree-Fock (RHF) determinant.

We define binding energies of the system, ∆E, as the difference between the CaHe van der Waals molecule electronic energy in its ground state( 1 Σ + ) and these of the atoms in their ground states ( 1 S 0 for both atoms). ∆E is defined as follows :

∆E = E CaHe -E Ca -E He
The next sections describe the influences of various parameters ( basis size, bond functions, ...) on the CaHe ground state calculations.

Influence of basis sets size

The following figure 4.1 shows the potential energy curves as a function of the interatomic distance for the CaHe molecule determined at the CCSD(T)/cc-pVDZ and CCSD(T)/cc-pVTZ levels. It is well known, especially in the case of variational methods, that the electronic energies increase in absolute value with respect to the range X of the Dunning type basis set, (aug)-cc-pVXZ. As the variable X which equals the number of exponents ζ in the definition of the gaussian basis set increases, the number of basis functions becomes larger. Thus the energy tends to the exact value, which is simply implied by the variational principle. The purpose is then to see if there is any binding energy at small nζ. It is also interesting to know how the potential binding energy varies against the basis set used. So we have made calculations at the CCSD(T) level of theory first using Secondly, cc-pVTZ basis for the helium atom and the cc-pCVTZ one for the calcium atom have been utilized in combination with the same set of bond functions. The potential energy curves are depicted in figure 4.1. During these calculations, all the core (1s 2 2s 2 2p 6 3s 2 3p 6 ) electrons of Ca have been kept frozen. The double ζ curve is characterized by a well depth of 0.4 cm -1 and an equilibrium distance of 13.16 a 0 , whereas examination of the triple ζ curve gives 1.16 cm -1 as well depth and 12.27 a 0 as equilibrium distance. The conclusion is that the double or triple ζ basis sets, especially when diffuse function are absent, are not able to properly describe the electronic interaction for this type of system, since the final CCSD(T) calculations (where aug-cc-pV5Z and aug-cc-pCVQZ have been used for respectively He and Ca in addition of the 33211 bond functions) give 3.27 cm -1 and 11.26 a 0 for respectively the well depth and the equilibrium distance of the CaHe ground state. In the other hand, we have performed standard calculations at the same level of theory (CCSD(T)) with the same basis functions except that the interatomic set of functions were absent and the active space was constructed using only valence electrons of He and Ca. But in this case, calculations have not converged in a reasonable threshold in energy which was equal to 10 -7 a.u. These results should not be interpreted as very surprising since our system is a very weakly bound van der Waals systems.

Bond-functions role

.2 graphically illustrates the influence of the mid-bond functions on the CaHe ground state binding energy. Bond functions are known to improve the description of the electronic interaction when one deals with very weakly interacting systems. Although modest basis sets such as aug-cc-pVDZ+BF and aug-cc-pVTZ+BF basis sets yield reasonable estimates of bond energies, in most cases, these results cannot be considered highly accurate [67]. On the other hand bond functions (BF) with large and diffuse basis sets with Nζ = 4 or 5, for example, tend to give accurate binding energies. In order to see how the potential evolves upon changing the number and type of bond functions, we have made calculations at the CCSD(T) level using aug-cc-pV5Z for helium atom and aug-cc-pCVQZ for calcium atom. As we can notice in table 5.2, the well depth increases when increasing the number of bond functions (listed as BFs). Thus without BF, the equilibrium potential energy is -14.90 µHartree (3.27 cm -1 ) at bond length of 11.26 a 0 (5.96 Å ) whereas a complete set of functions 3s3p2d1f1g (listed as 33211) gives -15.09 µhartree (3.31 cm -1 ) for a binding distance of 11.25 a 0 (5.95 Å ). Thus the potential well depth becomes deeper when adding bond functions. The use of an intermediate size of the bond function set (332) gives an equilibrium distance of 5.95 Å and a corresponding

-3.4 -3.3 -3.2 -3.1 -3 -2.

Comparison of methods

As seen in section 4.1, several electronic structure calculation methods have been used in order to determine the X ).

The notation aCVQZ and aV5Z are respectively for aug-cc-pCVQZ and aug-cc-pV5Z. The Ca core is 1s 2 .

equilibrium distances, they are in quite good agreement. In our case, the equilibrium distances and corresponding energies of the CaHe ground state complex are reported in table 4.2. One can notice, for the three methods (CCSD(T), MP2 and MP4), that the equilibrium distances stand in the range 5.83-5.95 a 0 , whereas the equilibrium energies stand in the range 3.31-4.24×10 -6 a.u. (in absolute value). This makes us think that the true potential is situated in the ranges previously mentioned. We can advance this assumption because the results deduced from all the methods are consistent with each other and do not differ by an important factor.

On the other hand, in the case of the most accurate method, at least in our opinion, ie, CCSD(T), the potential energy surface's characteristics (r e and ǫ) are in excellent agreement with those reported in literature (see Section 4.5). On the other hand the most accurate ab initio energy points are given in table 4.3.

Determination of dispersion coefficients

When we deal with extremely weak interaction as it is the case for the CaHe van der Waals complex, we have to keep in mind that the results we get must be verified in any case. One way to verify the coherence of the ab initio potential energy curves given by the different methods used here is to extract from them the dispersion coefficients and study the difference between them and those given in the literature. It is assumed that interaction between neutral closed shell atoms is mainly governed by dispersion energy [69,70] which can be expressed as :

V (r) = - ∞ k=3 C 2k r 2k
where r describes the interatomic distance and C 2k is the 2k th order dispersion coefficient.

If we just consider the first three terms in the equation above and neglect all others, we can rewrite the equation as :

V (r) = - C 6 r 6 - C 8 r 8 -
C 10 r 10 -... Deriving dispersion coefficients from ab initio energy points consists in fitting these points following the equation above by means of the non-linear least squares algorithm introduced by Levenberg and Marquardt. Although, the function V (r) is linear in C 2k thus a simpler algorithm is sufficient to evaluate the C 2k coefficients, we have preferred to develop a general non-linear fitting code which is useful in fitting ab initio data into a HFD functional form, as we will see later. Thus the strategy first consists in choosing an interval where dispersion interaction dominates, otherwise, C 6 , C 8 and C 10 will be affected by other phenomena such as overlaps between atomic electronic clouds. This will involve repulsive forces and then a perturbation of the dispersion interaction which are strictly attractive by definition. Concretely, we have defined our dispersion domain as starting from 20a 0 to up to 40a 0 . We think it is a reasonable choice since the minima of all potential energy curves are located around bond lengths values of 11 a 0 (see table 4.2). So we hope from r = 20a 0 the exchange repulsion is negligible. Generally, for the determination of C 6 and C 8 constants one just needs the final two calculated points of the potential energy curve. However, we have preferred to fit more that two final points. This is justified by the fact that energy point calculations become inaccurate for long distances. This well known effect is first intrinsic to the basis set used because long bond length imply very diffuse basis functions that is to say with very small exponents. Secondly, these diffuse functions often suffer from linear dependence. Finally when the two atom centers are placed too far from each other, numerical integration (among other) problems can arise in the sense that for example integrals become too small to be reliably computed numerically.

Comparison with literature

Most of the previous calculations were in conflict on the determination the CaHe potential energy curves. They predict binding energies that differ by up to a factor of 3 or more, [71][72][73]. Fortunately, the recent theoretical calculations which are summarized in table 4.4, are in relatively close agreement [66,74] with those of this work which are listed in table 4.2.

Methods r e / Å ǫ/cm -1 CCSD(T), ref [66] 5.95 3.33 Surface integral, ref [76] 5.10 10.3 CCSD(T), ref [74] 5.85 4.2 CCSD(T), ref [78] 6.04 3.14 CCSDT, ref [77] 5.9 3.43 MRCI, ref [75] 5 As can be noticed from the table 4.2, the present results for CaHe are in agreement with the results of Partridge et al. [66], obtained with the same approach (CCSD(T)). However, our results differ from those reported by Stienkemeier et al. [75] and Kleinekathöfer [76]. This proves that a contradiction still exists about the CaHe ground state energy curve. Hinde [77], using the CCSDT method, ie coupled cluster theory including single, double, and triple excitations, reported, for the CaHe complex, an equilibrium distance of 5.9 Å and a corresponding energy of 3.43 cm -1 (refer to table 4.4). These results nicely agree with our CCSD(T) results listed in table 4.2 since the equilibrium energies only differ by 0.12 cm -1 . In reference [68] C.C. Lovallo and co-workers have developed a new parameterization of the model core potential method (MCP) [79], developed by Huzinaga and co-workers, the wtMCP [80] expected to reproduce all-electron calculations nearly exactly by using a large valence basis set. Using this wtMCP pseudo-potential in combination with the CCSD(T) method, they calculated a potential well depth of 3. 4.6 as -0.69 cm -1 when the CCSD(T) interaction potential is used. The mean equilibrium distance is, for the CCSD(T) method, 7.43 Å . Table 4.6 also provides an estimation of the rotational constant which is found to be around 0.09 cm -1 for all methods.

Conclusion

Accurate binary potentials were calculated for the interaction between helium and calcium atoms using appropriate and very sensitive ab initio methods. Potential curves parameters ǫ and r e are in good agreement with the most recent results found in the literature. As a proof of a good accuracy, the dispersion parameters are in excellent agrement with those calculated in Ref. [66,68]. The radial Schrödinger equation has been solved with the CCSD(T), MP2 and MP4 aug-cc-pCVQZ+BF pair potentials and evidence of at least one vibrational level was found.

Since the coupled-cluster method with single and double excitations and a perturbational treatment of the triple excitations -CCSD(T)-is widely used to accurately calculate the interaction energies for weakly bound systems, we will therefore use the CCSD(T) potential surface as our CaHe binary potential in the Diffusion Monte Carlo simulations of Ca in helium nanodroplets.

Chapter 5

The MgHe 1 Σ + state

Introduction

In this section we present several ab initio calculations for the MgHe potential energy curves. These calculations employ large basis function sets centered on helium atom and magnesium atom. Furthermore diffuse mid-bond functions are used to better describe the interatomic electronic interactions. Appropriate ab initio methods, such as MP2, MP4, CCSD(T) and CCSDT, are employed to determine the van der Waals potential energy curves. The contributions of core-valence correlation effects and mid-bond functions have been studied in detail. Also, the basis set superposition errors have been corrected via the counter poise method of Boys and Bernardi [56]. Finally, we have tried to apply the complete basis set extrapolation (CBS) to this molecule in order to determine the binding energy for infinite ζ basis sets. Since the binding energy in the MgHe complex is mainly due to the dispersion interaction, this system belongs to the van der Waals ones and the binding energies are expected to be extremely weak. Taking into account this energy weakness, a great problem is the determination of the potential energy interaction function versus the internuclear distance. Furthermore, as, we will see in Chapter 9, particle solvation in helium droplets among other properties -in particular for alkaline earth atoms -, strongly depends on their binary interaction potential with helium atoms. For this reason MgHe potential curves must be calculated as accurate as possible.

Computational details

The binding interaction between the closed shell Mg (1s 2 2s 2 3p 6 3s 2 ) and He (1s 2 ) atoms belongs to the range of very weak van der Waals interactions. Because of this very small interaction the MgHe ground state system requires accurate and appropriate ab initio methods. As will be discussed later, the use of the CCSDT method seems to be very suitable for the MgHe system. All calculations have been carried out in the C 2v symmetry point group by means of the Molpro suite of programs [START_REF] Werner | Molpro 2006.10 is a package of ab initio programs[END_REF], except calculations which involve the use of the CCSDT method. In order to perform the latter calculations, the ACES2 code [82] has been used instead. We first performed electronic structure calculations at the coupled cluster level including explicit single and double excitations and a perturbative treatment of triple excitations (commonly abbreviated CCSD(T)) [52] in combination with several types of Dunning basis sets [55,64] listed in Table 5.1. Then we calculate the interaction energy of our system (MgHe) utilizing the second and fourthorder Møller-Plesset perturbation method, respectively symbolized by MP2 and MP4 [48].

All ab initio electronic potential energy points stand in the distance range from 7.0 a 0 (Bohr) to 26.0 a 0 with variable step so that the PES is well represented around the minimum. After this we treat the MgHe molecule with the coupled cluster method taking into account explicit treatment of triple excitations (ie, CCSDT) combined with large and diffuse basis set : aug-cc-pV5Z for the helium atom and the appropriate aug-cc-pCVQZ for the magnesium atom. Furthermore, in order to improve the bond description, the complete basis set of bond functions (33211), given in Table 4.1, has been used.

The binding energies of the system, ∆E, are defined as the difference between the MgHe van de Waals molecule electronic energy in its ground state ( 1 Σ + ) and these of the atoms in their ground states ( 1 S 0 for both atoms). ∆E is defined as follows :

∆E = E M gHe -E M g -E He (5.1)
All the ab initio energies presented in the following are counter-poise corrected. Moreover, if no specified, all the Mg electrons have been correlated except the 1s core electrons. The reference determinant for the correlation calculation methods was the restricted Hartree-Fock (RHF) determinant.

Results and discussions

In the following sections, the influence of various technical parameters (basis set sizes, core correlation effects, ...) on the geometry and energy of the MgHe X

Basis set

Table 5.1 shows, r e , the classical equilibrium bond length and the minimum of the potential energy curves of the MgHe molecule ǫ. The fundamental vibrational energy level and the mean value distance of the first vibrational level, r 0 , are also reported. All these values are estimated for a given level of theory, that is to say, an ab initio method and a set of basis functions for Mg and He completed in most cases with a set of interatomic functions; this is made in order the study the influence of atom-centered basis sets and that of the bond function (BF) as well as the consistency of the different methods used here.

Since calculations with double and triple ζ are insufficient (Chapter 4), we have directly started with quadruple ζ type basis sets. One can notice that at the "CCSD(T), cc-pVQZ + 33211 BF" level, the binding energy of the system (-16.31 × 10 -6 a.u) is not saturated since increasing the basis set for both atoms to cc-pV5Z + 33211 BF, we get -18.85 × 10 -6 a.u as potential energy minimum, that is to say a gain of -2.54 × 10 -6 a.u which could seem insignificant, but not negligible relatively to the well depth. The equilibrium bond length also changes passing from r e = 9.77 a 0 at the cc-pVQZ level to 9.63 a 0 at the cc-pV5Z level. On the other hand, we have performed calculations with the CCSD(T) method respectively using aug-cc-pV5Z and aug-cc-pVQZ for He and Mg atoms with the complete set of interatomic bond functions (33211). This gives -21.55 × 10 -6 a.u as a value for ǫ and 9.64 a 0 for r e . The same calculations have been repeated taking the cc-pVQZ basis for the Mg atom. The results were small variations of ǫ (-0.13 × 10 -6 ) and r e (0.02 a 0 ). Hence, when the complete set of bond functions is used, the contribution of diffuse basis function (aug-) becomes insignificant in front of the non-augmented basis set (cc-pVQZ). Table 5.2 provides ground state energies of MgHe determined at the CCSD(T) level of theory in combination of various basis sets.

Influence of core correlation effect

Calculations at the CCSD(T) level with the cc-pV5Z basis set for both helium and magnesium atoms, with addition of the 33211 bond function set have been carried out. First, only the 1s core orbital of the magnesium atom has been kept frozen, we found a well depth of -18.85×10 -6 a.u at the equilibrium geometry of 9.69 a 0 . Then we have only let the 3s orbitals participating to the valence space, this gives -17.60×10 -6 a.u as well depth for the potential energy curve and an equilibrium distance of about 9.81 a 0 (see Table 5.1). One can notice in absolute value a difference between the two minima of 1.25×10 -6 a.u. Although this difference may seem very weak, in case of very weakly bound system, electronic core correlation effect is relatively important to take into account. The core correlation has a less important influence on the equilibrium radius, but not negligible (around 0.1 a 0 ), however. Fig. 5.1 pictures the potential energy curves obtained by using the CCSD(T) method in combination with the complete set of bond functions and two types of basis sets : cc-pVQZ for both He and Mg and aug-cc-pCVQZ ( aug-cc-pV5Z) for Mg (He), the latter method is termed the b1 method. In both cases, the Mg core is considered to be 1s 2 or 1s 2 2s 2 2p 6 (real core). In the cc-pVQZ case, the following characteristics of the PES are determined (when the Mg core is only represented by the 1s electrons) : -16.31 × 10 -6 a.u for ǫ and 9.77 a 0 for r e whereas considering the real core for Mg, lightly different results are obtained (-15.60 × 10 -6 a.u for ǫ and 9.85 a 0 for r e ). On the other hand, similar approach have been made using the b1 method (see Fig. 5.1), the difference between the PES characteristics is clearly visible in Fig. 5.1. Therefore, it is important to correlate core electrons for the Mg atoms to better saturate the active space so that the van der Waals bond will be well described by molecular orbitals descended from atomic orbitals.

Influence of bond functions

Bond functions are known to improve the description of the electronic interaction when one deals with very weakly interacting systems [83]. Although modest basis sets such as aug-cc-pVDZ+BF and aug-cc-pVTZ+BF basis sets yield reasonable estimates of bond energies, in most cases, these results cannot be considered highly accurate [67]. On the other hand bond functions (BF) in combination with large and diffuse basis sets with (at least) Nζ = 4 or 5 for example give accurate binding energies [83]. In order to see how the potential evolves changing the number and type of bond functions, we have made calculation at the CCSD(T) level using aug-cc-pV5Z for helium atom and aug-cc-pCVQZ for Mg. As we can see in Fig. 5.2, the well depth increases (in absolute value) when increasing the number of bond functions. This influence can be quantitatively verified analysing Table 5.3. Thus without bond functions, the equilibrium potential energy is -21.38 ×10 -6 a.u at bond length of 9.62 a 0 whereas a complete set of functions 3s3p2d1f1g (listed as 33211) gives -21.72 ×10 -6 a.u for a binding distance of 9.63 a 0 . Therefore the potential well depth becomes deeper when adding bond functions. Similar changes are also noticed for the average bond length distances (r 0 ) and zero vibrational levels (E(v = 0)).

Although not necessary in case of usual systems (strongly bound), bond functions are revealed to bring an important contribution to the interaction energy (-0.31 ×10 -6 ).

PES characteristics : r 0 and ǫ

As one can see in Fig. ) a.u. for respectively MP2, MP4, CCSD(T) and CCSDT methods considering aug-cc-pV5Z basis set for the helium atom and the aug-cc-pCVQZ one in the case of the magnesium atom (see Table 5.1), in addition of the 33211 bond functions. On the other hand, all our curves converge to the value of the first atomic asymptote Mg( 1 S 0 )+ He( 1 S 0 ) which is taken as the zero energy reference. In our opinion, the suitable method is the CCSDT one : as the MP2 level of theory is known to underestimate electronic correlations [84,85] and the MP4 method over-estimate it [86]. Perturbative methods are not reliable enough for such a system whereas coupled cluster methods are to give reasonable results especially taking into account the monoconfigurational character of the MgHe ground state [86]. These methods are known to stay self-consistent, this is an advantage in determining very low binding energy.

The energy points used to picture the MgHe ground state potential energy curves (shown in Fig. 5.3) are listed in Table 5.4 of Section 5.4.

Difference between basis and C-basis set

Final calculations have been done using the aug-cc-pCVQZ basis set which contains more flexible valence orbital (sp), more polarization (dfg) and core polarization functions (spdf) and a better description of the core (s) than the aug-cc-pVQZ version. Referring to Table 5.1 and Fig. 5.4, one can notice that the CCSD(T) method in combination with aug-cc-pCVQZ ( aug-cc-pV5Z for He) basis set for Mg (for He) gives a larger well depth (-21.72 ×10 -6 at r e = 9.62 a 0 ) than that from the use of aug-cc-pVQZ basis set for Mg ( aug-cc-pV5Z for He) which is -21.55 ×10 -6 at r e = 9.64 a 0 . Therefore, using the aug-cc-pCVQZ basis set for the magnesium atom describes more accurately the potential energy of MgHe molecule since this basis contains functions which take into account the core-valence (and subvalence) electronic interactions. And this is important in this very weakly bound system.

Basis set superposition error (BSSE)

The graph 5.5 pictures the BSSE correction in the case of CCSD(T) level of calculations.

The basis sets used are aug-cc-pV5Z for He and aug-cc-pCVQZ for Mg. One can notice that this correction goes to zero for large value of the interatomic distances. Around the minimum of the potential, the BSSE correction is about 2.82 ×10 -6 , at a bond length of 9.75 a 0 (see figure 5.5). This correction could not be neglected taking into account the case of our very weak interaction. The BSSE energy correction is defined as :

E BSSE (r) = E CP (r) -E noCP (r) (5.2)
where E BSSE (r) represents the pure BSSE contribution which is defined as the difference between binding energies respectively with (E CP (r)) and without (E noCP (r)) the counter poise (CP) correction. The function E BSSE (r) is plotted in Fig. 5.5 with stars. This curve presents a smooth monotonic shape. It decreases from short bond lengths to larger ones. It tends to zero for large r values and relatively important values could be noticed for short r values. In fact this phenomenon could be easily explained considering the fact that at short interatomic distance, physically the overlap between electronic cloud of helium and magnesium is stronger than at larger distances. Since electronic density is simulated by atomic orbitals the overlap between the two atoms follows the same tendency.

Fit quality

We have fitted our CCSDT energy points by a least square algorithm we coded in fortran 90 language. The fit is very accurate since it reproduces dispersion coefficient C i ,i = 6, 8, 10 in good agreement with experimental ones as we can see in the following. Using the CCSDT energy points (see Table 5.4), we have fitted these points following an HFD-B shape. In this case the analytical potential energy curve is well expressed by the function:

V (r) = Ae -βr - 5 i=3 C 2i r 2i
(5.3)

In order to quantify the quality of the fitting process, we have evaluated the root mean square (RMS) error which can be expressed as follows [87] :

RMS = N p k=1 (V f it k -V meth. k ) 2 Np (5.4)
where Np is the number of ab initio points. For all methods, we obtain a reasonable RMS (0.03 -0.30 ×10 -6 a.u.), especially in the CCSDT calculation case. This last method presents a RMS of only 0.03 ×10 -6 a.u. This can be noticed observing Table 5.5 which also gives the fitting parameters for all the highest level of calculations performed. Table 5.4 also provides an indication on the fit quality since the relative errors between the CCSDT energy points and the corresponding fitting function, given in the last colmn, are very small. Another advantage of our fitting process is its capacity to reproduce the dispersion coefficients C 6 , C 8 and C 10 which stay in good agreement with those reported in the literature, given in Table 5.6, although in the fitting process all the fitting parameters are let free. In the case of the C 6 and C 8 coefficients, our values from the most accurate CCSDT method are respectively 1.889×10 1 and 1.279×10 3 in atomic units. The same conclusion can be hold for MP2, MP4 and CCSD(T) calculations : according to Table 5.6, all the dispersion coefficients we obtain are consistent and generally in good agreement with the literature. For comparison, Hinde [77] reported 2.0 ×10 1 and 1.1×10 3 respectively as C 6 and C 8 values. The dispersion coefficients which we have determined also nicely agree (see Table 5.6) with those determined by Standard et al. [81]. However, our C 10 is out of the range given by Standard.

Conventional CBS approximation

In this section we will study in detail the application of the Complete Basis Set (i.e CBS) approach to the MgHe complex. Logically, ab initio expectation values for any hermitian operator can be exact value only if the number of basis functions tends to infinity. And this is true when considering any method (RHF, CCSD(T), ...) and any system. Of course it is reasonable to suggest that this ideal limit, without being very pessimistic, will never be reached and will stay infeasible whatever the future computer power is. Furthermore it stays quite interesting making calculations using "infinite" basis set. The more interesting thing is, in my opinion, performing accurate calculations using a very small number of basis functions but at present accurate data are obtained using at least triple or quadruple ζ basis sets. Another mean to treat is to extrapolate energies to the complete basis set limit as we will discuss in detail in the following. The aim of this section is to determine if the CBS approximation can still be applied for extremely weakly bound systems like Mg@He n which are the main subject of this work.

r/a 0 V M P 2 V M P 4 V CCSD(T ) V CCSDT V CCSDT f it (r) ∆V (
From here, all the ab initio energies presented do correspond to the MgHe geometry of 9.5 a 0 (5.03 Å ). This distance is chosen because it is close to the equilibrium distance of 5.01 Å .

The total CCSDT energies of the MgHe molecule can be decomposed as a sum of a Hartree-Fock contribution and a correlation part : E CCSDT = E HF + E corr. . In the following sections, we present and discuss the different ways to perform the CBS approach for both the HF part and the correlation part of the energy. During the fitting process of the CBS approximation, an arbitrary weight proportional to n 3 has been applied, with n the number of ζ in the basis sets. This is done in order to reproduce as well as possible the ab initio points with high n ζ by the fitting functions.

Fitting of the HF energies

The Hartree-Fock energy contribution obtained from calculations with several basis sets, i.e aug-cc-pCVnZ for Mg and aug-cc-pVnZ for He where n takes 2-4 values and 2-6 values respectively for Mg and He, has been fitted to the following Dunning-Feller function [53,57,58] E HF (n) = E ∞ HF + αe -βn .

(5.5)

where α, β and E ∞ HF are fitting parameters. The important quantity to be remembered is E ∞ HF which represents an extrapolation the HF energy for an infinite basis.

Fitting of the correlation energies

In the system we are studying (MgHe), since its atomic components, ie He and Mg, are both close shell atoms, the major part of the binding energy is induced by the dispersion forces (London forces). Taking this fact into account, it reveals extremely important to correctly fit the correlation part of the binding energy. For this aim a choice of an appropriate fitting function is crucial, that is why we have first decided to use the conventional fitting function [60] described in the following

E cor1 (n) = E ∞ cor1 + C 1 n 3 . (5.6)
This method will be abbreviated as cor1 in the next sections. Then, we used the function reported in Ref. [61] which is said to take into account, via the 1 4 factor, higher order of inverse polynomial than 1 n 3 . This function (denoted in the following sections as cor2) is written as

E cor2 (n) = E ∞ cor2 + C 2 (n -1 4 ) 3
(5.7)

Finally, we have introduced a new fitting function, called cor3 in the next section, more flexible since it has one more adjustable parameter than the two previously presented (Eqs. 5.6 and 5.7) . This function takes the following form

E cor3 (n) = E ∞ cor3 + C 3 n β 3 (5.8)
Mg and He without BF neither BSSE (method I )

First of all, ab initio calculations have been performed for the MgHe complex, at a geometry of 9.5 a 0 close to the equilibrium bond length, using the CCSDT method in combination of several basis set sizes (aug-cc-pCVnZ for Mg and aug-cc-pVnZ for He where n takes 2-4 values and 2-6 values respectively for Mg and He). Fig. 5.6 depicts the fitting of the correlation energies. In all these figures, curves with solid lines draw the results of the cor1 method (Eq. 5.6), those with dots are from the cor2 method (Eq. 5.7) and those with dashed lines are derived from the cor3 method (Eq. 5.8). Observing the corresponding horizontal curves, one can notice that the lowest are obtained by the mean of the cor1 method (dashed lines). Referring to Table 5.7, one can notice that the most reasonable extrapolated binding energy is obtained when fitting ab initio points from method I via the cor3 fitting method. In fact an extrapolated energy of -2.65×10 -5 ± 5.85 × 10 -5 a.u. is found which can be a reasonable extrapolated value since the quadruple ζ points evaluated at the method I level gives -2.482×10 -5 a.u as binding energy (Table 5.8). However the extrapolated value is statistically insignificant because the error is too large (5.85 × 10 -5 a.u.).

Mg and He with BF and without BSSE (method II ) c) solid lines draw the results of the cor1 method (Eq. 5.6), dots those from the cor2 method (Eq. 5.7) and dashed lines those from the cor3 method (Eq. 5.8).

fitting of the correlation energies. In all the figures, curves drawn with solid lines are from the cor1 extrapolation method, curves with dots are derived from the cor2 method and finally the dashed curves are from the HF like fitting function 5.5. Generally, what can be noticed observing the horizontal curves is that the conventional cor1 method gives the best results since lowest correlation energies are obtained with it.

However no bound fundamental state is found (by fitting ab initio points from the method II ) because when observing Table 5.7 all the binding energies are positive with moreover very large error bars and this is true whatever the correlation fitting method is.

Mg and He with BF and BSSE (method III ) Fig. 5.8 picture the fitting processes of the correlation energies obtained with several basis set qualities. Curves with solid lines and dots take into account the energy points obtained using the double ζ basis sets in the fitting process. Dots and bold dashed curve are from the cor2 method that is to say that the fitting function is that described in equation 5.6. Whereas dashed lines and solid lines depict curves obtained by the cor1 method. Dashed and bold dashed lines represent the results without considering correlation energies given at the double ζ level. While horizontal curves represent the CBS extrapolated energies in each case.

The general tendency for both the atomic fragments (He and Mg) and the MgHe complex is that the lowest correlation energy is obtained when the double ζ energies are not used in combination of the cor2 method. However the balance sheet of the total correlation is not encouraging since the extrapolated correlation energies are around -6.3 × 10 -5 ± 0.008 a.u for both cor1 and cor2 methods. The first problem is that value has huge error bars, moreover, it can not compensate the HF energies (7.405 × 10 -5 a.u) enough to give a bound state. As a result, observing Table 5.7, one can see that the only extrapolated energy value (-2.65×10 -5 ± 5.85 × 10 -5 a.u.) is obtained by method I in conjunction of the correlation method cor3 (Eq. 5.8). This value is in fairly good agreement with that obtained by the various methods at the quadruple ζ level (see Table 5.8). However, the relative huge error of 5.85 × 10 -5 a.u. gives no sense to the extrapolated value (of -2.65×10 -5 a.u.). The conclusion is that the CBS schemes applied here do not provide convincing results for MgHe. And this may be generalized to very weakly bound systems. c) solid lines draw the results of the cor1 method (Eq. 5.6), dots those from the cor2 method (Eq. 5.7) and dashed lines those from (Eq. 5.8).

Correlation methods

Methods : cor1 cor2 cor3 method I 0.000196 ± 0.0065 0.00034 ± 0.005 -2.65×10 -5 ± 5.85 × 10 -5 method II 0.000109 ± 0.006 0.00032 ± 0.005 0.00052 ± 3.6 × 10 -5 method III 1.089×10 -5 ± 0.0075 1.084×10 -5 ± 0.0078 - Method III evaluated at an interatomic distance of 9.5 a 0 in combination with various basis sets.

Non conventional CBS approximation

In the section above, we have tried to extrapolate the CCSDT potential energy for the MgHe system for a bond length of 9.5 a 0 without real success in the sense that all energies for the different methods have given unreasonable values. So we discuss below another approach which consists in directly extrapolating binding energies at a given basis set level n. As in the section above, we consider that the binding energies consist of two parts : E bind. CCSDT = E bind. HF + E bind. corr. The strategy thus consists in fitting the Hartree-Fock contribution, E bind.

HF , following the same expression as in the case of atoms

E bind. HF (n) = Ae -Bn + E ∞ bind.
(5.9) Equation 5.9 supports some approximations because if the atomic energies can be fitted by a function like that of eq. 5.9, the binding energy will be a sum of exponential functions but this sum is not an exponential type function. However, in our special case the use of eq. 5.9 can be justified because the helium atom energies are negligible versus those of Mg and MgHe systems. This approximation involves that the Mg and MgHe energies are nearly identical (ie, the fitting function are not so far from each other), this allows us to fit the binding energies following expression 5.9.

Correlation methods

Methods cor1 cor2 cor3

method I -2.1891×10 -5 ± 3.9×10 -6 -2.3545×10 -5 ± 3.7×10 -6 9.7680×10 -3 method II -3.4555×10 -5 ± 5.5×10 -5 -4.3023×10 -5 ± 5.2×10 -5 3.2210×10 -2 method III -2.3298×10 -5 ± 6.5×10 -8 -2.3234×10 -5 ± 3.4×10 -8 -2.3129×10 -5 Table 5.9: CBS extrapolated binding energies. All energies are in atomic units.

On the other hand, the fitting processes of the correlation contribution to the binding energies can be made following equations 5.10 and 5.11 without any approximation. For clarity, we give again the equations that we use for the fitting process. The cor1 method uses the following equation in the fitting process :

E bind. cor1 (n) = E ∞ bind1 + C 1 n 3 . (5.10)
whereas the cor2 methods accounts for the following equations :

E bind. cor2 (n) = E ∞ bind2 + C 2 (n -1 4 ) 3
(5.11)

The third correlation method (termed cor3) is now chosen to be fitted with the same function (Eq. 5.9) than the HF binding energies.

Table 5.9 summarizes the extrapolated binding energies. Observing this table, especially the last line, one notices that the extrapolated binding energies are 5.1(1), -5.10(7) and -5.08 cm -1 for respectively the cor1, cor2 and cor3 fitting methods of the correlation part of the binding total energy, when considering the ab initio method III , ie, the method which uses the BSSE correction and BFs. These values are very similar and in good agreement with those calculated at the method III level of theory which is listed in Table 5.8 as equal to -22.96 µHartree (=-5.03 cm -1 ). However, in the case of the two other methods (method I and method II ) presented in Table 5.9, the resulting binding energies are far from the ab initio points (-5.03 cm -1 ). Moreover, large error bars are observed. For example, the method II gives an energies of -7.58 ± 12.07 cm -1 (-3.4555×10 -5 ± 5.5 × 10 -5 a.u).

Thus, the appropriate CBS method seems to be the non conventional method developed in the current section in combination with method III . J=: 24 Table 5.10: Rotational levels for main magnesium isotopes, for v = 0 and possible rotational constants for the three main magnesium isotopes. Energy units are cm -1 .

Vibrational level of MgHe ground state

Numerically solving the radial Schrödinger equation by a variational method [START_REF]The Fortran program (Vibvar) we used to integrate Schrödinger equation was writ[END_REF] using the CCSDT pair potential for the MgHe system gives us three bound rovibrational levels.

For the 3 main isotopes, Table 5.10 presents rovibrational levels with v=0 and J=0 to 2, it also provides the expectation value of the rotational constants (B 0 ) for each isotope. Compared to Ref. [START_REF] Funk | [END_REF], our rovibrational energies seem to be greater than those reported by between 20 % to more than 130 %. This could be interpreted by their shallower MP4 potential which only has 21 µHartree as well depth at an equilibrium geometry of 9.75 a 0 .

The rotational constant values (B 0 ) do not have a large dependence on the Mg isotopes and all are about 0.12 cm -1 . The MgHe ground state seems to only support a unique vibrational level (v = 0) whose energy is -0.90 cm -1 .

Conclusion

Along this chapter, we have shown that the MgHe ground state energy curve is particularly difficult to determine. Since this system belongs to the very weakly bound van der Waals systems. However, we have used the accurate CCSDT method in combination of large and diffuse basis set functions to perform calculations of our potential energy curve (PEC). Furthermore, we have taken into account several effects which influence the system energy such as the core-valence correlation (cv), the bond function (BF), etc. Doing this, we have obtained a PEC whose characteristics are in good agreement with those in literature.

On the other hand, dispersion coefficients obtained from our PEC are also in excellent agreement with those reported in literature. Finally, the rovibrational calculations carried out with the highest PEC (from CCSDT) have shown the existence of only three bound rovibrational states (v = 0, j = 0, 1, 2). Furthermore, various CBS approaches have been attempted without real success; although, the "non conventional" technique seems to give relatively satisfying results, the large error bars presented do not allow any conclusion.

Chapter 6

Introduction to quantum Monte Carlo methods

Among quantum theoretical methods, Quantum Monte Carlo methods (QMC) are becoming very interesting for highly accurate quantum calculations of both electronic and bosonic systems. QMC is a family of very diversified methods. Some of them are based on the variational principle like Variational Monte Carlo (VMC) [START_REF] Hammond | Monte Carlo Methods in ab initio Quantum Chemistry[END_REF], others on the diffusion equation such as Diffusion Monte Carlo (DMC) [START_REF] Hammond | Monte Carlo Methods in ab initio Quantum Chemistry[END_REF] and others on the path integral formalism like Path Integral Monte Carlo (PIMC) [START_REF] Ceperley | [END_REF]. By means of these methods, a large number of atomic and molecular systems have been studied [92][93][94][95][96], as well as rare-gas bosonic clusters [97][98][99]. QMC methods can almost scale linearly [100] or as a low-order polynomial with system size, N, in practice not worse than order N 3 . Thanks to this low scaling rate1 , it is possible to study large systems which are impossible to approach with other methods. This is the case of doped helium clusters which are the subject of this work. QMC methods have a wide domain of application, they are useful in physics as well as in chemistry to perform molecular electronic structure calculations. In this chapter, we will orient our attention in briefly introducing two main QMC methods: Variational Quantum Monte Carlo (VMC) and Diffusion Quantum Monte Carlo (DMC).

Variational Quantum Monte Carlo

The VMC method [101] is based on the variational principle in the sense that a wave function which depends on a set of arbitrary parameters p is formulated : Ψ(r, p). The set of parameters is then adjusted so that the energy (or its variance) of the system is minimal. If the energy derived from the previous wave function is E(p) and the exact energy is E 0 for a given state of a system, the variational theorem proves that for a set of optimum parameters p opt the energy, E(p opt ), may be a good approximation of the exact energy E 0 . In the sense that both values are linked by this relation : E(p opt ) ≥ E 0 , the equality is verified in the case where Ψ(r, p) is the exact wave function. In practice minimizing the variance of the energy expectation value rather than energy itself seems to be less troublesome and is a common way to efficiently optimize wave function.

Energy point calculation

For any hermitian operator Â, the expectation value (reported as A ) is given by the following integral equation :

< A >= Ψ| Â|Ψ Ψ|Ψ (6.1)
where Ψ is the wave function governing the system. Now, if we are interested in the ground state energy of the system, which will always be the case in our study, the interesting operator will be the Hamiltonian and equation 6.1 takes the following form

< E 0 >= Ψ| Ĥ|Ψ Ψ|Ψ (6.2)
which is equivalent in more explicit notation to

< E 0 >= Ψ * (r) ĤΨ(r)dr 3n Ψ * (r)Ψ(r)dr 3n (6.3)
where Ψ(r) is, in our case, a bosonic wave function 2 , Ĥ is the Hamiltonian operator for the system, n is the number of particles (helium atoms) in the system, and r is a 3n-dimensional vector containing the positions of all n particles. Some manipulations of expression 6.2 yield

< E 0 >= Ψ * (r)Ψ(r)E l (r)dr 3n Ψ * (r)Ψ(r)dr 3n (6.4)
We can rewrite the last equation as:

< E 0 >= p vmc (r)E l (r)dr 3n (6.5)

2 As, in our case, Ψ(r) is a real function so that Ψ * (r) = Ψ(r), we will write indifferently Ψ * (r) or Ψ(r) for the complex conjugate of Ψ(r).

where: p vmc (r) is the normalised probability for particles to occupy positions r and E l (r) = Ĥψ(r)

ψ(r)
is the corresponding energy (the local energy) for the system at these positions. Many approaches to solve equation ( 6.5) exist. One of them is the standard grid integration. But this method suffers from a high scaling of the computational cost since it scales as 2 n . Fortunately the high scaling of standard integration methods can be overcome when VMC is used. In fact VMC employs Monte Carlo integration to evaluate equation 6.5. It can be shown that Monte Carlo integration is faster than standard integration algorithms when the integral's dimensionality is greater than about 7 [START_REF] Press | Numerical Recipes in Fortran 77: The Art of Scientific Computing[END_REF].

In Monte Carlo integration, N random vectors, r k , distributed with respect to p vmc (r) are generated. The energy expectation value is then found to be:

< E 0 >≈ 1 N k=N k=1 E l (r k ) + O 1 √ N (6.6)
One can notice in the previous equation ( 6.6) that the standard deviation in the expected energy decreases with the square root of the number of samples (N). An advantage of this integration scheme is that the error is independent of the system dimensionality.

On the other hand, a disadvantage of the VMC method is that optimizing the wave function parameters is a difficult task since expectation values (energy or its variance) are stochastic quantities.

VMC wave functions

In principle any symmetric wave function may be used as the bosonic wave function for VMC calculations of helium clusters( 4 He). Nevertheless, the closer the wave function is to the exact wave function, the faster the VMC calculations converge. Furthermore an appropriate wave function should have a restricted number of parameters to simplify its optimization in order to obtain the optimal state. Generally Jastrow type wave functions may serve as one and two-body components wave function φ jas ij and the total wave function, Ψ, is constructed following the form :

ψ(r) = i<j φ jas ij (r ij ) i φ jas/f er i (r i ) (6.7)
where φ Jas ij (r ij ) is a Jastrow type function and explicitly takes the form :

φ jas ij (r ij ) = e -( c 1 r 5 ij + c 2 r 2 ij +c 3 r ij +c 4 ln(r ij )) (6.8)
and φ jas/f er i (r i ) is a Jastrow/Fermi type function which is given as :

φ jas/f er i (r i ) = e -( a 1 r 5 i + a 2 r 2 i )
1 + e a 3 (r i -a 4 ) (6.9)

In the previous equations r ij refers to the distance between atom i and atom j, r i to the distance of particle i from the overall center of mass. The a i and c i are constant parameters to be optimized.

Metropolis algorithm

It is almost impossible to dissociate the Metropolis algorithm from QMC methods because of its important application in the acceptance probability of a move attempt. This algorithm [START_REF] Metropolis | [END_REF] starts by the following equation for the time evolution of a density

∂ρ(r, t) ∂t = [T (r ′ → r)ρ(r ′ , t) -T (r → r ′ )ρ(r, t)]dr ′ (6.10)
In equation 6.10, the function ρ(r, t) is any probability distribution at time t. The function T (r → r ′ ) represents the transition probability for moving from position r to position r ′ and is normalized so that T (r → r ′ )dr ′ = 1. Furthermore, if the system is in equilibrium we have a stationary density with ∂ρ(r,t) ∂t = 0. The equilibrium assumption implies a time independent character of the system so that ρ(r, t) is just r dependent (ρ(r, t) = ρ(r)). Taking into account these two assumptions transforms equation 6.10 into the simpler one below : 0 = [T (r ′ → r)ρ(r ′ ) -T (r → r ′ )ρ(r)]dr ′ (6.11) Equation 6.11 has an infinite number of solutions but the simplest and most used solution invokes detailed balance. This solution is based on the fact that an equivalence exists between f (r)dr = 0 and f (r) = 0. If we let f (r) = T (r ′ → r)ρ(r ′ ) -T (r → r ′ )ρ(r), we obtain the following relation, also known as the detailed balance solution :

T (r ′ → r)ρ(r ′ ) = T (r → r ′ )ρ(r) (6.12)
From equation 6.12, the probability for accepting a move attempt from position r to position r ′ is given by

A(r ′ , r) = min 1, T (r ′ → r)ρ(r ′ ) T (r → r ′ )ρ(r) (6.13)
In eq. ( 6.13), it should be noted that in most implementations of the Metropolis algorithm, the transition3 probability T is chosen so that T (r ′ → r) = T (r → r ′ ). Now, we are able to produce random numbers distributed with respect to any given distribution ρ(r) even if the distribution is complicated. The random numbers distributed with respect to ρ(r) are the samples r i . The algorithm can be summarized by the following steps :

• choose an initial point r 0

• choose a random point r ′ i to generate the (i + 1) th random sample.

• generate a uniform random number ξ

• if ξ > A(r ′ i , r i ), r i+1 ← r i else r i+1 ← r ′ i
Repeating this process a certain number of times finally gives the desired random numbers distributed according to ρ(r). The initial point r 0 should, however, be chosen carefully so that the ρ(r) distribution is rapidly sampled.

Diffusion Monte Carlo

The disadvantages of the VMC method can be overcome by the diffusion quantum Monte Carlo (DMC) method at the cost of a more expensive simulation. The DMC method [104] has been extensively discussed in the literature [START_REF] Hammond | Monte Carlo Methods in ab initio Quantum Chemistry[END_REF]105,106]. In this section we aim at summarizing the main features of the method.

Why diffusion?

First of all, we can wonder about the meaning of the word "diffusion". In fact the key idea of the DMC method is the analogy between the time dependent many body Schrödinger equation and a classical diffusion equation with anisotropic diffusion coefficients. The time dependent many body Schrödinger equation for a system composed of n particles is written as

ih ∂ψ( r, t) ∂t = - h2 2 n j=1 1 m j ∇ 2 j + V ( r) ψ( r, t) (6.14)
where r is a 3n-dimensional vector which stores the system coordinates and m j the mass of particle j. Now, if we introduce an imaginary time variable τ = it/h and shift the absolute energy scale by an arbitrary quantity E ref , we obtain the following equation :

∂ψ( r, τ ) ∂τ = n j=1 h2 2m j ∇ 2 j -(V ( r) -E ref ) ψ( r, τ ) (6.15)
The multidimensional reaction diffusion phenomenon for n particles is governed by the following relation

∂C( r, t) ∂t = n j=1 D j ∇ 2 j -k( r) C( r, t) (6.16)
If we identify the inverse mass terms with diffusion coefficients D j and the shifted potential V ( r) -E ref with the position dependent rate terms k( r), the analogy between equation 6.15 and 6.16 becomes obvious.

DMC method

Diffusion Quantum Monte Carlo (DMC) can give exact expectation values for N-body quantum problems within statistical errors. The DMC method solves the time dependent Schrödinger equation exploiting the isomorphism between the latter and the multidimensional diffusion equation. The Hamiltonian operator, Ĥ, for a system of n particles is described by the formula below

Ĥ = - h2 2 n i 1 m i △ i + { V -E ref } (6.17)
where V is the potential energy and E ref a shifting energy constant. The derivatives involved by the Laplacian operator, △, are with respect to the 3n coordinates. The Schrödinger eigenvalue equation satisfies the following equation

Ĥ|Φ k (r) = E k |Φ k (r) (6.18)
E k and Φ k (r) are the eigenvalue and eigenfunction pairs. From the Hamiltonian in equation 6.17, it is possible to construct a new Hamiltonian, H, which has eigenvalue-eigenfunction pairs of E k and ϕ k = Φ k (r)Ψ trial (r) and where Ψ trial (r) is a guiding wave function.

H = n i -1 2m i △ i + 1 m i ∇ i • ∇ ln |Ψ trial (r)| + E l (r) (6.19)
The local function E l (r) = ĤΨ trial (r) Ψ trial (r) has the same meaning as in the VMC case. It is the local energy of the system for given positions of particules. The time-dependent Schrödinger equation for H in imaginary time τ is

- ∂ ∂τ |ϕ(τ ) = ( H -E ref )|ϕ(τ ) (6.20)
Assuming that the Hamiltonian H is time-independent, equation 6.20 has the formal solution :

|ϕ(τ ) = e -( H-E ref )τ |ϕ(0) (6.21)
Expanding the formal solution in terms of the eigenfunctions of H gives

ϕ(r, τ ) = k c k e -(E k -E ref )τ Φ k (r)Ψ trial (r) (6.22)
where c k = Φ k (r)Ψ trial (r)|ϕ(r, τ = 0) . One can notice that the high-energy components in expansion 6.22 die out exponentially with τ and : which is proportional to |Φ γ Ψ trial . Therefore, random coordinates generated with the distribution ϕ(r, τ → ∞), are also distributed with respect to |Φ γ Ψ trial . This is what is required to evaluate the DMC energy using Monte Carlo integration.

lim τ →∞ |ϕ(r, τ ) = lim τ →∞ c γ e -(Eγ -E ref )τ |Φ k (r)Ψ trial (

Time evolution and Green's function

Equation 6.21 can not be exploited by the DMC algorithm in its present shape and must be expressed in terms of a Green's function :

ϕ(r ′ , τ ) = G(r ′ , r, τ )ϕ(r, τ = 0)d 3n r (6.25)
where

G(r ′ , r, τ ) = r ′ |e -( H-E ref )τ |r (6.26)
is the Green's function for the problem. In most interesting physical problems, it is impossible to correctly evaluate G(r ′ , r, τ ) for arbitrary τ . Only a small amount of physical problems has an analytical Green's function solution. However, in the small time step approximation, dτ , G(r ′ , r, dτ ) can be factored 4 into two independent parts, via the Trotter formula [107], easy to evaluate.

G(r ′ , r, dτ ) = G dif f (r ′ , r, dτ )G branch (r ′ , r, dτ ) + O(dτ 2 ) (6.27)
G dif f (r ′ , r, dτ ) is the "diffusion" function which describes the probability of moving from a point r to r' in dτ imaginary time [98]

G dif f (r, r ′ , dτ ) = j    m j 2πdτ 3 2 exp   - m j 2dτ r j -r ′ j - dτ 2m j F j (r) 2      (6.28)
where m j is the mass of particle j and F j (r) = Ψ trial ∇ ln Ψ trial is a quantum drift force. The branching part of the Green's function, G branch (r ′ , r, dτ ), is a factor evaluating how the value of ϕ changes in going from (r, τ ) to (r ′ , τ + dτ )

G branch (r ′ , r, dτ ) = exp -dτ ef f E l (r) + E l (r ′ ) 2 -E ref (6.29)
Under the small time step approximation, equation 6.25 is transformed into a recursive form ϕ(r ′ , (k + 1)dτ ) = G(r ′ , r, dτ )ϕ(r ′ , kdτ )dr 3n + O(dτ 2 ) (6.30)

Iteratively applying equation 6.30 offers a good approximation of ϕ(r, τ ) for large imaginary time τ . Equation 6.30 is 3n-dimensional, which represents an equation of about 600 dimensions if a helium cluster composed of 200 atoms is simulated. Thus Monte Carlo integration is a reasonable way to determine the integral. In the stochastic evaluation of this integral, a correspondence can be established where ϕ(r, τ ) is represented by

ϕ(r, τ ) ≈ j ω j (τ )δ(r -r j (τ )) (6.31)
where ω j (τ ) represents a statistical weight and δ(r -r j (τ )) is the well-known Dirac delta function centered at r j (τ ) . The couple {r j (τ ), ω j (τ )} is known as a walker. Combination of equations 6.30 and 6.31 leads to ϕ(r ′ , (k + 1)dτ ) ≈ j ω j (kdτ )G dif f (r ′ , r j (kdτ ), dτ )G branch (r ′ , r j (kdτ ), dτ ). (6.32) 4 Other schemes exist for factoring the Green's function and for recovering the delta function representation of ϕ(r ′ , (k + 1)dτ . The details of these algorithms and their advantages and disadvantages are covered in the literature [START_REF] Hammond | Monte Carlo Methods in ab initio Quantum Chemistry[END_REF]93]. Umrigar's algorithm [93] seems to be the most used since it is stable, and has a small time-step bias.

Equation 6.32 can be transformed into the delta function form. To do this, each of the new delta function positions, r j ((k + 1)dτ ) , is randomly chosen from the distribution G dif f (r ′ , r j (kdτ ), dτ ). The new weights are then given by ω j ((k + 1)dτ ) = G branch (r ′ , r j (kdτ ), dτ )ω j (kdτ ) (6.33) This new set of walkers is a stochastic representation of ϕ(r, (k + 1)dτ ). The new set of random points, r j ((k + 1)dτ ), given the appropriate statistical weights, ω j ((k + 1)dτ ), are random points distributed with respect to r j ((k + 1)dτ ). By choosing ϕ(r, τ = 0) to be |Ψ trial (r)| 2 , a stochastic representation of ϕ(r, τ = 0) can be generated by setting ω j (τ = 0) = 1 and r j (τ = 0) equal to random points generated with respect to |Ψ trial (r)| 2 using the Metropolis algorithm. After many applications of equation 6.30, the walkers will provide a stochastic representation of ϕ(r, τ → ∞), which equation 6.24 showed to be proportional to the distribution we are trying to sample, Φ 0 (r)Ψ trial (r). This produces the samples required to evaluate the DMC energy.

Move acceptance

As implied above accepting a move attempt is made by the Metropolis criterion following the probability :

P (r → r ′ ) = min{1, A(r → r ′ )} (6.34)
where

A(r → r ′ ) = |ψ trial (r ′ )| 2 G(r ′ → r) |ψ trial (r| 2 G(r → r ′ ) (6.35)
The asymmetric transfer function G(r ′ → r) has to be explicitly taken into account in this acceptance5 decision. In the branching part of this function, eq. ( 6.29), the effective time step δτ ef f is defined through the ratio of accepted displacements and attempted displacements according to

δτ ef f = < δ r acc > 2 < δ r att > 2 (6.36)

DMC wave function

A common ansatz for atomic clusters and bulk systems (see for example [98]) expresses Ψ trial as a product over a set of one-dimensional functions φ jas ij defined over all pairs of particles. Just like in VMC this result is multiplied by the product of φ jas/f er i

Ψ trial (r) = n i<j φ jas ij (r ij ) n i φ jas/f er i (r i ) (6.37)
where r ij is the distance between particles i and j, r i is the distance between atom i and the center of mass.

In our case of an alkaline-earth (Ca or Mg) atom embedded in an atomic helium environment, we used a guiding function (Ψ trial ) determined in equation 6.37. We write Ψ trial as a product over pair wise radial functions φ jas as above connecting the helium atoms and a product over isotropic functions φ jas/f er which make hold the cluster together by an attractive long range tail.

DMC Energy Evaluation

The DMC energy can be evaluated by two main way: the local energy estimator which will be presented first and the growth estimator.

Local energy estimator

The DMC energy, E dmc , is evaluated using a mixed estimator :

< E dmc >= Φ 0 (r)| Ĥ|Ψ trial (r) Φ 0 (r)|Ψ trial (r) (6.38)
< E dmc >= Φ 0 (r) ĤΨ trial (r)dr 3n Φ(r)Ψ trial (r)dr 3n (6.39)

In the equation above, |Ψ trial (r) should be an approximation to the desired ground state |Φ 0 and r is a 3N-dimensional vector containing the coordinates of the N particles. Since Ĥ is an Hermitian operator and both |Φ 0 and |Ψ trial are real, Φ 0 | Ĥ|Ψ trial = Ψ trial | Ĥ|Φ 0 , it can be shown that E dmc = E 0 . Rearranging equation 6.38 yields to :

< E dmc >= Φ 0 (r)Ψ trial (r)E l (r)dr 3n Φ 0 (r)Ψ trial (r)dr 3n (6.40)
where E l (r) = ĤΨ trial (r) Ψ trial (r) is the local energy of the particle positioned at coordinates r. In practice, as shown in section 6.3.3, the instantaneous weights are a stochastic representation of the ground state wave function. Since in the DMC method, only expectation values of local operators are directly accessible, the integration reduces to an average over operator values. Here we are interested in the ground state energy (E dmc ) which is approximated by :

< E dmc >= N i=1 ω i E l (r i ) Ω (6.41)
where E l is the local Hamiltonian operator (the local energy previously defined), ω i the current relative weight of the random walker i or its weight multiplied by the sum of descendant weights [108] and Ω is the total sum of walkers' weights

Ω = N i=1 ω i (6.42)
Use of the random walkers' instantaneous weights amounts to a "mixed" expectation value Ψ trial | Ô|Ψ . From the mixed DMC expectation values, one can compute the extrapolated quantity

Ψ| Ô|Ψ = 2 Ψ trial | Ô|Ψ -Ψ trial | Ô|Ψ trial + O(∆Ψ 2 ) (6.43) ∆Ψ = Ψ -Ψ trial
but this approach can give bad DMC expectation values and its results are not taken into account in our final results. We rather use the collection of descendant weights to generate a statistical estimate of the true |Ψ| 2 distribution [108].

Growth estimator

Another way to determine the DMC energy is the use of the growth estimator method. The latter has the advantage of presenting a smaller time step error [98] compared to local energy estimator method, discussed in the previous section (6.3.6). The imaginary time evolution of the walker i is determined by the following equation

ω i (r i , τ + δτ ) = exp{-(E l (r i ) -E ref )δτ }ω i (r i , τ ). (6.44)
The growth energy, E growth , is proportional to the logarithmic derivative of the total weight with respect to imaginary time

E growth = E ref - d ln Ω(τ ) dτ (6.45)
where Ω(τ ) is the sum over all the weights of walkers Ω(τ ) = i ω i (6.46)

Error analysis

Error analysis of data is central to Monte Carlo simulations. For this reason, statistical techniques used to compute the expectation values and their error bars are presented in the following sections.

Correlated samples

Quantum Monte Carlo simulations (VMC as well as DMC) calculate expectation values of physical systems, < X >, usually using a high-dimensional probability distribution function, p(x), in combination of a multidimensional integral.

< X >= p(x)φ(x)dx (6.47)

Although the previous expression seems simple and compact, in most physical interesting cases, it is too time consuming to directly solve equation ( 6.47) using modern computers. Instead, it is typically useful and more convenient to calculate the (time) average of X, denoted as X.

X = 1 N i=N i=1 φ(x i ) (6.48)
In equation 6.48, the index i is related to the Monte Carlo step number, and x i is a 3n-dimensional coordinate vector which is sampled from the distribution p(x). Assuming ergodicity, < X > and X can be related through the following relationship

X = lim N →∞ X = 1 N i=N i=1 φ(x i ) (6.49)
To approach equality in equation 6.49, very large sample sets should be used. But, in any case, since the number of samples (N) is necessarily finite, X will fluctuate as the calculation progresses because of its non-zero variance V ar( X). This variance can be written as

V ar( X) = 1 N N i,j=1 [ φ(x i )φ(x j ) -φ(x i ) φ(x j ) ] (6.50)
In the simple case of uncorrelated data ({φ(x i )}), the covariance terms are zero, and equation 6.50 reduces to the well-known variance relation :

V ar( X) = X 2 -X 2 N = σ 2 (X) N (6.51)
Since DMC (as well as VMC) calculations are based on a random walk [START_REF] Metropolis | [END_REF] they produce {φ(x i )} with non-zero covariances which means a serial correlation between data.

In fact the probability of picking x i depends on the value of x i-1 . (see section 6.2). The Metropolis algorithm generates by this way a Markov chain [START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF]). If one utilizes equation 6.51 to calculate the variance of such markovian systems, the result will be incorrect and statistical error bars are underestimated because the covariance between samples is missing.

Correlation analysis

To decorrelate data given by a Markov process, one has to take into account the correlation length of data.

Autocorrelation function

Since of the Markov chain data are serial correlated, the determination of the correlation time, τ c , is crucial. The latter gives an estimate of how long the system takes to evolve between states that are practically independent. The autocorrelation function of the data, γ(t), which is given here in its normalized form,

γ(t) = +∞ -∞ f (t ′ )f (t ′ + t)dt ′ +∞ -∞ f 2 (t ′ )dt ′ (6.52)
is required to define the correlation time. For discrete data, the autocorrelation function, given by Eq. 6.52, is estimated by the following expression :

γ(t) = γ(jδτ ) = Cov(x i , x i+j ) V ar(x i ) (6.53)
which is explicitly computed as

γ(jδτ ) = M -j i=1 (x i -x)(x i+j -x) M -j i=1 (x i -x) 2 (6. 54
)
where δτ is the DMC time step, M is the number of samples and x is the average of the

x i x = 1 M M i=1
x i (6.55)

The important quantity to be computed in order to correct the DMC energy variance is the correlation length, that is to say τ c . This is defined as :

τ c = 2 ∞ 0 γ(t) 2 dt (6.56)
In our case, the autocorrelation function is defined such that τ c = 1 for uncorrelated data and τ c > 1 for correlated data. The size of the correlation length indicates the degree of correlation between data. In practice, since the autocorrelation function γ(t) is computed following its discrete formula (equation 6.53), τ c is evaluated via equation 6.56 using a numerical integration method such as Simpson's formula.

The DMC case

Considering a DMC expectation value < X >, typically energy, for a given system, in practice, for each block, we calculate < X > as a mean over time steps. The standard deviation of the observable < X > is then given by the following relation [START_REF] Newman | Monte Carlo Methods in Statistical Physics[END_REF]

σ b = 1 M b -1 M b j=1 (< X 2 j > -< X > 2 b ) , (6.57)
In Eq. 6.57, M b is the total number of blocks. < X j > and < X > b are respectively averages over block j and over all blocks. The standard deviation, σ b , is expected to increase with block length until a limiting value which corresponds to the true standard deviation if samples were not correlated. Since the latter are necessary correlated, σ b is finally corrected by multiplication by a factor which is the root square of the correlation length, τ c to give the "true" standard deviation of the expectation value, σ σ = √ τ c σ b . (6.58)

Statistical errors

In the previous section, we have presented the way to properly determine error bars, ie, statistical errors committed during simulation. In the following we list the main origin of statistical errors encountered during our DMC simulations :

• insufficient number of walker.

• insufficient number of blocks.

• bad guiding function (trial wave function)

Systematic errors

The statistical techniques previously presented are unable to detect nor correct any systematic errors below, except in the first case:

• Average with non-relaxed data.

• Time step error due to the approximation for the Green's function.

• Too short simulation.

• not enough random walkers.

Calculation of main properties

The current section presents the different properties which are necessary in visualizing the helium droplets. We also discuss technical details which allow their computation. All histograms shown in the figures contained in the chapter on DMC results have been computed with descendant weights [108].

Radial distribution

The descendent weights approach is in particular applicable to the positional correlation function which is very useful in visualizing the structure of the clusters. The radial distribution of helium atoms relative to the center of mass of the whole cluster is computed as

P rad (r) = 1 n n i δ(r i -r) r 2 walk (6.59)
This quantity is very interesting in judging the density of a helium cluster compared to the bulk helium density. We can also qualitatively determine the magnesium or calcium position in the helium droplet when visualizing the density histogram of the Mg (or Ca) atom.

Pair correlation function

In a similar way we compute the discrete version of the helium pair correlation function according to The radial distribution function can be easily converted to the spherically averaged radial density distribution ρ(r) using

P pair (r) = 2 n(n -1) n i<j δ(r ij -r) r 2
n = 4π ∞ 0 ρ(r)r 2 dr (6.62)
and then ρ(r) = 4π n P rad (6.63)

Two-dimensional histograms

In a similar way we compute two-dimensional histograms in cylinder coordinates to analyze the density distribution ρ(r, z) of helium around the magnesium or calcium atom. The z-axis is defined by the center of mass of the droplet and the alkaline-earth dopant.

On the other hand the perpendicular distance of helium atoms to this axis defines the polar radius r. The origin coincides with the center of mass of the droplet and Mg or Ca is on the positive z-axis. The density distribution is computed by the following expression

ρ(r, z) = n 2π n i δ(r i -r) r δ(z i -z) walk (6.64)
The quantity ρ(r, z) is accumulated on a grid which is equidistant in z and r 2 which eliminates the need to take square roots during the data collection.

The number of helium atoms, n, can be found by integration of the density distribution ρ(r, z) following the expression

n = 2π ∞ 0 ∞ -∞ ρ(r, z)rdrdz (6.65)
and was used as a check for the completeness of the histogram.

Pseudo-codes

In the following paragraphs, we present Fortran pseudo-codes to illustrate the implementation of VMC and DMC methods adapted from [START_REF] Shumway | [END_REF]. These codes aim at presenting the main steps of both methods. In reality, the implementation is a little bit more difficult while the algorithms are easy to understand.

Conclusion

In short, we have seen in the current chapter that the VMC method can reasonably determine some properties of the system (energy, structure,...) with a modest calculation cost, but the main disadvantage is the strong dependence of the results on the quality of the wave function. In the sense that a bad wave function systematically provides bad expectation values and this, whatever the simulation length is. These inconveniences can be overcome by the DMC method which exactly solves the many-body Schrödinger equation (within statistical errors). But this is done at the cost of longer simulations because DMC is more time consuming than VMC.

Chapter 7 Doped helium nanodroplets 7.1 Introduction

Helium nanodroplets present interesting characteristics. Some of the most important questions to be clarified are : Are the droplets superfluid? What is the temperature of these systems in the extreme conditions where they are produced? Some elements of answers are detailed in the next sections. However, an interesting question can be about the shape of the droplets which are generally said as spherical but since the droplets are produced with non-zero angular momentum, the spherical form is not especially preserved. The main properties as well as some applications are given in the following.

4 He nanodroplet properties 7.2.1 Superfluidity

One of the first interesting experiments on doped helium clusters was performed by Scoles et al. [START_REF] Goyal | [END_REF] during the year 1992. They mesured infrared absorption spectra of the SF 6 monomer and noticed a line width much smaller than in other rare gas matrices. After this, Toennies et al. [12] measured a similar spectrum of the same system and it appears that SF 6 freely rotates in the droplets (a well defined rotational structure was observed). The same conclusion was found by Grebenev et al. [16] for the OCS molecule.

The final demonstration of the superfluid character of the helium ( 4 He) droplet is attributed to Grebenev during the year 1998. In fact Grebenev et al. [16] performed an IR spectroscopic experiment of both 3 He and 4 He systems doped with the OCS molecule. The spectra of OCS in 3 He and 4 He nanodroplets show very different shapes in these two systems. In Fig. 7.1, from Ref. [16], the spectra measured present a well resolved rotational structure (spectra (A) corresponding to the bosonic system) whereas the spectra implied by the 3 He system only shows a broad band. The difference demonstrates the superfluid character of the ( 4 He) droplets since 3 He is not superfluid at this temperature (0.4 K).

The superfluidity of the helium droplet is confirmed by theoretical calculations. Pure and doped 4 He nanodroplets are supposed to be superfluid at a temperature of about 0.4 K for a minimal number of 4 He atoms [START_REF] Szabo | Modern Quantum Chemistry[END_REF] [17,18] or at least have a significant superfluid component.

Moreover, Grebenev et al. give an indication on the minimum size of the helium 4 droplet which make it superfluid. It appears that quasi-free rotation is observed (see figure 7.2) when about 60 4 He atoms solvate the OCS molecule. So manifestation of superfluidity is provided by at least 60 helium 4 He at this operating temperature (0.4 K).

Temperature of the droplets

Hartmann et al. [12] experimentally studied the rotational structure of the infrared spectra of the SF 6 molecule. Thanks to this spectrum the temperature of the droplets was determined to be 0.37 ± 0.05 K. The same temperature (0.37 ± 0.02 K) has been determined by Grebenev et al. [15] using the IR spectra of the OCS molecule embedded in large 4 He droplets.

On the other hand, Nauta et al. [14] deduced a temperature of 0.34 ± 0.01 K from the rotational spectra of the HF dimer in 4 He clusters and Callegari and co-workers found a similar temperature of 0.348(2) K when they measured the IR spectra of cyanoacetylene in 4 He droplets.

The light difference between the temperature experimentally obtained can be explained by the fact that the temperature depends on the size of the droplets [13].

The experimental values for the droplet temperature are in good agreement with the theoretical prediction of 0.32 K [START_REF] Gspann | Physics of Electronic and Atomic Collisions[END_REF]11].

Experimental aspects 7.3.1 Production of helium nanodroplets

Actually, many ways exist for the formation of helium droplets with sizes starting from a few helium atoms to macroscopic aggregates. These helium droplets are often formed by supersonic expansion of helium gas but many other approaches can be applied in order to accomplish this aim. Producing the helium droplets is a real challenge, at least at the beginning, because of the extreme working conditions : very low temperatures, high helium pressures, etc. However, nowadays, the supersonic expansion technique can be divided into two categories : -continuous helium droplet beam [START_REF] Goyal | [END_REF]5,[START_REF] Scheidemann | [END_REF] -pulsed helium droplet beam [113]. In both cases physical characteristics of the apparatus must be well controlled in order to obtain droplets in a given size range. The helium droplet formation is well discussed in Ref. [19] and references inside. The droplet size depend on the initial conditions, pressure and temperature, at the nozzle (respectively P 0 and T 0 ). The nozzle diameter also has an effect on the droplets sizes.

Doping of droplets

Doping the helium droplets is generally achieved by means of collisions between atomic (or molecular) impurity and the helium droplet beam in a scattering cell. This method is known as the pick-up technique [START_REF] Scheidemann | [END_REF]114]. In the case of helium clusters doped with Mg atom, the doping process is performed by heating metallic Mg until a partial pressure is reached in magnesium (typically between 0.01 to 0.1 Pa [28,115], depending of the mean number of Mg wanted in the droplets) within the doping cell (see figure 7.3). This cell is crossed by the nanodroplet beam. The droplets are then doped by inelastic collisions with the Mg atoms. The number of the Mg atoms in a droplet depends on the partial pressure of Mg and the cross section of this droplet. The mean number of the dopant in the droplet follows the Poisson statistics and is given by

P n = µ n n! e -µ
where n is the number of atoms picked up and µ is the mean number of n, which is proportional to pressure in the pick-up cell. The pressure dependence has been studied for many impurities (see for example [116]). The mean number of dopants, µ, can be approximated as µ = σρL where ρ is the density of Mg in the pick-up cell, L is the cell length and σ is the cross section of the droplet.

Applications of helium nanodroplets 7.4.1 Helium Nanodroplet Isolation Spectroscopy

The group of Toennies was one of the first to perform Helium Nanodroplet Isolation Spectroscopy experiments (HENDI) [118]. Then Scoles et al. [4] experimentally studied the SF 6 (monomer and dimer) embedded in helium nanodroplets. Since this time, the field has continuously progressed. Spectroscopy in He nanodroplets presents numerous advantages. The main ones can be listed as :

• Resolved rotational spectra [12][13][14][15]119].

• Controlled size of the droplets [120].

• Fast cooling of the impurities (ability to stabilize reactive species) [20,21]. HENDI appears to be a very promising technique since it uses an ideal isolation matrix.

Other applications

All the applications of the helium nanodroplets are far to have been discovered. We just list in the following the main field of applications :

• Production of very pure magnesium clusters [START_REF] Reinhard | Introduction to Cluster Dynamics[END_REF][START_REF] Diederich | [END_REF] • Chemistry of cryogenic medium [20,21,23] • Kinetics at very low temperature (0.4 K) [22] Another important application we imagine is the possibility to use helium droplets in order to separate species whose separation is infeasible by other ways. Since some species are completely solvated in the helium droplets and others are not, it may be possible to separate them using helium droplets. For example, if we have to separate a mixture of Ag and Na atoms, we can pick-up these species by a beam of helium droplets. The silver atoms will be solvated and the sodium atoms will stay at the surface. Then, after the pick-up cell, if we are interested in the silver atoms, we can make leave sodium atoms by ionizing them and redirecting them by an electric field for example.

Chapter 8 DMC computational details 8.1 Introduction

In the current section we will deal with the manner to properly make DMC calculations. As we have seen in section 6, observables obtained from the DMC algorithm should not depend, at least in theory, on the wave function parameters or the random number suite used to advance the walkers. Furthermore, observables are by definition exact in the limit of an infinite number of walkers and a time step equal to zero. But in practice these ideal limits can never be achieved. Thus, in concrete terms DMC observables such as energies or dopant positions are expected to depend on the wave function parameters. Bad results could also be produced if, for example, the imaginary time step is too large. On the other hand, an insufficient number of walkers may affect the accuracy of the calculations. Statistical errors increase when an insufficient number of blocks is used. So, the calculation strategy consists in optimizing all the technical parameters as well as possible. For this, we introduce in the following sections a systematic study of all the technical parameters. Of course, when we study a particular parameter all the others are kept constant. That is to say, if we are interested, for example, in the effect of the number of walkers, N w , all the simulations we will run to test this variable will have the same technical parameters, except for N w .

Influence of the number of walkers

It can be demonstrated that the DMC mean energy, E dmc mean , for a given number of walkers N w has a strictly positive error proportional to 1 Nw . This is true even when δt, the time step is zero. Thus E dmc mean can be written, in the limit δt = 0 as :

E dmc mean = E ∞ + C N w . (8.1)
Where E ∞ is the exact energy for N w → ∞ and C a positive constant related to the variance of the local energy and depending on the nature of the system. Thus the "exact" energy is reached when an infinite number of walkers is used. Of course, this is a physical and philosophical limit which is not realizable in practice. Nevertheless, we want to get the best possible honest estimate of the exact energy and for a reasonable accuracy, we can run calculations with a finite number of walkers. Fig. 8.1 depicts the DMC energy as a function of the number of walkers N w , points with error bars are DMC mean energies for the MgHe 12 cluster at a given ensemble size. The solid curve represents the fitting function :

E f it (N w ) = E ∞ + C N w (8.2)
where E ∞ expresses the DMC energy extrapolated to an infinite number of walkers. The general tendency is that the computed energy is well fitted by the formula above and the error bars are significatively reduced when the population size (number of walkers) is important. The DMC points of Fig. 8.1 have been fitted following the formula above (Eq. 8.1) and the asymptotic energy is found to be -18.76(1) cm -1 . The latter is in excellent agreement with that obtained using 1600 walkers (-18.75(2) cm -1 ). It is then not necessary to systematically make extrapolation to extract the DMC energy since a finite number of walkers is usually sufficient for this aim.

On the other hand, the accuracy of properties like the helium density in the droplet is also dependent on the population size. Figs 8.2 and 8.3 depict the helium density respectively for simulations with 100 and 12800 walkers. One can notice that in the first case (100 walkers), the 2d and 3d density profiles are very noisy whereas in the case of 12800 walkers, the situation is strongly different : the density profiles have become much smoother. But this has a cost, the last simulation is 130 times more expensive in time. Thus, once again, the good decision, at least in my opinion, is to adapt the number of walkers to the system size. For small clusters, we can use large number of walkers, typically 4000 to 6000, and for larger ones a number of walkers varying from 1000 to 2000. In any case, for our final DMC calculations, whatever the droplet size is, we use at least 1000 walkers to stay in a reasonable accuracy.

Influence of the time step

DMC observables are exact in the limit of a time step, δτ , which equals zero (see section 6.3.3). Bearing this in mind, we are aware that using a non zero δτ inevitably introduces systematic time step errors. So we have to find a compromise between a small time step (10-25 a.u) which make simulations long in the case of large systems (N atoms > 50) and a larger one (25-100 a.u) which is likely to cause more errors. This because, for example, the DMC mean energy, E dmc mean , can formally be expressed in a polynomial expansion of δτ :

E dmc mean = ∞ i=0 C i δτ i (8.3)
If we truncate the previous equation into the second order, we can write

E dmc mean = E δτ =0 + C 1 δτ + C 2 δτ 2 + O(δτ 2 ) (8.4)
One can observe in the previous equation, an error in energy which is proportional to δτ 2 that is why a large time step may produce more errors. One way to reasonably evaluate the DMC exact energy, E δτ =0 , is running calculations for at least 3 different time steps and then fitting the data to the previous expansion. This has been done for the MgHe 25 cluster : the curve resulting from the fitting process (see Fig. 8.4) gives an absolute energy of -49.2 ± 0.1 cm -1 as δt equals zero. Another important property is the magnesium atom radial probability density, p(R com ): Fig. 8.5 shows this quantity (p) for two distinct time steps. Circles and crosses respectively represents the helium density for a time step of 25 and 50 (a.u). It is important to underline that in both cases, the total length of simulation is the same, ie 10000 a.u. The two curves generally coincide well at large value of R com (which is the radial distance of the particles from the center of mass),ie, for R com ≥ 6 Å . Whereas at shorter distances (for R com ≈ 2 Å ), the two histograms are distinguishable. However, although the two distribution present different shapes, they are statistically equivalent since the error bars generally overlap.

Influence of the duration of the simulation

Among influencing parameters, the simulation length takes an important position. A small number of time steps may give false expectation values. In fact, considering for example the Mg@He n system, an insufficient simulation length will not let us determine correctly the position of the magnesium atom because walkers have not received enough "time" to explore the droplets and will consequently stay close to its initial coordinates.

Another main aspect of the simulation length resides in the fact that a short length gives bad energy values for the ground state of the system since statistical errors are increase by the shortness of the simulation. Figure 8.6 deals with a concrete case : the MgHe 12 cluster. It presents the variation of energy versus the number of time steps, N steps . One can notice that the energy becomes "stable" from a sufficient simulation length which is here around 100 time steps. This represents simulation time length of 2500 a.u. As a result, we can say the longer the simulation is the smaller the error bars are. Thus it is essential to choose an optimum number of time steps because a too long simulation is much more expensive and on the other hand a too short one gives noisy expectation values. So the optimum simulation length must be chosen so that the smallest overall brownian walk radius is satisfied. In this case all the particles have been given the opportunity to travel through the droplet.

The smallest overall brownian walk radius is deduced from the Einstein diffusion relation: where r is the brownian radius of the diffusing particle, D its diffusion constant and t the time. Another consequence of the shortness of the simulation can be a bad helium probability density p, as we can notice in figure 8.7 which shows two helium probability density distributions for a simulation length of 250 a.u. (square points) and 12500 a.u. (triangles). Curve presenting the short simulation shows a variation around the maximum probability density whereas the latter is expected to be smoother. This can be noticed in the second case where the distribution is smooth.

r 2 = 6Dt (8.5)

Influence of the number of blocks

Another important parameter which may strongly affect expectation values or histogram properties is the number of blocks which is used in the simulation to determine averaged quantities. Figure 8.8 shows the evolution of the DMC mean energy as a function the number of blocks in the case of a 12 helium atom cluster doped with a magnesium atom. Each point corresponds to an entire simulation. From one point to another, all parameters are identical except the number of blocks, N b . One can notice that generally the error bars decrease when N b increases. The energy tends to an asymptotic limit (-18.74(2) cm -1 ) for a number of blocks greater than 450. The horizontal line represents the energy mean value estimated from N b =450 to 1700. This shows that from 450 blocks in the simulation, all error bars overlaps. We can conclude that for such a typical system size, one has to use at least a number of blocks greater or equal to 450. Figure 8.9 presents two particle pair correlation function, π, between the magnesium atom and helium atoms. Curve with triangles shows π for a simulation of 50 blocks whereas the other curve is derived from a calculation where 1700 blocks were used. We can notice that these two curves are similar but do not overlap for all R values. The overlap can be seen for short and large R values. However, around the maximum of the distributions, these latter differ. In fact the π function in the case of a 50 blocks simulation seems to overestimate the number of helium atoms surrounding a given He atom since its maximum is higher than the 1700 blocks distribution. However, the two distributions are statistically equivalent since error bars overlap for both distributions.

Influence of the random number seed

As implied above DMC expectation values (energies, etc) (lightly) depends on the technical parameters. When changing the random number seed, fluctuations in energy can be noticed. Nevertheless, as one can notice observing figure 8.10, the energy value fluctuates around a mean value of -18.77(4) cm -1 . The fluctuation should be, however, reduced if a longer simulation is made. The fluctuation may be explained by two main factors : firstly the simulation of the current nanodroplet is made with a modest population ensemble size (2000 walkers and 200 blocks), secondly the randomness of the random numbers is not perfect. Thus we should better call them pseudo-random numbers.

Now, if we consider other properties such as the particle probability densities (here the particles in question are helium and magnesium atoms), we can say as expected that these properties also depend on the random number seed utilized for the calculations for a short simulation. If we refer to figure 8.11, one can remarkably notice a sensitive difference for the probability density curves of the magnesium atom when two different random number seeds are used. However, the two curves are statistically equivalent if one takes into account the error bars. In fact, the latter overlap for both distributions.

Trial wavefunction and parameters

Optimization of a trial wave function with many nonlinear parameters is a difficult task, especially in the field of quantum Monte Carlo simulations, where usually one uses trial wave functions to sample the system. To optimize our trial wave functions we use the variational Monte Carlo method (VMC). But in our study, we do not consider the VMC method as a final method since the DMC method is used to numerically evaluate expectation values such as energy or other properties of the system. Since VMC is directed by the variational principle, there is an optimal set of parameters a for which the energy is minimal and the trial wave function, Ψ(R, a), the closest to the exact wave function Ψ exact (R). Thus, the methodology consists in running the DMC program (DMCCOM [START_REF]The Fortran program we used to perform DMC calculations[END_REF]) in which all weights (of walkers) are set to one in running the DMC program and large time steps are used to better sample the helium droplets, changing the trial wave function parameters. This is done several times until a minimal energy is found. Furthermore others properties are verified such as the helium density or pair correlation functions. The set of parameters a which gives the minimum energy is then used in DMC calculations. The many body wave functions we have utilized is a pair product of a Jastrow type wave functions and Jastrow/Fermi ones (see example in figure 8.12). As an example of the trial wave function parameter influence, we depict in figure 8.13 the helium density in cases of optimized and non optimized wave functions for a pure helium droplet containing 12 atoms. Both curves show the helium radial density, but one can see that the density calculated with bad wave function parameters (curve with triangles) is broadly overestimated since its maximum is very close to the bulk helium density (horizontal curve).

Conclusion

Although the DMC method algorithm is quite easy to implement (at least in case of a serial code), in practice we have seen during the previous sections that several technical parameters have to be optimized by test calculations. This is because the DMC method uses random numbers but in any case we do not want our results to be random. In this way, we hope to have demonstrated that optimization parameters importance to properly perform accurate calculations. The technical parameters such as time step, simulation length among others, as well as trial wave function should be optimized for each typical system size.

Chapter 9

DMC results for MgHe n and CaHe n clusters

Introduction

While some atoms (Ag, Xe, ...) [30,124] or molecules (SF 6 , OCS, ...) [4,16] are completely solvated in the helium nanodroplet environment, other species like alkali atoms are found to reside at the surface of the droplet [25][26][27]. However, the situation of the alcaline-earth atoms is more interesting, in particular that of the Mg atom whose position with respect to the helium droplet is more ambiguous since the experiments are still in conflict. In fact, Reho et al. [28] have studied the 3 1 P 0 1 ← 3 1 S 0 transition of Mg in helium nanodroplets. They compared their spectrum with that obtained for the Mg atom in bulk helium (where the Mg atom is solvated) by Moriwaki et al. [29]. The good overlap between the two spectra (see Fig. 9.1) demonstrates, according to the authors, that the Mg atom is also solvated in the helium droplets. On the other hand, Ren et al. [30], in their recent ionization experiments, find evidence for the solvation of all the alkaline-earth atoms.

From the theoretical point of view the same contradiction persists. In fact, DMC calculations performed by Mella et al. [31] show solvation of Mg only for large clusters whereas DFT calculations produced by Hernando et al. [32] clearly identify a central position of the Mg atom.

In the following, we present the main results from our DMC calculations for both calcium and magnesium atoms. External positions are found for both calcium and magnesium. The Mg position is not influenced by replacing the helium pair potential of Aziz et al. [34] by that recently published by Jeziorska et al. [35].

Ancilotto's model 9.2.1 Principle

A simple and approximate way to estimate whether an impurity is solvated or not is given within the Ancilotto model [24]. The essential idea of this model is based on the fact that the gain in interaction energy due to dopant-helium interactions is balanced with the cost in energy needed to produce a cavity inside the helium droplet for the dopant. The balance is quantified by a dimensionless parameter λ. The latter allows the determination of the surface location or its complete solvation. The adimensional parameter λ based on the well depth ǫ (in cm -1 ) and the corresponding equilibrium distance r e (in Å ) of the interaction potential between the impurity and the helium atom is defined as :

λ = ρǫr e /(2 1 6 σ)
where ρ = 0.022 Å -3 and σ = 0.179 cm -1 Å -2 are respectively the density and surface tension of bulk liquid helium. The threshold for solvation in 4 He is [24] λ ≈ λ 0 = 1, 90. If λ < 1.90, an energetically favoured surface position of the impurity is expected, otherwise (if λ > 1.90) the impurity solvation is energetically favoured. Species with small λ(<< 1.90) values, like alkali atoms, are effectively on the surface of the droplets whereas atoms of silver which presents a large (compared to 1.90) λ value (5) are solvated.

Limits of the model

Since the model is based on a Lennard-Jones type interaction potential, the shape of the real interaction potential is not taken into account. This is not a real problem for species characterized by λ values far from the threshold value (1.90). However, the shape of the binary potential between impurities and He remains important when the λ value are close to the threshold value. This situation is encountered when the impurity are Mg or Ca atoms. The model uses the density and the surface tension of bulk liquid helium to determine the λ parameter but we know that the density of helium in the droplets is not uniform and the surface tension of the latter may depend on their sizes. All these approximations make the Ancilotto model applicable to delicate cases like alkaline-earth doped helium droplets.

The alkaline earth case

Alkaline earth atoms are very interesting system because of their extremely weak interaction with helium atoms. Because of this weakness, the λ parameter of Ancilotto is close to These values, as one can notice, are very close to the critical value. This is true for all the ab initio methods used. Thus, we can not use the Ancilotto model in the specific case of Mg@He n and Ca@He n 9.3 Pair potential of the He 2 , MgHe and CaHe

In chapter 9, we have shown that the DMC method needs interaction potentials for each type of particle pairs. In principle, the pair potential should take into account two-body terms, three-body terms and terms of higher order. Since it has been demonstrated that, in the case of the interaction between helium atoms, the three body-terms are negligible versus the two-body ones [125,126], we have only constructed our potential for the helium atoms as a sum over two-body interaction :

V He = i>j v(r ij ) (9.1)
where v(r ij ) is the analytical pair potential of Aziz et al. [34] or the more recent one of Jeziorska and co-workers [35]. and in the same way the contribution due to the XHe (X = Ca or Mg) complex is written as

V M g = i v CCSDT (r i ) (9.2) V Ca = i v CCSD(T ) (r i ) (9.3)
where v CCSDT (r ij ) is the pair potential calculated at the CCSDT level of theory and is mentioned in chapter 5. In Eqs. 9.1 and 9.3, r ij is the distance between particle i and particle j and r i the distance between particle X and particle i. The total CaHe potential energy is determined following Eq. 9.3 and using the CCSD(T) pair potential reported in chapter 4. Within this approach the interaction energy in the ground-state clusters is predicted simply by summing up all the pair wise HeHe (= V he ) and XHe (= V X ).

Energy calculation

For a doped helium droplet with an X atom (X = Mg or Ca), we define the binding energy as the energy difference between the doped droplet and the pure helium droplet

E bind = E XHe N -E He N (9.4)
Using Eq. 9.4, we produce the binding energies reported in table 9.2 which are plotted in Fig. 9.2. These energies are obtained considering the Jeziorska potential for the helium dimer. Fig. 9.3 represents the absolute energies of the helium clusters (circles) and the absolute energies of the same clusters doped with a magnesium atom(triangles). One can notice that both curves have the same general monotonic shape and present a negative curvature for small cluster sizes (for N < 25). For clusters with N up to around 25 helium atoms, the absolute energies become quasi-linear versus N. This phenomenon can be explained in both cases (doped and pure clusters) by a "saturation effect". This means, in case of doped clusters, that from a certain size the magnesium atom does not feel the presence of more helium atoms. By a similar assumption, it is possible to explain the shape of the curve drawing the absolute energies of pure helium clusters. In fact each helium atom will feel the interactions of the first shells because the potential interaction rapidly decreases with the distance between two particles. That is why, for large pure helium clusters, the absolute energy is almost linear in N. The same conclusion can be done for the energy point presented in Fig. 9.4 where the points are obtained using the helium pair potential of Aziz al. [34]. Absolute and binding energies, in the latter case, are given in Table 9.3. The binding energies, obtained according to the latter helium binary potential, are reported in Fig. 9.5 and the binding energy limit is close to that in Fig. 9.2 (about -20.25 cm -1 ).

On the other hand, in Fig 9 .6, for the Ca doped clusters, the binding energies (b) and absolute DMC mean energies (a) are plotted. The pair potential used for helium atoms is from Ref. [34] and that between Ca and He is the CCSD(T) pair potential calculated in Chapter 4. The binding energy of the CaHe N (when N → ∞) is around -16.0(3) cm -1 . Binding energies of the Mg@He N complexes when using helium pair potential of [35]. Solid line is the fitting function from Eq. 9.17. 9.2: DMC mean energies for pure and corresponding doped helium clusters. The last two columns give the binding energies for a given N. All energies and standard deviations are expressed in cm -1 . The helium pair potential of Jeziorska [35] was used in combination with the CCSDT potential for Mg-helium interaction. .5: Binding energies of the Mg@He N complexes when using helium pair potential of [127]. The last two columns give the binding energies for a given n. All energies and standard deviations are expressed in cm -1 . The data are obtained using the helium pair potential of Ref. [34] in combination with the CCSDT potential for the MgHe pair interaction. The last two columns give the binding energies for a given n. All energies and standard deviations are expressed in cm -1 and the time step, δt in atomic units. The helium pair potential of Ref. [34] have been used in the DMC calculations. The Ca-helium pair interaction was described by the CaHe CCSD(T) potential.

n δt E(M gHe n ) δE(M gHe n ) E(He n ) δE(He n ) E n δE n 4
Since no model is presented in the literature (at our best knowledge) for fitting the binding energies, we have developed a model in next section. This model is based on an analytical expression of the pair potentials.

Binding energy model

In the aim to determine an analytical fitting function for the binding energies, we have proposed a new method which is based on pair interactions between particles. We suppose that the interaction between atom i and atom j is given by The total potential energy for a pure helium cluster is then described by

V ij (r ij ) = A exp(-βr ij ) - C 6 r 6 ij ( 9
V He N = N i>j V ij = V 1 (9.6)
whereas for an X doped helium cluster Eq. 9.6 becomes

V XHe N = N i>j V ij + N i V iX (r iX ) = V ′ 1 + N i V iX (r iX ) (9.7)
and the total absolute energies of the pure helium cluster is written

E He N = V He N + T He N (9.8)
where the new term in Eq. 9.8, T He N , represents the total quantum kinetic energy of the pure helium cluster. Whereas in the X doped helium clusters, we obtain the following formula

E XHe N = V XHe N + T XHe N (9.9)
where T XHe N is the total quantum kinetic energy of the X doped helium cluster. If we remember that the binding energy is defined as

E bind. (N) = E XHe N -E He N (9.10)
Taking into account that the difference between T XHe N and T He N we can be supposed constant, Eq. 9.10 can be rewritten as

E bind. (N) ≈ V XHe N -V He N (9.11) E bind. (N) = V ′ 1 + N i V iX (r iX ) -V 1 (9.12)
Now, if we assume that the impurity does not really "disturb" the helium environment, especially true for low interacting surface located impurities or large size system, we can then suggest that V ′ 1 ≈ V 1 so that equation 9.12 becomes

E bind. (N) = N i=1 A exp(-βr iX ) - C 6 r 6 iX (9.13)
Transforming the sums in equation 9.13 into integrals the following relation is obtained :

E bind. (N) = rmax r 0
A exp(-βy) -C 6 y 6 y 2 dy (9.14)

E bind. (N) = N 1/3 N 0 A exp(-βt) - C 6 t 6 t 2 dt (9.15)
Finally after integration of equation 9.15, the binding energies could be fitted to the following formula which has a physical sense :

E(N) = -a 1 ( N 2/3 a 2 + 2N 1/3 a 2 2 + 2 a 3 2 ) exp(-a 2 N 1/3 ) + a 3 C 6 N + E ∞ (9.16)
In the previous relation the a i are positive real parameters C 6 is the the dispersion coefficient and E ∞ the extrapolated binding energy. An equivalent form which is more convenient for the fitting process is given below:

E(N) = (b 1 N 2/3 + b 2 N 1/3 + b 3 ) exp(-b 4 N 1/3 ) + b 5 C 6 N + E ∞ (9.17)
where b i are real parameters and E ∞ has the same meaning as previously discussed (in Eq. 9.16). The fitting model is also based on the approximations that the helium density is constant and the droplets are spherical. Obviously, this simple model can be improved by explicitly introducing a modelling function for the real density of the droplets. It can also be generalized to the many-doped droplet case.

In order to verify the validity of our model (Eq. 9.17) we have summarized in table 9.5 the difference between the DMC data and the corresponding energies derived from Eq. 9.17. One can hence remark a very good agreement between DMC values and the fitting model. Following the fitting function presented in Eq. 9.17, we obtain an extrapolated binding energy (for N → ∞ ) of -20.9 ± 0.8 cm -1 . This value is already reached for MgHe 100 whose binding energy is -21.51 cm -1 (±0.56).

Comparison with literature

The structural problem of Mg doped helium clusters has been treated in the literature. In Ref. [31], a DMC calculation has been done in order to determine the magnesium location in the droplet. The authors (Mella et al) found that the Mg atom is fully solvated in case of helium clusters composed whith more than 25 helium atoms whereas for smaller helium clusters an opposite situation is predicted for the Mg atom. A transitional situation was found for the MgHe 25 cluster (Fig. 9.7). Performing these calculations the authors have used as a pair wise interaction potential between He and Mg the CCSDT pair potential determined by Hinde [77]. What can be surprising in the results given by Mella et al. is that the position of the Mg atoms depends on the cluster size. Secondly, the DMC method used by Mella seems to fail to produce bound states for very small Mg doped helium clusters. Table 9.5: Comparison between the DMC binding energies (E n ) and the extrapolated energies (E f it ) from the fitting process according to Eq. 9.17. δE n and ∆E n are respectively DMC error bars and absolute errors of the fitting function with respect to the DMC points. Energy units are cm -1 . Another approach has been applied to the Mg doped helium cluster. Hernando et al [32] have used the DFT method in combination with the MgHe pair potential of [77]. They have studied clusters whose sizes were between 300 and 5000 helium atoms. They found a completely solvated Mg atom whatever the cluster sizes are. Fig. 9.8(b), from [32], depictes the helium density for clusters with size from N=300 to N=5000 and figure 9.8(a) represents the helium probability density for MgHe 220 which we computed. Comparing these two figures we notice that the helium density profile computed by means of the DFT method in Ref. [32] shows a lack of helium density at the origin of the diagram, where our helium density profile (Fig. 9.8(a)) does not. Hence the DFT method places the Mg atom at the center of the droplet whereas we find a near surface location for this atom. Since the DFT method is a semi-empirical method, its results may strongly depend on the choice of the functional. Hence, the DFT method does not appear to be the method of choice to determine the location for such a subtle impurity. We are then not very confident in the results given by this method. We are more confident in the DMC method since it solves the Schrödinger without significant approximations (see Chapter 6 for further details).

In the experimental way, Ren et al. [30] have measured the electron energy dependence of the ionization yield of alkaline-earth and xenon picked-up atoms. They found a qualitative shape difference between the yield curves of species solvated in the middle of the droplet (Xe atom) and species located in the surface region (alkaline-earth atoms). Their measurements, which can be seen in figure 9.9 demonstrate that all the alkaline-earth atoms reside at or near the droplet surface. These experimental results all agree with our theoretical results in the sense that a surface location is found for both Ca and Mg for any helium pair potential.

Ca and Mg positions on the droplets

Except the DMC energy of the doped helium clusters, an important expectation value to consider is the structure of the system that is to say the position of the impurity with respect to the center of mass of the droplet. This can be determined by observing the radial probability density of the alkaline earth atoms or the helium density in cylinder coordinates. The system of coordinates is defined so that the origin coincides with the center of mass of the system and the impurity with the positive Z-axis.

Radial probability densities

Figures 9.10 and 9.11 respectively depict the radial probability densities of Ca and Mg atoms for several cluster sizes. One can notice independently from the cluster size, for both Calcium and magnesium atoms, that the radial probability densities of the dopant take maximum values far from the center of the droplets. This is a proof of a surface location for both Mg and Ca.

Helium densities in cylinder coordinates

Now if we orient our attention into the helium density for several sizes of doped helium clusters, one can clearly notice that Mg resides at the surface of the droplets considering Fig. 9.12. A lack of helium density can be observed for external positions which imply the presence of the impurity (Mg) at these positions. The contour plots of the helium densities (figures 9.15 and 9.13) better shows the lack of helium density and these figures Radial probability densities of Mg in helium clusters whose size are respectively, from the left to the right of the figure, 18, 25, 32, 53, 83 and 170. These distributions are obtained with the helium pair potential from Ref. [35].

show that both Mg (Fig. 9.13) and Ca (Fig. 9.15) reside at the surface of the droplets. In the interesting case of the Mg atom, the helium pair potential of Jeziorska et al. [35] or that of Aziz et al. [34] does not influence the position of the foreign atom. In this way, Fig. 9.14 also demonstrates an external position for Mg when the Aziz potential is used for helium in the DMC simulations in case of MgHe 20 and MgHe 100 . Fig. 9.12 shows the three-dimensional helium density for several MgHe N clusters and Fig. 9.13 the corresponding projections where the Mg position is manifested by the lack of helium density on the bottom left corner of each subfigure. Fig. 9.13 clearly supports a surface location of the Mg atom whatever the cluster size is.

Structural relaxation of the MgHe N cluster

A concrete view of the MgHe N clusters can be imagined. Let us consider the specific cases of a small and a relatively large cluster (MgHe 32 and MgHe 170 ). Fig. 9.16 shows a example of the evolution of the MgHe 32 in the imaginary time τ . The DMC calculation has started with an initial configuration, optimized at the classical dynamic level, with Mg near the center of mass. One can qualitatively notice on the bottom left corner that the cluster presents a very organized structure which is not foreseen for a quantum system whereas the typical DMC state shows a very disordered helium cluster where Mg has moved to the surface. This difference between classical dynamic structure optimization and the DMC resulting state is noteworthy and can be seen for the MgHe 170 droplet in Fig. 9.17.

Pair density distributions

Figure 9.18 shows the pair density distributions of helium atoms in case of doped helium cluster with different sizes. The maximum density corresponds to the nearest neighbour distance between 2 helium atoms. In figure 9.18, all curves show a coincidence for the first maximum which is situated at a distance of 3.6 Å . However the distance corresponding to the maximum of the distributions decreases when the size of the helium clusters increases. This phenomenon can be explained by the fact that the density of small helium cluster is below that of bulk helium. The bigger the clusters are the closer to the bulk helium density are their densities. In the other hand, figure 9.19 presents the same type of curves than previously, for magnesium doped helium cluster. In figure 9.19 the maxima of the Mg-helium pair particle distributions are all situated at a distance of ≈ 5.4 Å . The evolution of the curve maxima can be explained as just discussed above. In order to determine the energy profile of the MgHe N system, a "pseudo-shake" algorithm has been used in the DMC code. This gives us the possibility of fixing the distance between an atomic impurity and the center of mass of the helium atoms. The only constraint concerns this distance and the impurity has, however, the ability to make free rotations.

Applying this constraint, we have performed DMC calculations for four typical systems: MgHe 50 , CaHe 50 , NaHe 50 and AgHe 50 . For all these systems, the DMC mean energy has been calculated in a range started from 0 to up to 20 a 0 , steps of 0.5 a 0 of the constraint. Figures 9.20 and 9.21 depict the DMC mean energy as a function of the length of constraint for the four systems. As expected, Fig. 9.20, shows that the energy of the AgHe 50 system (curve with circles) increases with respect to the constraint distance (R) and the minimum is found for a distance of zero. Recalling the definition of the constraint, the silver atom is then energetically favourable to a complete solvation. This fact is in excellent agreement with both theoretical [124] and experimental [128] results. The same conclusions are valid for the calcium atom since in Fig. 9.21 the energy profile is minimum for a distance corresponding to a surface location of this dopant (curves with bold circles).

In the other hand, the sodium dopant shows a completely opposite energy profile compared to the the previous cases (Ag and Ca) : the DMC energy (curve with circles) generally decreases with respect to the constraint length (Fig. 9.21). Sodium atom is hence energetically favoured by a surface location. This is experimentally well known [24,25]. Now, if we consider the case of the MgHe 50 system, observing figure 9.20 (bold circles), it presents a noticeable aspect in the sense that the energy curve is very flat in the interesting domain. This demonstrates that the magnesium atom is energetically indifferent about its position in the droplet. However, observing figure 9.20, the energy curve of MgHe 50 shows "some structure". In fact in the area of small constraint distances the curve presents local minima and maxima. But this structure may not be interpreted as a real one if statistical errors are taken into account. Moreover, applying the constraint algorithm between the center of mass of the helium atoms and the impurity does not lead to a clear dissociation balance because the dopant leaves the droplet with generally some rare gas atoms, depending on the pair interaction between the dopant and the helium atoms. The balance for a large presents some strange aspect in the sense that when progressively moving Mg away from the center of the droplet by varying the constraint, this density shows irregularities. This can demonstrate that the droplet does not keep a spherical shape and Mg moves keeping with it some helium atoms. The same observations have been made in the MgHe 50 cluster. Fig. 9.23 and 9.24 depict the 2d-projection of the helium density when MgHe 20 and MgHe 50 are respectively considered. One can clearly notice the evolution of the Mg atom in the droplet when the constraint distance increases. This shows that the constraint algorithm works well.

Rovibrational calculation in the constrained potential

In this section, we performed rovibrational calculation using the potential descended from the energy scanning to determine the fundamental rovibrational level of the MgHe 50 complex. For this purpose the radial Schrödinger equation is solved using a variational program.The radial Hamiltonian is written as :

Ĥ = - h2 2µ ∂ 2 ∂r 2 + V (r)
where h is the reduced Planck constant and V (r) the potential energy operator formed by the DMC energy of the MgHe 50 scanned over the distance between the center of the helium atoms and the single magnesium atom. And µ is the reduced mass of the system defined as :

µ = m M g 50 i=1 m He m M g + 50
i=1 m He The radial Schrödinger equation is variationally solved in a set of 200 Laguerre basis functions. A rovibrational level whose energy is -119.74 cm -1 is found. The average distance between the helium group and the Mg atom is 4.48 Å . According to these results the Mg atom is situated at an intermediate distance between the surface and the center of the droplets since the classical radius of this droplet is about 8 Å . For comparison, a DMC calculation (without constraint) of the same system gives an energy of -117.0 ± 0.1 cm -1 which is sensibly different from that previously found for the rovibrational state. As expected the Mg atom is located near the surface as one can notice observing Fig. 9.25.

However, in the rovibrational calculations the problem is how to choose the reduced mass of the system since the separation between the helium group and the Mg atom is not clear as previously discussed. In fact when Mg moves in the droplet, some helium atoms accompany its motion. Hence there is no simple way to determine the reduced mass. And the latter can change with the Mg motion.

Conclusion

DMC calculations based on very accurate ab initio, ie, CCSDT potential for the Mg-He pair interactions, CCSD(T) potential for the Ca-He pair potential, in combination with two accurate helium pair potentials [34,35], have demonstrated a surface location for both Mg and Ca atoms. In particular the position of the Mg atom, as well as its binding energy to the helium droplet, is very robust with respect to the choice of a particular helium potential. Although the recent Jeziorska helium pair potential [35] is deeper than the Aziz potential [127], the latter is sufficient for studying such systems and its main advantage is that the calculations become less time consuming.

In comparison with Mella et al. who have performed calculations of the MgHe N clusters with N up to 50 and found a surface location of Mg for small clusters and solvated state for larger one, our calculations of the same system (with N up to 220) always give a surface position of Mg whatever the cluster size is. Furthermore our predictions for both Mg and Ca are in good agreement with the recent ionization experiments of Ren et al..

From an energetically point of view, our constrained DMC methods has shown a very flat potential for both MgHe 20 and MgHe 50 . The Mg atom seems to have no energy preference when travelling in the droplets. The binding energies are about -16 and -21 cm -1 for respectively calcium and magnesium doped clusters.

However, the Mg@He N system is far from being completely understood since for example, when helium nanodroplets are doped with more than one Mg, the presence of single Mg atomic transition is detected [115].

Chapter 10

Dynamics of Mg doped Helium Clusters 10.1 Introduction

The main goal of this dynamics study of Mg doped helium clusters (with one and two Mg atoms) is an attempt to understand the experiment performed by Przystawik et al [115] which shows evidence of single atoms in multiply Mg doped helium droplets. The key idea of our dynamic simulations is to know whether two Mg atoms rapidly collapse or slowly diffuse in a droplet composed of around 2000 helium atoms. Dynamic calculations based on the ZPAD (Zero Point Averaged Dynamics, explained in Ref. [129]) have been carried out and methods and main results are detailed in the following.

Potential energy curves

10.2.1 Mg 2 (X 1 Σ + g )
For the dynamics simulations of helium clusters doped with more than one magnesium atom in their ground states, the Mg 2 pair potential is needed. One of the most accurate Mg 2 ground state potential available in the literature is the RKR one [130]. Unfortunately, the RKR potential is restricted to relatively short interatomic distances. For this reason, we had to fit the RKR data into an HFD-B type function, V (r), where dispersion forces have only been taken into account by the C 6 coefficient : With the model of Eq. 10.1, we have obtained long range points which allow us to perform dynamic simulations of multiply Mg doped helium clusters. The fitting function is in relatively good agreement with the RKR points as can be seen in Fig. 10.1. The RKR Mg 2 potential has an equilibrium geometry equal to 3.89 Å and a well depth of 424.00 cm -1 whereas our value from the fitting functions are lightly different, however in reasonable agreement (r e =3.88 Å and D e = 427.07 cm -1 ). This difference is not really important since the r e and D e values present larger uncertainties in the literature [132][133][134].

V (r) = A exp(-βr) - C 6 r 6 (10.1)

MgHe (X 1 Σ + )

The ground state MgHe pair potential used in the dynamic simulations is that calculated and fully detailed in chapter 5. We just recall that this potential energy curve is obtained at the coupled clusters level of theory with explicit treatment of double and triple electronic excitations (CCSDT). The CCSDT method has been used in combination with large basis sets for both helium (aug-cc-pV5Z) and magnesium (aug-cc-pCVQZ) with additional set of bond functions (BF). According to the heavy ab initio treatment, the X 1 Σ + MgHe state is depicted in Fig. 10.2. This potential will be used in the dynamic simulation.

Effective He 2 potential

Helium clusters are known to be dominated by quantum effects which influence the binding energy per particle as well as the helium density. For example the ground state binding energy of He 2 contains more than 99.99% of the well depth as zero point energy (ZPE).

Because of these quantum effects, dopant embedding dynamics will differ from classical dynamics. In the present section, we describe an approximate technique which attempts to include part of the quantum effects for the helium clusters through an effective potential for the helium pair interaction. Basically, each particle is represented by a single particle wave function. Particle positions evolve according to classical dynamics on the effective potential. This iterative approach is described in Ref. [129] and is basically composed of 5 steps. Conceptually, each particle is replaced by a probability distribution Φ 2 (r) centered on its classical position, and assumed to have spherical symmetry. Starting from the original pair potential V cl (R) = V q,0 (R), where R = |R| represents the interatomic distance, a delta distribution Φ 2 0 (r), and known masses, the construction is based on the following sequence of calculations, which is repeated until the nth-order quantum effective potential V q,n and all distributions have reached self-consistency:

1. construction of the pair correlation function P n (R) between classical particle positions from a classical molecular dynamics simulation at temperature T with the current potential V q;n (R).

2. convolution of the pair correlation function with the current single particle distribution Φ 2 n (s) according to

P q;n (|R|) = P n (|R|)Φ 2 n (|R -R ′ |)dR ′ (10.2)
3. construction of the radial potential V rad;n (r) experienced by each particle in the "cage" formed by the other particles by integration over the pair potential and the particle distribution following V rad;n (|r|) = V cl (|r -R|)P q;n (|R|)dR (10.3) 4. solution of the radial Schrödinger equation in the radial potential V rad;n (r) in order to find Φ n+1 (r) for each particle in the mean field of the others.

-h 2 2m 4 H e d 2 dr 2 + V rad;n (r) Φ n+1 (r) = E n+1 Φ n+1 (r) (10.4)
5. construction of the next generation effective pair potential V q;n+1 (R) by the convolution

V rad;n+1 (|R|) = V cl (|R + r -r ′ |)Φ 2 n+1 (|r|)Φ 2 n+1 (|r ′ |)drdr ′ (10.5)
Following the previous five step algorithm an effective potential for the He-He interaction has been constructed. It is also possible to construct an effective potential for the MgHe interaction in a similar fashion but we have judged that it was not primordial since the ZPE of this system is not as extreme as that of helium. Interpolation by cubic splines is used for the pair potentials computed on a grid for the dynamic simulations (step 1) and for the solution of the Schrödinger equation (step 4). Five or six iterations are required to reach convergence. The convergence of the effective helium pair potential can be seen in Fig. 10.3. In this figure, the deepest curve is the classical potential. After five iterations, the effective potential seems to approach convergence. Observing Fig. 10.3 one notices that the well depth corresponding to the converged effective potential (-1.55 cm -1 ) is much shallower than that of the classical potential (-7.6 cm -1 ). The classical equilibrium distance (2.96 Å ) is shifted to a larger value of 4.19 Å . These changes in well depth and equilibrium distance lead to an approximately correct amount of zero point energy in the system and a correct bulk helium density (0.022 Å -3 ).

Dynamic results

The classical trajectories, based on the effective potential previously discussed, have been integrated using a velocity Verlet algorithm [135]. For both MgHe 1998 and Mg 2 He 1997 , simulations have been carried out using a time step of 20 fs which enough accurate since the relative energy fluctuation is only 4.6 × 10 -7 for both systems. Four hundreds consecutive simulations have been made, each of them has a length of 200 ps. Each simulation waits 2000 fs before starting the computation of histograms and expectation values. The total simulation length, when cumulating all the sub-simulations, is 80 ns. This duration is sufficient to allow to the Mg atoms a brownian radius comparable to the droplet diameter of 32 Å . In the case of the MgHe 1998 system, the initial position of the Mg atom has been set near the center of the droplets whereas for the Mg 2 He 1997 cluster the two Mg atoms were solvated in the helium droplet and the initial distance between them was 28.7 Å . The latter value is sufficiently large so that the two Mg atoms do not immediately recombine since Przystawik et al. predict a metastable state of the Mg 2 dopant dimer with a equilibrium length of 10 Å .

MgHe 1998

Thanks to the Zero Point Averaged Dynamics, abbreviated as ZPAD (classical dynamics based on an effective potential), we have a quasi-stationary position for the magnesium atom in the droplet. In fact, in Fig. 10.4(a) is drawn the evolution of Mg projected onto the (XY) plane. The initial coordinates of Mg have been set to the origin of the graph. Each point in Fig. 10.4(a) is the result of a 200 ps simulation. It can be noticed in the latter figure that the Mg atom has a brownian diffusive motion but it stabilizes in an area near the center as if it was caged. These effects are also visible in Fig. 10.4(b) which shows the MSD (mean square displacement) of Mg (curve with oscillations). The latter curve shows a stationary position of Mg since the MSD does not grow with respect to time but only presents some oscillations. Whereas the MSD of He is a monotonic growing function with respect to time; this is a proof of a diffusive motion of helium atoms and equally demonstrates the liquid character of the droplet at the temperature of 0.38 K. The liquid character of the droplet at T=0.38 K is also clearly visible when observing Fig. 10.5 which depicts the radial pair particle distribution of helium. In fact, the P rad curve only shows three main peaks which correspond to the three shell neighbours. This is a typical feature of a disordered system, hence a liquid phase.

In order to understand how Mg moves inside the helium droplets, we have plotted in Fig. 10. 6(a) the distance between the initial Mg position and its position at time t. The results obtained in this figure show a maximum distance of Mg of about 3 Å from its initial position at a simulation length of 50 ns. Beyond this time, the distance starts decreasing. This tendency is confirmed observing Fig. 10.6(b) which represents the Mg radial density for a long simulation of 80 ns (curve with triangles) and a short one (0.2 ns). The two distributions almost overlap. In brief, Mg seems to be not very mobile in this medium. This fact is quantified by the weak average diffusion constant of Mg which is 7.5(2) × 10 -11 m 2 s -1 . The self-diffusion constant of helium averaged over four hundred simulations is found to be 4.1(7) × 10 -9 m 2 s -1 . This value is in excellent agreement which those reported in the literature : 5.7-5.8×10 -9 m 2 s -1 in Ref. [136] and 4.5-5.8×10 -9 m 2 s -1 in Ref. [137] at a temperature of 4 K. The difference between our result and those reported in the literature is probably due to the temperature at which the diffusion constant are calculated. In our case self-diffusion constant of helium are extracted from a 0.38 K simulation. The diffusion constant of Mg is equal to 7.5 × 10 -11 m 2 s -1 at the same temperature. The weak constant can explain the slow diffusion of Mg in helium droplets. On the other hand, as announced in the introduction, our effective potential fairly well reproduces the helium density as can be seen in Fig. 10.7. One can see that our computed density is fairly close to the bulk one which is 0.022 Å 3 .

Mg 2 He 1997

Concerning the interesting case of the Mg 2 He 1997 cluster, the dynamics simulations whose total length is 80 ns give a very surprising result. In fact, despite their relatively strong interaction, the two Mg atoms seem to not recombine and conserve the initial distance between them as can be noticed in Fig. 10.8(a). In the latter, the distance between the two Mg atoms is drawn as a function of the simulation length. This distance just oscillates around a value of 29 Å close to the initial separation distance of 28.7 Å . This fact is better shown in Fig. 10.8(b) which is a zoom of the former one. The quasi-constant distance between the two Mg atoms is also shown in Fig. 10.9(a). One can notice a stationary distance between the Mg atoms. The apparent rotation of Mg 2 is in reality due to the rotation of the whole system which has a non zero angular momentum. Fig. 10.9(b) is a zoom of the motion of the Mg atom located near the center of the droplet. It shows a similar feature to Fig. 10.4(a). The Mg atoms when sufficiently far from each other seem to independently diffuse and do not recombine after a simulation length of 80 ns. This can be an interpretation of evidence of the single Mg presence in multiply Mg doped helium droplets found in the literature [115].

Conclusion

Along this chapter, the dynamic results based on an effective helium pair potential have been presented for two typical systems : MgHe 1998 and Mg 2 He 1997 , typical because these systems are large enough to allow us comparisons with experiments. Classical trajectories have been integrated using the velocity Verlet algorithm [135] with a time step of 20 fs which leads to relative energy fluctuation of only about 4×10 -7 . The RKR potential [130] has been extrapolated to large interatomic distances via a HFD-B form and has been used as the Mg pair interaction. For the MgHe interaction, our calculated CCSDT potential has been used in the simulations in combination with the Jeziorska pair potential for helium.

As results from the dynamic simulations, the self-diffusion constant of He, at 0.38 K, has been found to 4.1(7) × 10 -9 m 2 s-1 in good agreement with the literature [136,137]. For Mg, the diffusion constant has been found equal to 7.5 × 10 -11 m 2 s -1 at the same temperature but no data are available in the literature at our best knowledge. Finally this dynamics study the Mg 2 He 1997 cluster shows a diffusive brownian motion for the two Mg atoms. The Mg atoms do not recombine after 80 ns, this fact can explain the experiment of Przystawik et al [115] which shows evidence of single atoms in multiply Mg doped helium clusters.

Chapter 11

General conclusions

In this thesis, we have presented a deep theoretical study of calcium and magnesium helium doped clusters. In particular the MgHe n systems with n up to 220 helium atoms have been highlighted. The study has been performed using the DMC method in combination with accurate interaction potentials. Since the latter are known to be of extremely weak van der Waals interaction type, we have carried out intensive ab initio calculations at high level of theories such as MP2, MP4, CCSD(T) and CCSDT in combination with large and diffuse basis sets and a large bond function set.

Concerning the Mg-helium ground state (X 1 Σ + ) pair interactions, the well depth and corresponding equilibrium distance are respectively evaluated to -5.06 cm -1 and 5.07 Å , using the highest level of theory (CCSDT). These results nicely agree with the recent theoretical predictions [77] and are also consistent since CCSD(T) gives -4.77 cm -1 and 5.09 Å , MP2 -4.38 cm -1 and 5.16 Å , and MP4 -5.70 cm -1 and 5.03 Å , respectively for the depth well and the equilibrium distance. On the other hand, for the CaHe ground state PES, we have obtained from the CCSD(T) methods a well depth of -3.31 cm -1 and an equilibrium bond length equals to 5.83 Å , results which are in good agreement with those presented in the literature [66,68,77]. Moreover, for both CaHe and MgHe ground states, accurate dispersion coefficients (C 6 , C 8 and C 10 ) have been obtained. The latter are consistent for all the ab initio approaches used. Furthermore, for the two complexes (CaHe and MgHe) the dispersion coefficients are in good agreement with both theoretical predictions [66,68,77] as well as with experimental values [81]. According to the previously discussed informations (C i , ǫ and r e ) we think to have determined accurate PES for MgHe and CaHe ground states.

The highest level PES therefore calculated have been used in the DMC method as a sum over pair interactions in combination with accurate analytical helium pair poten-tials [34,35]. For both, helium potentials and the CCSDT ground state pair potential of MgHe, DMC calculations systematically place the Mg impurity near the surface of the helium droplet whatever the droplet size is. The binding energies of the MgHe n systems (with n from one to up to 220) have been computed as a function of n. For both helium potentials used, the binding energies converge to a value of around -21 cm -1 . The latter is only reached from a hundred atom helium droplet. In order to better predict the asymptotic binding energies (n → ∞), a fitting function model based on HFD pair interactions has been developped. The latter seems to be promising in extrapolating binding energies (from helium droplets to bulk helium).

Concerning the DMC results of the CaHe n clusters, the ground state CaHe CCSD(T) pair potential has been used to carry out DMC calculations with n up to 190. The results obtained by using the CCSD(T) CaHe PES in the DMC calculation have been reported and provide a surface location of the calcium atom and a binding energy equals to around -16 cm -1 .

Although our DMC results for both Ca and Mg doped clusters are in good agreement with the recent experiment of Ren [30], we have computed the MgHe n energy as a function of the Mg position in the droplets in order to understand the conflicting structural results in the published literature. In order to determine the energetically favoured position of Mg in the helium clusters, we have used a constrained DMC algorithm. Energy calculations via the latter have been carried out scanning the distance between the impurity and the helium center of mass. Four typical systems (AgHe 50 , CaHe 50 , Na 50 and MgHe 50 ) have been treated by the constrained DMC method. All systems present a discriminating energy profile, except the MgHe 50 cluster which gives a very flat energy profile. This, therefore, means that the Mg atoms does not seem to have energetically preferred position in the droplets.

Finally, in order to compare with experiments, classical dynamics trajectories based on an effective potential have been computed for MgHe 1998 and Mg 2 He 1997 . The effective potential have well reproduced the helium density as well as the droplet energy. The Mg 2 He 1997 study has been motivated by the fact to know whether the two Mg atoms rapidly meet thanks to their relatively strong interaction (with a well depth of -424 cm -1 ) or they adopt a purely brownian moves before meeting. The results we get is that they do not rapidly meet. This can be a reasonable explanation for both experiments and theoretical calculations which find a metastable state (with a large equilibrium distances) of Mg 2 in helium nanodroplets. This work is a first step of a wide project and can be completed by studying dynamic and structural properties of multiply Mg (and Ca) doped helium clusters. It is also interesting to perform constrained DMC simulations with more than an alkaline-earth atom. It also appears interesting to investigate mixed alkaline-earth reactions in helium droplets by mean of our ZPAD method. -are counter-poise corrected with respect to the (33211) bond functions. The aug-cc-pV5Z and aug-cc-pCVQZ basis sets were respectively used for He and Mg in combination with the bond functions (33211).

Appendix A Electronic energies

Table A.2: CCSD(T) energies for the MgHe complex and its components. The binding energies are counter-poise corrected with respect to the (33211) bond functions. The Mg basis was aug-cc-pCVQZ+BFs and the He one was aug-cc-pV5Z+BFs. A.3: CCSDT energies for the MgHe complex and its components. The binding energies are counter-poise corrected with respect to the bond functions (BFs=33211). The Mg basis was aug-cc-pCVQZ+BFs and that of He was aug-cc-pV5Z+BFs. A.9: CCSD(T) energies for the MgHe complex. Cp means that energies are counter-poise corrected. All energies are expressed in micro-hartree. The Mg basis was aug-cc-pCVQZ and the He one was aug-cc-pV5Z (see Fig. Radial probability density of Mg in the MgHe 53 cluster for different interaction potentials. The reference potential is the CCSDT one. Increasing the CCSDT potential by only 10% is sufficient to give a fully solvated state of the impurity.

4. 1

 1 Exponents of the gaussian bond function basis from Ref. [65]. . . . . . . . 4.2 Equilibrium distances and corresponding potentials for the CaHe ground state 1 Σ + . Distances are in a 0 ( Å ) whereas energies are expressed in µhartree (cm -1 ). The notation aCVQZ and aV5Z are respectively for augcc-pCVQZ and aug-cc-pV5Z. The Ca core is 1s 2 . . . . . . . . . . . . . . . . 4.3 CaHe interaction energies calculated at different level of theories in combination of the aug-cc-pV5Z basis set for He and aug-cc-pCVQZ basis set for Ca. The complete set of bond functions (33211) has been used except in case a where no bond functions are used and case b where only the 332 bond functions have been used. All energies are in 10 -6 a.u. and are counter-poise corrected. All the Ca electrons are correlated, except the 1s electrons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Potential energy curve parameters, ie r e and ǫ . . . . . . . . . . . . . . . . 4.5 Dispersion coefficients in atomic units. . . . . . . . . . . . . . . . . . . . . 4.6 First vibrational energy levels for various ab initio methods. . . . . . . . . 5.1 Equilibrium distances, r e , and the corresponding potential energy values, ǫ, using different basis sets and bond functions. E 0 is the ground state energy and r 0 the ground state distance expectation value. * Only the 2 valence electrons of Mg have been correlated in the calculations. . . . . . . . . . .

(

  r) and V CCSDT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 HFD-B fit parameters and dispersion coefficients from the highest levels of calculation. Use of ab initio methods in combination of the largest basis sets for both helium ( aug-cc-pV5Z) and magnesium (aug-cc-pVQZ) with addition of the complete set of bond functions (33211). All data are in atomic units. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Dispersion parameters in atomic units. . . . . . . . . . . . . . . . . . . . . 5.7 Conventional CBS extrapolated binding energies of MgHe at an internuclear distance of 9.5 a 0 . All energies are in Hartree. . . . . . . . . . . . . . 5.8 Energies (in a.u.) of the MgHe ground state from Method I , Method II and Method III evaluated at an interatomic distance of 9.5 a 0 in combination with various basis sets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 CBS extrapolated binding energies. All energies are in atomic units. . . . . 5.10 Rotational levels for main magnesium isotopes, for v = 0 and possible rotational constants for the three main magnesium isotopes. Energy units are cm -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.1 λ values corresponding to different ab initio methods used for MgHe and CaHe. Large basis sets have been used in combination with bond functions for both systems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2 DMC mean energies for pure and corresponding doped helium clusters.

5. 3 47 5. 4 48 5. 5

 3474485 Graphic representation of the ground state potential energy surfaces of the MgHe molecule at the MP2-4, CCSD(T) and CCSDT level of calculations associated with the aug-cc-pV5Z basis set for He and aug-cc-pCVQZ for Mg and completed with the complete set of bond functions. . . . . . . . . Comparison between the aug-cc-VQZ and aug-cc-pCVQZ for the magnesium atom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Basis set superposition error as a function of interatomic distance at the CCSD(T) level of theory. The atomic basis sets used are aug-cc-pV5Z for He and aug-cc-pCVQZ for Mg. CP means that energies are counter-poised corrected. Numerical data are given in table A.9. . . . . . . . . . . . . . . 49 5.6 Complete basis set extrapolation for the correlation energy of the He (a), Mg (b), MgHe (c) using for each system the three fitting methods (cor1, cor2 and cor3). Ab initio points are performed at the method I level of theory. Fig.(d) shows the fitting of HF energies of He. In Fig. (a,b and

7. 3

 3 Schematic representation of the laser-depletion apparatus used for the pickup and depletion spectroscopy of molecules inside helium droplets (from Ref.[117]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .8.1 Evolution of the MgHe12 energy versus the number of walkers N w . Circles with error bars are the DMC points. The solid curve represents the fitting function of Eq. 8.2. The extrapolated energy is found to be -18.76(1) cm -1 8.2 3d and 2d helium density in the MgHe 12 droplet for N w = 100 . . . . . . . 8.3 3d and 2d helium density in the MgHe 12 cluster for N w = 12800 . . . . . . 8.4 Evolution of the MgHe 25 mean energy versus the time step length. . . . . . 8.5 Density (in Å -3 ) of He atoms in the MgHe 25 cluster for two time steps (in atomic units): 25 (circles) and 50. . . . . . . . . . . . . . . . . . . . . . . . 8.6 Evolution of the MgHe 12 mean energy versus the simulation length, ie, the number of steps, N steps . For all simulations, the time step, δτ was 25 a.u. . 8.7 Evolution of the helium probability density in the MgHe 12 droplet for N steps = 10 (squares) and 500 (triangles) . . . . . . . . . . . . . . . . . . . 8.8 Evolution of the MgHe 12 absolute energy versus number of blocks N b . . . 8.9 Evolution of the helium pair particle distribution function, π, in case of 50 blocks (triangles) and 1700 blocks (squares) in the MgHe 12 cluster. . . . . 8.10 Evolution of the MgHe 12 energy against seed numbers. . . . . . . . . . . . 8.11 Radial probability density of particles in the MgHe 12 cluster using different random number generator seeds. . . . . . . . . . . . . . . . . . . . . . . . 8.12 Typical pair particle trial wave function components : function drawn with circles represents φ jas/f er i .The curve with crosses shows a representation of φ jas i . Dashed line is the product φ jas/f er i × φ jas i which is the trial wave function for a couple of binary particles, and solid line is the square of the trial wave function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Fig. (b) shows the DMC binding energies corresponding to those potentials. . . . . . . . . 123 9.7 Fig.(a) : Mg probability density in the Mg@He n clusters. Fig.(b) : Density

131 9 .

 9 11 Radial probability densities of Mg in helium clusters whose size are respectively, from the left to the right of the figure, 18, 25, 32, 53, 83 and 170. These distributions are obtained with the helium pair potential from Ref.[35]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .132 9.12 Helium density in the Mg@He n clusters. In figures a, b, c and d n takes respectively 4, 6, 15 and 18. The helium pair potential of Jeziorska et al.[35] have been used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 9.13 Helium contour density in the Mg@He n clusters. In figures a, b, c and d n takes respectively 4, 6, 15 and 18. . . . . . . . . . . . . . . . . . . . . . . . 9.14 Helium density in cylinder coordinates for MgHe 12 and MgHe 100 clusters. The contours of the helium density are also plotted on the right of the figure. r and z are in Å and the helium density in Å -3 . . . . . . . . . . . . 9.15 Projection of the helium densities on the (rz) plane for different CaHe n clusters. In all figures, r and z are in Å . . . . . . . . . . . . . . . . . . . . 9.16 Schematic evolution of the position of the Mg atom in the He 32 cluster. Horizontal axis represents the evolution in imaginary time τ . Figure below this axis picture the projection of the helium density in the (rz) plane. Representative corresponding particle arrangements are also shown. . . . . 9.17 Structural difference between classical dynamic optimization of the MgHe 170 and a typical DMC configuration (picture on the right). . . . . . . . . . . . 9.18 Pair particle distance distributions for helium atoms in case of Mg doped with 18, 25, 32, 53, and 83 helium atoms. All curves are normalized to one. 9.19 Helium-Mg pair particle distance distributions for clusters with 12, 25, 32, 53, and 83 helium atoms. All curves are normalized to unity. . . . . . . . . 9.20 Evolution of the MgHe 50 and AgHe 50 system energies as a function of the distance between Mg (or Ag) and the helium atom group center. Bold circles and circles respectively depict the DMC absolute mean energies for MgHe 50 and AgHe 50 for a given constrained distance. Solid line represents the helium density for the He 50 cluster in particles per Å 3 . . . . . . . . . . 9.21 DMC absolute mean energies of the CaHe 50 cluster (bold circles) and the NaHe 50 system with respect to the distance between Ca (or Na) and the helium atom group center of mass. Solid line represents the helium density for the He 50 cluster in particles per Å 3 . . . . . . . . . . . . . . . . . . . . . 9.22 3d density helium density as a function of the fixed distance, d cst , for the MgHe 20 droplet. The density is expressed in particles per Å 3 . . . . . . . . 9.23 Projection of the helium density of MgHe 20 in the (rz) plane as a function of the fixed distance, d cst . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.24 Projection of the helium density in the (rz) plane for the MgHe 50 cluster as a function of the fixed distance, d cst . . . . . . . . . . . . . . . . . . . . . 9.25 Three-dimensional helium density (in Å -3 in the MgHe 50 cluster (figure a) and the rz projection of the density (figure b). . . . . . . . . . . . . . . . . 9.26 Fundamental rovibrational state (horizontal line which is shifted by 115 cm -1 ). The thick curve shows the fundamental rovibrational wave function (×15) and the thin curve represents the "constrained" potential shifted by 115 cm -1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 10.1 Mg 2 X 1 Σ + g potential energy surface. Solid curve represents the fitting function 10.1. Squares are the experimental RKR data [130]. Data are obtained from NIST [131]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 10.2 MgHe X 1 Σ + potential energy surface. Solid curve represents the fitting function (Eq. 5.3). Triangles our CCSDT potential. . . . . . . . . . . . . . 155 10.3 Convergence of the effective He 2 potential at T = 0.38 K. The deepest well curve represents the classical potential of Ref. [35], the successive lines the first iteration steps. Convergence is reached after 6 iterations. . . . . . . . 157 10.4 (a) : MSD (in Å 2) of helium depicted by the curve with points and the second curve is that of Mg. (b) : Brownian move of Mg in He 1998 projected in the (XY) plane. The coordinates have been shifted so that the initial position of Mg is the origin of the sytem axies. . . . . . . . . . . . . . . . . 159 10.5 Radial pair particle distribution of helium at T=0.38 K, convoluted with Φ 2 5 (r). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 10.6 (a) : Distance (D) of Mg from its initial position with respect to the time t. D is in Å . (b) : Density (ρ in Å -3 ) of the Mg atom in the Mg@He 1998 complex : curve with triangles is for a simulation length of 0.2 ns and the second one is from a simulation of 80 ns. r represents the distance from the center of mass (in Å ). . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 10.7 Helium density (ρ) in the MgHe 1998 cluster, convoluted with Φ 2 5 (r). ρ is in Å -3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 10.8 (a) : Distance (R in Å ) between the two Mg in the Mg 2 He 1997 cluster with respect to the time t. (b) : Zoom of the beginning of figure (a). . . . . . . 164 10.9 (a) : Motion of the two Mg atoms in Mg 2 He 1997 projected onto the (xy) plane. (b) : Projected motion of the Mg atom located near the center of the doplet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 A.1 Potential energy curves of the MgHe

Figure 4 . 1 :

 41 Figure 4.1: Potential energy curves of the CaHe1 Σ + state obtained at the CCSD(T) level of theory : circles picture the curves obtained by the use of cc-pVDZ and augcc-pVDZ respectively for calcium and helium. The second curves represents the use of cc-pVTZ for calcium and aug-cc-pVTZ for helium at the same level of theory. For both curves, the complete set of mid-bond functions (33211) has been used.

Figure 4 . 3 :

 43 Figure 4.3: X 1 Σ + CaHe potential energy curves derived from different ab initio approaches in combination with aug-cc-pV5Z for He and aug-cc-pCVQZ for Ca plus bond functions (33211): curve with triangles are the CCSD(T) results, squares are for MP2, stars for MP4. Dots are results from Ref. [68] (calculated at the CCSD(T) level with the wtMCP pseudo-potential for Ca) and smooth curve represents the unpublished data from W. Meyer.

Figure 5 . 1 :

 51 Figure 5.1:The core-valence correlation influence studied at the CCSD(T) level of theory in combination with the following basis sets : VQZ for both He and Mg, and aug-cc-pCVQZ ( aug-cc-pV5Z) for Mg (He). All curves are calculated taking into account the complete set of bond functions (33211). Notice that the CCSD(T) level with aug-cc-pCVQZ ( aug-cc-pCV5Z) for Mg (He) + BF is abbreviated as b1.

Figure 5 . 2 :

 52 Figure 5.2: Influence of bond functions considering the CCSD(T),b1 method.

Figure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.3: Graphic representation of the ground state potential energy surfaces of the MgHe molecule at the MP2-4, CCSD(T) and CCSDT level of calculations associated with the aug-cc-pV5Z basis set for He and aug-cc-pCVQZ for Mg and completed with the complete set of bond functions.

Figure 5 . 5 :

 55 Figure 5.5: Basis set superposition error as a function of interatomic distance at the CCSD(T) level of theory. The atomic basis sets used are aug-cc-pV5Z for He and aug-cc-pCVQZ for Mg. CP means that energies are counter-poised corrected. Numerical data are given in table A.9.

Fig. 5 .

 5 Fig. 5.7(d) shows the fitting of the HF energies obtained by method II , curve depicted by the solid line takes into account the double ζ energies whereas that represented by dots does not. The insert in figure 5.7(d) better shows the light difference between these two curves. This lightness is due to the fact that the fit is a weighted fit and the weight is very low for the double ζ energies (see section 5.5.2). Fig. 5.7(a,b and c) represent the

Figure 5 . 6 :

 56 Figure 5.6: Complete basis set extrapolation for the correlation energy of the He (a), Mg (b), MgHe (c) using for each system the three fitting methods (cor1, cor2 and cor3). Ab initio points are performed at the method I level of theory. Fig.(d) shows the fitting of HF energies of He. In Fig.(a,b and c) solid lines draw the results of the cor1 method (Eq. 5.6), dots those from the cor2 method (Eq. 5.7) and dashed lines those from the cor3 method (Eq. 5.8).

Figure 5 . 7 :

 57 Figure 5.7: CBS extrapolation for the correlation energy of the He (a), Mg (b), MgHe (c) using for each system the three fitting methods (cor1, cor2 and Eq. 5.5). Ab initio points are performed at the method II level of theory. Fig.(d)shows the fitting HF energies : Curve with dots does not take into account the double ζ energies. In Fig.(a,b and c) solid lines draw the results of the cor1 method (Eq. 5.6), dots those from the cor2 method (Eq. 5.7) and dashed lines those from (Eq. 5.8).

Figure 5 . 8 :

 58 Figure 5.8: CBS extrapolation for the correlation energy of the He (a), Mg (b), MgHe (c) using for each system the two fitting methods (cor1, cor2). Ab initio points are performed at the method III level of theory. In all figures, curves with dots and solid lines do not take into account the double ζ energies. In Fig.(a,b and c) solid and dashed lines draw the results of the cor1 method (Eq. 5.6). The other curves are from the cor2 method (Eq. 5.7).

walk ( 6 0 P rad r 2 dr = ∞ 0 P

 600 .60) and the He-Mg or He-Ca pair correlation function. Both expressions are normalized such that ∞ pair r 2 dr = 1(6.61)

Figure 7 . 1 :

 71 Figure 7.1: An expanded view of the OCS IR spectrum in pure 4 He droplets with N4 He = 6000 atoms (A) and in pure 3 He droplets with N3 He = 12000 atoms (B). The depletion is plotted versus the change in wave number with respect to the origin of the spectrum in (A) at ν 0 = 61.64 cm -1 . From Ref. [16].

Figure 7 . 2 :

 72 Figure 7.2: A series of OCS IR spectra in 3 He similar to Fig. 7.1 (B) but with increasing average numbers N4 He of added 4 He atoms [ N4 He =0 (A), 7 (B), 25 (C), 35 (D), 60 (E), and 100 (F)]. The change in wave number in the abscissa is with respect to ν 0 = 2061.71 cm -1 . From Ref. [16].

Figure 7 . 3 :

 73 Figure 7.3: Schematic representation of the laser-depletion apparatus used for the pickup and depletion spectroscopy of molecules inside helium droplets (from Ref. [117]).

Figure 8 . 1 :

 81 Figure 8.1: Evolution of the MgHe 12 energy versus the number of walkers N w . Circles with error bars are the DMC points. The solid curve represents the fitting function of Eq. 8.2. The extrapolated energy is found to be -18.76(1) cm -1.

Figure 8 . 2 :Figure 8 . 3 :

 8283 Figure 8.2: 3d and 2d helium density in the MgHe 12 droplet for N w = 100 .

Figure 8 . 4 :Figure 8 . 5 :Figure 8 . 6 :

 848586 Figure 8.4: Evolution of the MgHe 25 mean energy versus the time step length.

Figure 8 . 7 :

 87 Figure 8.7: Evolution of the helium probability density in the MgHe 12 droplet for N steps = 10 (squares) and 500 (triangles)

Figure 8 . 8 :

 88 Figure 8.8: Evolution of the MgHe 12 absolute energy versus number of blocks N b

Figure 8 . 9 : 1 Figure 8 . 10 :

 891810 Figure 8.9: Evolution of the helium pair particle distribution function, π, in case of 50 blocks (triangles) and 1700 blocks (squares) in the MgHe 12 cluster.

Figure 8 . 11 :

 811 Figure 8.11: Radial probability density of particles in the MgHe 12 cluster using different random number generator seeds.

Figure 8 . 12 :

 812 Figure 8.12: Typical pair particle trial wave function components : function drawn with circles represents φ jas/f er i .The curve with crosses shows a representation of φ jas i . Dashed line is the product φ jas/f er i × φ jas i which is the trial wave function for a couple of binary particles, and solid line is the square of the trial wave function.

Figure 8 . 13 :

 813 Figure 8.13: Radial helium densities (ρ) of the He 12 cluster : curves with triangles and diamonds show the densities respectively obtained with a non optimized wave function and an optimized one. Horizontal solid line is for the bulk helium density. ρ is in Å -3 and R com in Å .

Figure 9 . 1 :

 91 Figure 9.1: Comparison of the 3 1 P 0 1 ← 3 1 S 0 transition of Mg atoms picked up by He nanodroplets (thin line) and solvated in bulk liquid helium (thick line). This figure is reported in Ref. [28].

Figure 9 . 2 :

 92 Figure 9.2: Binding energies of the Mg@He N complexes when using helium pair potential of[35]. Solid line is the fitting function from Eq. 9.17.

Figure 9 . 3 :Figure 9 . 4 :

 9394 Figure 9.3: Absolute DMC energies obtained by using the Jeziorska pair potential for the helium atoms[35]. Circles are the energies of pure helium clusters (He n ) and triangles are that of the MgHe n clusters.

Figure 9

 9 Figure 9.5: Binding energies of the Mg@He N complexes when using helium pair potential of[127].

Figure 9 . 6 :

 96 Figure 9.6: Absolute DMC energies (a) of the CaHe N clusters (crosses) obtained by using the Aziz pair potential for the helium atoms [34] and our CaHe CCSD(T) potential. Circles depict the helium energies. Fig. (b) shows the DMC binding energies corresponding to those potentials.

Figure 9 . 7 :

 97 Figure 9.7: Fig.(a) : Mg probability density in the Mg@He n clusters. Fig.(b) : Density of the Mg atom in the Mg@He n complexes according to Ref. [31].

Figure 9 . 8 :

 98 Figure 9.8: Radial helium density (ρ in Å -3 ) in the MgHe 220 droplet (a) and density profiles for Mg@ 4 He N drops for N = 300, 500, 1000, 2000, 3000, and 5000 (from [32])(b).

Figure 9 . 9 :Figure 9 . 10 :

 99910 Figure 9.9: From [30]: ion yield curves for Ca, Sr, Ba, and Xe atoms picked up by a helium nanodroplet beam (on the left). The shape difference reflects the fact that the metal atoms are located at the droplet surface and the xenon atom is solvated inside the droplet. The same conclusion is true in the case of the Mg atom (the right of the figure).

Figure 9

 9 Figure 9.11:Radial probability densities of Mg in helium clusters whose size are respectively, from the left to the right of the figure, 18, 25, 32, 53, 83 and 170. These distributions are obtained with the helium pair potential from Ref.[35].

Figure 9 . 12 :Figure 9 . 13 :Figure 9 . 14 :Figure 9 . 15 :

 912913914915 Figure 9.12: Helium density in the Mg@He n clusters. In figures a, b, c and d n takes respectively 4, 6, 15 and 18. The helium pair potential of Jeziorska et al. [35] have been used.

Figure 9 . 16 :

 916 Figure 9.16: Schematic evolution of the position of the Mg atom in the He 32 cluster. Horizontal axis represents the evolution in imaginary time τ . Figure below this axis picture the projection of the helium density in the (rz) plane. Representative corresponding particle arrangements are also shown.

Figure 9 . 17 :Figure 9 . 18 :Figure 9 . 19 :

 917918919 Figure 9.17: Structural difference between classical dynamic optimization of the MgHe 170 and a typical DMC configuration (picture on the right).

Figure 9 . 20 :Figure 9 . 21 :

 920921 Figure 9.20: Evolution of the MgHe 50 and AgHe 50 system energies as a function of the distance between Mg (or Ag) and the helium atom group center. Bold circles and circles respectively depict the DMC absolute mean energies for MgHe 50 and AgHe 50 for a given constrained distance. Solid line represents the helium density for the He 50 cluster in particles per Å 3 .

Figure 9 . 22 :- 15 Figure 9 . 23 :- 15 Figure 9 . 24 :Figure 9 . 25 :

 9221592315924925 Figure 9.22: 3d density helium density as a function of the fixed distance, d cst , for the MgHe 20 droplet. The density is expressed in particles per Å 3 .

Figure 9 . 26 :

 926 Figure 9.26: Fundamental rovibrational state (horizontal line which is shifted by 115 cm -1 ). The thick curve shows the fundamental rovibrational wave function (×15) and the thin curve represents the "constrained" potential shifted by 115 cm -1 .

Figure 10 . 1 :

 101 Figure 10.1: Mg 2 X 1 Σ + g potential energy surface. Solid curve represents the fitting function 10.1. Squares are the experimental RKR data [130]. Data are obtained from NIST [131].

Figure 10 . 2 :

 102 Figure 10.2: MgHe X 1 Σ + potential energy surface. Solid curve represents the fitting function (Eq. 5.3). Triangles our CCSDT potential.

Figure 10 . 3 :

 103 Figure 10.3: Convergence of the effective He 2 potential at T = 0.38 K. The deepest well curve represents the classical potential of Ref. [35], the successive lines the first iteration steps. Convergence is reached after 6 iterations.

Figure 10 . 4 :

 104 Figure 10.4: (a) : MSD (in Å 2 ) of helium depicted by the curve with points and the second curve is that of Mg. (b) : Brownian move of Mg in He 1998 projected in the (XY) plane. The coordinates have been shifted so that the initial position of Mg is the origin of the sytem axies.

Figure 10 . 5 :

 105 Figure 10.5: Radial pair particle distribution of helium at T=0.38 K, convoluted with Φ 2 5 (r).

Figure 10 . 6 :

 106 Figure 10.6: (a) : Distance (D) of Mg from its initial position with respect to the time t. D is in Å . (b) : Density (ρ in Å -3 ) of the Mg atom in the Mg@He 1998 complex : curve with triangles is for a simulation length of 0.2 ns and the second one is from a simulation of 80 ns. r represents the distance from the center of mass (in Å ).

Figure 10 . 7 :

 107 Figure 10.7: Helium density (ρ) in the MgHe 1998 cluster, convoluted with Φ 2 5 (r). ρ is in Å -3 .

Figure 10 . 8 :

 108 Figure 10.8: (a) : Distance (R in Å ) between the two Mg in the Mg 2 He 1997 cluster with respect to the time t. (b) : Zoom of the beginning of figure (a).

Figure 10 . 9 :

 109 Figure 10.9: (a) : Motion of the two Mg atoms in Mg 2 He 1997 projected onto the (xy) plane. (b) : Projected motion of the Mg atom located near the center of the doplet.

Figure A. 1 :

 1 Figure A.1: Potential energy curves of the MgHe1 Σ + state obtained at different levels of theories. For all methods, the basis sets used are aug-cc-pV5Z for He and aug-cc-pCVQZ for Mg in combination with the complete set of bond functions (33211). The BSSE has only been taken into account for the bond functions.

  5.5).

  Figure B.1:Radial probability density of Mg in the MgHe 53 cluster for different interaction potentials. The reference potential is the CCSDT one. Increasing the CCSDT potential by only 10% is sufficient to give a fully solvated state of the impurity.

  Radial probability density of Mg in the MgHe 53 cluster for different interaction potentials. The reference potential is the CCSDT one. Increasing the CCSDT potential by only 10% is sufficient to give a fully solvated state of the impurity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184

1

Σ + state obtained at different levels of theories. For all methods, the basis sets used are aug-cc-pV5Z for He and aug-cc-pCVQZ for Mg in combination with the complete set of bond functions (33211). The BSSE has only been taken into account for the bond functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 B.1

Table 4 . 1 :

 41 Exponents of the gaussian bond function basis from Ref.[65].

Table 4 . 2 :

 42 Equilibrium distances and corresponding potentials for the CaHe ground state 1 Σ + . Distances are in a 0 ( Å ) whereas energies are expressed in µhartree (cm-1

Table 4 . 3 :

 43 CaHe interaction energies calculated at different level of theories in combination of the aug-cc-pV5Z basis set for He and aug-cc-pCVQZ basis set for Ca. The complete set of bond functions (33211) has been used except in case a where no bond functions are used and case b where only the 332 bond functions have been used. All energies are in 10

-6

a.u. and are counter-poise corrected. All the Ca electrons are correlated, except the 1s electrons.

Table 4 . 4 :

 44 Potential energy curve parameters, ie r e and ǫ

	METHODS	C 6	C 8	C 10
	MP2 MP4 CCSD(T) calc. [68] exp. [81]	37.9 37.3 29.4 36.9	2.6 × 10 3 2.53 × 10 3 4.3 × 10 3 2.51 × 10 3	1.3 × 10 4 8.6 × 10 4 6. × 10 5 -

45.1-48.5 (1.48 -2.19) × 10 3 (1.08 -1.29) × 10 5

Table 4 .

 4

5: Dispersion coefficients in atomic units.

Table 4 . 6 :

 46 First vibrational energy levels for various ab initio methods.4.6 Vibrational levels of the CaHe1 Σ + state Variational vibrational calculations have been carried out considering the CCSD(T), MP2 and MP4 ground state potential energy curve of the CaHe molecule and using a basis set of 200 optimized Laguerre functions. There is evidence of, at least, one bound state whose energy is reported in Table

Table 5 .

 5

1: Equilibrium distances, r e , and the corresponding potential energy values, ǫ, using different basis sets and bond functions. E 0 is the ground state energy and r 0 the ground state distance expectation value. * Only the 2 valence electrons of Mg have been correlated in the calculations.

Table 5 . 2 :

 52 Ground state MgHe interaction energies calculated at the CCSD(T) level of theory in combination of the aug-cc-pV5Z basis set for He and the indicated basis set for Mg, except in cases1,2 where the basis indicated is also used for He. In all cases, only the Mg 1s electrons have been kept frozen in calculations, exception made for case1 where only the Mg 3s electrons have been correlated. The complete set of bond functions (33211) has been used except in case a where no bond functions are used and case

		0
		-5
	V/µhartree	-10 -15				Mg core = 1s 2 , VQZ Mg core = 1s 2 2s 2 2p 6 , VQZ Mg core = 1s 2 , b1
							Mg core = 1s	2 2s	2 2p 6 , b1
		-20
		-25
		8	9	10	11	12	13	14	15	16
						r/a 0

b

where only the 332 bond functions have been used. All energies are in 10 -6 a.u. and are counter-poise corrected.

Table 5 . 3 :

 53 Equilibrium distances, r e , and the corresponding potential energy values, ǫ, using different sets of bond functions. All the the electrons of Mg have been correlated in the calculations, except the 1s electrons. E 0 is the ground state energy and r 0 the ground state distance expectation value.

and -23.04 (×10

-6

Table 5 . 4 :

 54 Ab initio counter-poise corrected interaction energies : calculations have been carried out with an aug-cc-pV5Z (aug-cc-pVQZ) for He (Mg) basis set in combination of the complete set of bond function (33211). All energies are in ×10 -6 a.u. ∆V (r) is the relative error between V CCSDT

	Methods	A	β	C 6	C 8	C 10	RM S
	MP2	1.58×10 1 1.37 1.86×10 1 1.15×10 3 7.20×10 4	0.20
	MP4	1.90×10 1 1.39 2.09×10 1 1.10×10 3 9.14×10 4	0.30
	CCSD(T) 1.66×10 1 1.38 1.87×10 1 1.15×10 3 6.95×10 4	0.15
	CCSDT	1.60×10 1 1.38 1.89×10 1 1.28×10 3 7.03×10 4	0.03

f it (r) and V CCSDT .

Table 5 . 5 :

 55 HFD-B fit parameters and dispersion coefficients from the highest levels of calculation. Use of ab initio methods in combination of the largest basis sets for both helium ( aug-cc-pV5Z) and magnesium (aug-cc-pVQZ) with addition of the complete set of bond functions (33211). All data are in atomic units.

		C 6	C 8	C 10
	Standard et al. [81] (2.11 -2.21)×10 1 (0.843 -0.924)×10 3 (3.55 -4.19)×10 4
	R. Hinde [77]	2.0×10 1	1.1×10 3	-
	Partridge [66]	1.978×10 1	-	-
	This work :
	MP2	1.86×10 1	1.15×10 3	7.20×10 4
	MP4	2.09×10 1	1.10×10 3	9.14×10 4
	CCSD(T)	1.87×10 1	1.15×10 3	6.95×10 4
	CCSDT	1.89×10 1	1.28×10 3	7.03×10 4

Table 5 . 6 :

 56 Dispersion parameters in atomic units.

Table 5 . 7 :

 57 Conventional CBS extrapolated binding energies of MgHe at an internuclear distance of 9.5 a 0 . All energies are in Hartree.

						-0.2
		-0.034				-0.22
		-0.036				-0.24
						-0.26
	E/a.u	-0.038			E/a.u	-0.28
		-0.04				-0.3
						-0.32
		-0.042				-0.34
		-0.044				-0.36
		apVDZ	apVTZ	apVQZ	apV5Z
			n			apCVDZ	apCVTZ	apCVQZ	apCV5Z
			(a)			(b)
				-0.22
				-0.24
				-0.26
				-0.28
				-0.3
			E/a.u	-0.32
				-0.34
				-0.36
				-0.38
				-0.4
				-0.42
				apCVDZ	apCVTZ	apCVQZ	apCV5Z
					(c)

Table 5 . 8 :

 58 Energies (in a.u.) of the MgHe ground state from Method I , Method II and

  Mg 4 He 25 Mg 4 He 26 Mg 4

				He
	0	-0.902	-0.909	-0.915
	1	-0.651	-0.659	-0.666
	2	-0.176	-0.185	-0.193
	B 0	0.1269	0.1264	0.1260

  6.6.1 VMC call initstate(r old ) ! Initialize the ensemble of states Psi old =Ψ trial (r old ) * Ψ trial (r old ) ! Evaluate Ψ trial for i=1,nblocks ! Loop over blocks for j=1,nsteps ! Loop over time steps for k=1,nwalkers ! Loop over walkers call sample (r old , r new , T new ) ! Sample new state from drifted Gaussian Psi new=Ψ trial (r new ) * Ψ trial (r new ) ! Evaluate Ψ trial call sample (r new , r old , T old ) ! Find transition probability for going backward A = P si 2 initstate(r old ) ! Initialize the ensemble of states P si old = Ψ trial (r old ) ! Evaluate Ψ trial D old = Drif t(r old ) ! Evaluate ∇Ψ trial for i=1,nblocks ! Loop over blocks for j=1,nsteps ! Loop over time steps for k=1,nwalkers ! Loop over walkers call sample (r old , r new , T new , D old ) ! Sample new state from drifted Gaussian Psi new=Ψ trial (r new ) ! Evaluate Ψ trial D new = Drif t(r new ) ! Evaluate ∇Ψ trial call sample (r new , r old , T old , D new ) ! Find transition probability for going backward A = P si 2 old /T old ! Determine the probability of acceptance. if(A > rand () ) then r old = r new P si old = P si new naccept = naccept + 1 endif weight(r old ) = weight(r old ) * exp(-dτ * (E l (r old -E ref ))

	enddo
	call reweight(r old )
	call averages(r old )
	enddo
	call output()
	enddo	new /Tnew P si 2 old /T old ! Determine the probability of acceptance.
	if(A > rand () ) then
	r old = r new
	P si old = P si new
	endif
	enddo
	call averages(r old )
	enddo
	call output()
	enddo
	6.6.2 DMC
		new /Tnew P si 2

call

Table 9 .

 9 1: λ values corresponding to different ab initio methods used for MgHe and CaHe. Large basis sets have been used in combination with bond functions for both systems. the threshold value of 1.90. For Mg and Ca, the λ parameter values are given in table 9.1.

Table

Table 9 . 3 :

 93 DMC mean energies for pure and corresponding Mg doped helium clusters.

  CaHe n ) δE(CaHe n ) E(He n ) δE(He n )

	n	δt	E(E n	δE n
	4	25.00	-3.08	0.01	-0.40	0.00	-2.68 0.01
	6	25.00	-5.52	0.01	-1.65	0.00	-3.88 0.01
	8	25.00	-8.41	0.01	-3.58	0.00	-4.82 0.01
	12	25.00	-15.29	0.01	-8.85	0.03	-6.44 0.03
	15	50.00	-21.19	0.01	-13.89	0.03	-7.30 0.04
	18	25.00	-27.56	0.02	-19.48	0.02	-8.08 0.02
	25	25.00	-43.84	0.02	-34.20	0.05	-9.64 0.06
	32	25.00	-61.57	0.03	-50.93	0.04	-10.65 0.04
	53	50.00	-120.54	0.06	-107.74	0.05	-12.80 0.08
	75	50.00	-187.62	0.21	-173.42	0.05	-14.20 0.22
	83	50.00	-213.44	0.12	-199.02	0.13	-14.42 0.18
	100 100.00	-269.66	0.11	-253.38	0.17	-16.29 0.20
	133 100.00	-379.59	0.22	-365.59	0.18	-14.00 0.28
	170 50.00	-511.18	0.23	-494.71	0.32	-16.47 0.40
	190 100.00	-583.10	0.18	-566.80	0.32	-16.30 0.37

Table 9 . 4 :

 94 DMC mean energies for pure and corresponding Ca doped helium clusters.

Table A . 1 :

 A1 Hartree-Fock energies for the MgHe complex and its components. The binding energies -E Bind.

  Bind. /10 -6 a.u. 7.00 -199.955534662055 -2.903205125395 -202.858415185247 324.60220298968 7.50 -199.955477887023 -2.903204352680 -202.858594280167 87.959535982395 8.00 -199.955424584146 -2.903203752537 -202.858642421821 -14.085138015484 8.50 -199.955379614875 -2.903203276150 -202.858638679048 -55.788022999436 8.75 -199.955358587930 -2.903203070835 -202.858625814753 -64.155988004355 9.00 -199.955337613764 -2.903202880919 -202.858607147478 -66.652795024158 9.25 -199.955316588953 -2.903202701962 -202.858584291235 -65.000319994013 9.50 -199.955295788876 -2.903202530272 -202.858559201765 -60.882617003876 9.75 -199.955275669493 -2.903202363432 -202.858533709197 -55.676271985216 10.00 -199.955256718110 -2.903202200714 -202.858509232428 -50.313603983998 10.50 -199.955223936874 -2.903201893204 -202.858466620782 -40.790703990545 11.00 -199.955199386483 -2.903201628751 -202.858434042896 -33.027661991091 11.50 -199.955182477596 -2.903201426434 -202.858409968980 -26.064950001458 12.00 -199.955170627296 -2.903201283642 -202.858391495934 -19.584995987642 12.50 -199.955161069996 -2.903201176077 -202.858376278163 -14.032090012161 13.00 -199.955152260942 -2.903201075583 -202.858363149712 -9.8131870189810 14.00 -199.955136557298 -2.903200861506 -202.858342531413 -5.1126089966047 15.00 -199.955125928129 -2.903200706144 -202.858329937272 -3.3029989960908 16.00 -199.955120425703 -2.903200658254 -202.858323521320 -2.4373629874396 20.00 -199.955110521366 -2.903200562093 -202.858311763496 -0.6800369951598 26.00 -199.955107021423 -2.903200530034 -202.858307645879 -0.0944219911147

	R/a 0 E Table E CCSDT Mg /a.u. E CCSDT He /a.u. E CCSDT MgHe /a.u.

  a.u. E Bind. /10 -6 a.u.

	R/a 0 E	CCSD(T) Mg	/a.u. E CCSD(T) He	CCSD(T) /a.u. E MgHe
	6.00	-199.9557854	-2.90311127	-202.8569799	1916.69871316
	7.00	-199.9556320	-2.90310837	-202.8583107	429.69749106
	7.50	-199.9555626	-2.90310691	-202.8584933	176.20733587
	8.00	-199.9554927	-2.90310569	-202.8585421	56.32209377
	8.50	-199.9554397	-2.90310471	-202.8585406	3.78292432
	8.75	-199.9554156	-2.90310427	-202.8585288	-8.92175149
	9.00	-199.9553907	-2.90310388	-202.8585108	-16.28438787
	9.25	-199.9553649	-2.90310353	-202.8584885	-20.08607066
	9.50	-199.9553391	-2.90310322	-202.8584639	-21.56724813
	9.75	-199.9553144	-2.90310294	-202.8584389	-21.58345853
	10.00	-199.9552915	-2.90310268	-202.8584149	-20.71706228
	10.50	-199.9552530	-2.90310224	-202.8583730	-17.76313960
	11.00	-199.9552244	-2.90310189	-202.8583407	-14.45797110
	11.50	-199.9552037	-2.90310164	-202.8583168	-11.47596248
	12.00	-199.9551880	-2.90310147	-202.8582984	-9.00945938
	12.50	-199.9551749	-2.90310135	-202.8582833	-7.05210434
	13.00	-199.9551635	-2.90310125	-202.8582702	-5.53011091
	14.00	-199.9551452	-2.90310102	-202.8582497	-3.45354422
	15.00	-199.9551341	-2.90310085	-202.8582371	-2.21618834
	16.00	-199.9551285	-2.90310079	-202.8582308	-1.46303134
	17.00	-199.9551254	-2.90310077	-202.8582271	-0.99144577
	18.00	-199.9551227	-2.90310074	-202.8582241	-0.63654772
	19.00	-199.9551202	-2.90310070	-202.8582213	-0.43852737
	20.00	-199.9551180	-2.90310067	-202.8582190	-0.30475255
	22.00	-199.9551155	-2.90310064	-202.8582163	-0.14681694
	24.00	-199.9551147	-2.90310064	-202.8582154	-0.11288577
	26.00	-199.9551142	-2.90310063	-202.8582149	-0.06755062
	28.00	-199.9551136	-2.90310063	-202.8582143	-0.04159176
	30.00	-199.9551132	-2.90310063	-202.8582139	-0.02592751
	35.00	-199.9551129	-2.90310063	-202.8582135	-0.00782336
	40.00	-199.9551128	-2.90310063	-202.8582135	-0.00142214

/

Table A . 6 :

 A6 CCSD(T) energies for the MgHe complex and its components. The binding energies are counter-poise corrected. The Mg basis was aug-cc-pCVQZ and the He one was aug-cc-pV5Z. The complete set of bond functions (BFs) has also been used in combination of the atomic basis sets.

	R/a 0	E noCP	E CP	E BSSE
	6.00	1889.81862	1936.35382	46.5352
	7.00	423.252521	440.157937	16.905416
	7.50	171.553588	183.7718	12.218212
	8.00	51.3416022	61.8233424	10.4817402
	8.50	-1.05002917	7.76337396	8.81340313
	8.75	-13.3567173	-5.54162708	7.81509022
	9.00	-20.2008464	-13.4236438	6.7772026
	9.25	-23.4470477	-17.6720617	5.774986
	9.50	-24.406826	-19.5366587	4.8701673
	9.75	-23.9822784	-19.881007	4.1012714
	10.00	-22.777333	-19.2945593	3.4827737
	10.50	-19.4414472	-16.7801975	2.6612497
	11.00	-16.0209862	-13.7877397	2.2332465
	11.50	-13.0092176	-11.024215	1.9850026
	12.00	-10.475111	-8.70738566	1.76772534
	12.50	-8.37262274	-6.85118982	1.52143292
	13.00	-6.64590317	-5.39658208	1.24932109
	14.00	-4.13928427	-3.39377789	0.74550638
	15.00	-2.62046674	-2.18853907	0.43192767
	16.00	-1.76299457	-1.44978821	0.31320636
	17.00	-1.28494861	-0.985076269	0.299872341
	18.00	-0.991451572	-0.684701007	0.306750565
	19.00	-0.777463413	-0.485195926	0.292267487
	20.00	-0.602005429	-0.349752355	0.252253074
	22.00	-0.335068705	-0.189417051	0.145651654
	24.00	-0.172353708	-0.106810745	0.065542963
	26.00	-0.0862668799	-0.0615608258 0.0247060541
	28.00	-0.0437262315	-0.0356262149 0.0081000166
	30.00	-0.0223201084	-0.0199654551 0.0023546533
	35.00 -0.00193281835 -0.00186446414 0.0000683542
	Table

Equation 3.72 follows the convention which means that E bind is zero for a completely dissociated molecule.

Whereas, ab initio methods have a high-order polynomial scaling with system size, for example the coupled cluster methods, scale, at least as N 6 .

However more intelligent choices for T (r ′ → r) can be used to increase the probability of accepting an attempted move and, therefore, to improve the algorithm efficiency.

In this work, the acceptance rates were typically between 98% and more than 99%, depending on the droplet size and the time step used in our simulation.

Remerciements

Je veux remercier, tout d'abord, Mme Gilberte Chambaud de m'avoir proposé l'opportunité de faire une thèse au LCT.

Appendix B

Position of Mg for several MgHe potentials