J L Lemaire

G Vidali

S Baouche

M Chehrouri

H Chaabouni

H Mokrane

COMPETING MECHANISMS OF MOLECULAR HYDROGEN FORMATION IN CONDITIONS RELEVANT TO THE INTERSTELLAR MEDIUM

Keywords: astrochemistry -ISM, atoms -ISM, molecules -methods, laboratory -molecular processes

The most efficient mechanism of the formation of molecular hydrogen in the current universe is by association of hydrogen atoms on the surface of interstellar dust grains. The details of the processes of its formation and release from the grain are of great importance in the physical and chemical evolution of the space environments where it takes place. The main puzzle is still the fate of the 4.5 eV released in H 2 formation and whether it goes into internal energy (rovibrational excitation), translational kinetic energy, or heating of the grain. The modality of the release of this energy affects the dynamics of the ISM and its evolution toward star formation. We present results of the detection of the rovibrational states of the just-formed H 2 as it leaves the surface of a silicate. We find that rovibrationally excited molecules are ejected into the gas phase immediately after formation over a much wider range of grain temperatures than anticipated. Our results can be explained by the presence of two mechanisms of molecule formation that operate in partially overlapping ranges of grain temperature. A preliminary analysis of the relative importance of these two mechanisms is given. These unexpected findings, which will be complemented with experiments on the influence of factors such as silicate morphology, should be of great interest to the astrophysics and astrochemistry communities.

INTRODUCTION

As the most abundant molecule in the universe, hydrogen influences both the way stars are formed and the chemical evolution of interstellar clouds. Although it was recognized long ago that the formation of molecular hydrogen most likely occurs on dust grains (Hollenbach & Salpeter 1971), as the radiative recombination involves spin-forbidden transition from the protomolecule to the ground state, it is only in the last decade that experiments were done aimed at understanding the mechanism and measuring the rates of H 2 formation on dust grain analogs in simulated space environments (Pirronello et al. 1997).

The energy that the nascent molecule brings into the interstellar medium (ISM) is of interest in the overall evolution of space environments and in programs of detection of excited H 2 in space. Measurements of the translational energy of HD from amorphous water-ice showed that the just-formed HD can carry energy into the ISM and could be detected (Manicòetal.2001). The observed abundance of warm H 2 gas in outer layers of photodissociation regions (PDRs) is hard to explain with current models, and contributions of excited H 2 from grains have been invoked (Habart et al. 2004). Searches of rovibrational excitation of H 2 that could be ascribed to "UV formation pumping" in PDR and dense clouds yielded no positive detection (Tinéetal. 2003;Lemaire & Field 2001;Thi et al. 2009).

Laboratory studies can help observational searches by assessing how the energy is distributed and by characterizing rovibrational excitations upon H 2 formation so to distinguish it from other excitation mechanisms, such as UV excitation, collisions, and shocks.

Results of experiments of H 2 formation on surfaces of polycrystalline and amorphous silicates (Pirronello et al. 1997; 3 Current address: Physics Department, Syracuse University, Syracuse, NY 13244, USA V i d a l ie ta l . 2007,2009), amorphous water-ice (Roser et al. 2003;Hornekaer et al. 2003;Amiaud et al. 2007), highly oriented pyrolitic graphite (HOPG; Islam et al. 2007;Latimer et al. 2008), and amorphous carbon (Katz et al. 1999) are consistent with H 2 formation via the Langmuir-Hinshelwood (L-H) mechanism or possibly the Harris-Kasemo (H-K) "hot atom" mechanism (Harris & Kasemo 1981), at least when the number of atoms on the surface is small. While the L-H mechanism is known to occur in the vast majority of chemical reactions at surfaces, the other involves the superthermal motion of atoms as they land on the surface. In this case, atoms are not accommodated by the surface, and the new molecule can carry substantial energy in translational and rovibrational excitations. Another mechanism, the Eley-Rideal reaction, consists in the direct reaction, without prior energy accommodation, of an atom from the gas phase with one on the surface. This mechanism is operative in H 2 / HD formation in H-loaded amorphous carbon (Mennella 2008). Photodissociation of amorphous solid water also involves other mechanisms to yield both cold and hot H 2 products (Yabushita et al. 2008).

Most of the prior work on H 2 formation via the H-K mechanism was done on metals and semiconductors, where there is large condensation energy (Kammler et al. 2000). It has not been established whether or how the H-K mechanism operates in H 2 formation on dust grain analogs at low temperature (<20 K). Experiments showed that rovibrationally hot molecules formed on ice leave the surface de-excited (Congiu et al. 2009). On HOPG (Islam et al. 2007;Latimer et al. 2008;Zecho et al. 2002), a large fraction of molecules desorb in the ν = 4 state, but in this case the measurements were done at higher sample temperature (>15 K) than in the other studies. This leads to two pressing questions: (1) whether and how silicates, which are a major constituent of interstellar dust and shown to catalyze H 2 formation at low temperature via the L-H mechanism (Vidali et al. 2009), are also involved in the formation of highly excited H 2 ; and (2) over which grain temperature range the formation of excited H 2 occurs, since this is relevant to understanding H 2 formation in ISM regions where the L-H mechanism does not work because of the high temperature of grains (Habart et al. 2004;Cazaux & Tielens 2004).

We report here the detection of rovibrationally hot D 2 formed on the surface of a silicate (deuterium is used because hydrogen is the main background gas in a well-baked ultra-high vacuum apparatus); we show that it occurs in a surface temperature range that, while consistent with the weak interaction between D and the silicate, is much wider than thermal programmed desorption (TPD) experiments indicated (Perets et al. 2007). We propose that it must be ascribed to the H-K mechanism competing with the L-H mechanism in the formation of D 2 .

EXPERIMENTAL METHODS AND RESULTS

Experiments are performed using the setup of Figure 1.W e summarize here the major features; details are given in Baouche (2004), Amiaud (2006), andMatar (2009).

The apparatus consists of an ultra-high vacuum stainless steel chamber (<10 -10 mbar). At its center, a sample holder is thermally connected to the cold finger extension of a closedcycle He cryostat. The temperature (in the 5.5-350 K range) is measured with two calibrated Si diodes, one clamped on the sample and another on the tip of the second stage of the cryostat where a resistive heater is located. In a TPD, the combination of the length and cooling power of the cold finger yields a delayed desorption of gases adsorbed on the sample supports. This allows us to separate the contribution of the sample of the one from other parts of the cold finger. An amorphous olivinetype silicate film was kindly provided by Dr. D'Hendecourt (IAS Orsay) by thermal evaporation of San Carlos olivine (Djouadi et al. 2005) onto a gold-coated substrate. Its properties have been checked by infrared spectroscopy on a KBr substrate coated by a silicate film at the same time as the copper sample (Davoisne et al. 2008).

The atomic jet is prepared in a microwave cavity in a triply differentially pumped line. The beam is bent in a teflon tube connected to a thick aluminum cone-shaped nozzle attached to a closed-cycle He cryostat. In this way, we ensure that no light or excited molecules created in the plasma reach the surface. We checked that no excited molecules are detected in the beam, either at room temperature or at 50 K. The experimentally estimated D 2 beam flux is ∼9 × 10 12 molecules cm -2 s -1 .TheD 2 dissociation is checked regularly during the experiments (typically 80%), as measured by a quadrupole mass spectrometer (QMS). During the irradiation phase, the sample is exposed to ∼1.35 × 10 13 D-atoms cm -2 s -1 (Amiaud et al. 2007;Accolla et al. 2010). The QMS is moved above the sample to avoid interferences with the time-of-flight (TOF) spectrometer and to record the background pressure during both deposition and TPD experiments. TPDs are performed usinga10Kminute -1 linear heating ramp, and the D 2 , HD, and H 2 signals are recorded.

To measure the rovibrational excitation of the molecules leaving the surface, we used resonantly enhanced multi-photon ionization (REMPI) coupled with TOF mass spectroscopy. A UV laser beam is tuned to ionize and detect selectively a given rovibrational state via 2-photon absorption through a virtual state followed by ionization: (2+1) REMPI. In this way, we sample the population of D 2 formed in the v ′′ = 4, J ′′ = 2 rovibrationally excited state of the ground state. This transition has been chosen for its suitable Franck-Condon factor (Fantz & Wunderlich 2004). A tunable dye laser pumped with the second harmonic of an ND:YAG laser (2 J/20 Hz) followed by frequency mixing gives ∼222.6 nm photons with in red, shows the formation of D 2 in the v ′′ = 4, J ′′ = 2 excited state when irradiating the sample at 12 K with D atoms from time t = 0 to 600 s. Blue (scale offset): QMS D 2 signal; its rise at t = 0 is due to the increased gas load from the beam. Orange: sample temperature. During the temperature ramp of the TPD (beam valve closed), a sharp peak of D 2 molecules desorbing from the surface is seen whereas no signal appears on the REMPI (the second broad and smaller peak during the TPD is an artifact due the cryostat, see the text). When the temperature reaches 70 K, the irradiation resumes. The QMS signal goes to the same level as at t = 0 (which is due to the extra D 2 gas load from the beam) while the REMPI signal is almost negligible.

200-300 µJ pulse -1 . The laser beam is focused ∼4m ma w a y from the center of the surface. At the exit of the main chamber, the shot-by-shot laser pulse integrated power appears rather stable (∼7%) over a few hours. The REMPI-TOF signal is recorded and integrated. Correcting the REMPI-TOF data using a power law to account for its nonlinear dependence with laser power made little difference, as the main source of noise is the random desorption of excited molecules.

An ancillary setup is composed of a hot molecular hydrogen lamp and a TOF spectrometer with characteristics very similar to the one in the main apparatus (Malmasson 1994). The source provides rovibrationally excited states (Schermann et al. 1994) by dissociating hydrogen on a 2900 K tungsten filament followed by recombination and thermalization on the cool walls. The laser light is diverted to this setup for wavelength adjustment and REMPI signal intensity calibration.

The experiments are conducted as follows: once the silicate has reached a given temperature, it is irradiated with a beam of D and some undissociated D 2 .D 2 , HD, and H 2 signals are detected by the QMS while the REMPI-TOF simultaneously measures the population in a given D 2 state. After 10 minutes of irradiation (5 minutes for temperatures >16 K), the separation valve is closed and immediately after, a TPD is started using the same diagnostics. We investigated the formation and excitation of D 2 molecules at 13 sample temperatures between 5.5 and 70 K. The REMPI and QMS signals obtained when the sample is irradiated at 12 K and then brought to 70 K are shown in Figure 2. We observe the following main features: (1) rovibrationally excited molecules are formed as soon as the sample is irradiated and (2) desorption in the TPD experiment yields no REMPI signal in the probed state. This means that there are no highly excited molecules leaving the surface during the TPD. Figure 3 shows the REMPI and QMS signals during D irradiation with the sample at 36 K. Surprisingly, we still detect the formation of D 2 in an excited state at this temperature, while there is no QMS-TPD signal because the residence time of atoms is too short. This is in agreement with experiments of HD formation on amorphous silicates using H and D beams (Vidali et al. 2009), where there is little or no TPD signal of HD after irradiation at sample temperature >28 K. At 70 K, there 2 with the sample at 36 K during irradiation. There is no TPD-QMS signal of D 2 from the surface during the temperature ramp (t > 360 s) as no deuterium stick on the surface at such temperature.

is almost no detectable REMPI signal during irradiation (see Figure 2).

Figure 4 shows the integrated REMPI and TPD signals for D irradiation at different sample temperatures. The maximum formation efficiency in both cases is at T ∼ 10 K. From this figure we find that (1) the rovibrationally excited molecules are detected over a range of sample temperature that is considerably wider than the one for the TPD signal. The decay of the REMPI signal versus sample temperature at irradiation is markedly different from the TPD signal, pointing to different mechanisms of D 2 formation (see below). (2) For T < 10 K, there is a diminishing REMPI signal while the QMS signal is slightly lower only at 5.5 K. This corresponds to a more efficient accommodation of the incoming atoms and to a slower diffusion time on the surface. Both the QMS and REMPI signals increase over time since the beginning of the irradiation (see Figure 2). This is also observed at temperatures up to 18 K where the signal quickly reaches a constant value. This is due to the changing sticking coefficient of D atoms as the fraction of the surface covered by D and D 2 increases. At high temperature (see Figure 3), the REMPI signal reaches steady state immediately. Experiments of D 2 sticking on a single crystal forsterite show that the sticking probability has a sharp transition at the sample temperature of 15-16 K (P. Frank et al. 2010, private communication).

ANALYSIS OF THE EXPERIMENTAL RESULTS

At P = 4 × 10 -4 mbar in the test chamber, there are ∼1.1 × 10 13 molecules cm -3 in all vibrational and rotational states. Using 5.8 × 10 -3 as the ratio of the populations in v ′′ = 4, (all J ′′ ) and in all states (Schermann et al. 1994), we get 6.4 × 10 8 molecules cm -3 in v ′′ = 4. The ratio of TOF-REMPI v ′′ = 4, J ′′ = 2 signals in the ultra-high vacuum chamber (sample at 12 K) and in the test chamber is 6.3 × 10 -3 .

For the population of D 2 desorbing from the silicate at 12 K in v ′′ = 4, the J ′′ excited state has 6.4 × 10 8 × (6.3 × 10 -3 ) = 4 × 10 6 molecules cm -3 . Of these, 1.1 × 10 6 molecules cm -3 are in the v ′′ = 4, J ′′ = 2 state, where for lack of other data, we used the results for HD on graphite at 15 K (Islam et al. 2007): of the molecules in the (v ′′ = 4, all J) states 27.5% are in J ′′ = 2. The density in all (v, J) states in the incoming beam is given by flux/velocity (v = 5.6 × 10 4 cm s -1 at 50 K). The D 2 flux reported above gives, at 80% dissociation, 7.2 × 10 12 molecules cm -2 s -1 and a density of ∼1.3 × 10 8 cm -3 , where we assumed that all atoms stick and recombine on the sample. This is the highest number of molecules that can be formed by whatever mechanism. In the case where they all formed by the H-K mechanism and assuming that they leave the surface in any direction, the number density that could be detected depends on the size of the laser's active region for ionization. Assuming a Gaussian beam, the confocal parameter b = 7 cm and the beam waist 2w 0 = 100 µm. From the beam waist the surface is seen under an angle of ∼41 • , allowing us to detect 41/180 of the molecules and a density of ∼2.9 × 10 7 molecules cm -3 in all states (including v ′′ = 0). Referring again to the results of Latimer et al. (2008) at 15 K, assuming the same ratio of 46% of v ′′ = 4 versus all v ′′ = 1-7 for silicates, we get ∼8.7 × 10 6 molecules cm -3 in excited states (v ′′ = 1-7) and ∼2.1 × 10 7 molecules cm -3 formed by whatever mechanism in v ′′ = 0 of the ground state. Some of them might be released immediately in the gas phase and some stick to the surface. The net result is that ∼30% of the D 2 molecules are formed in excited states, 13.5% in the v ′′ = 4, all J ′′ state and 3.5% in the v ′′ = 4, J ′′ = 2 excited state immediately during deposition. Even if this is only a rough estimate, we find that more than 1/3ofD 2 molecules are formed on a silicate surface at 12 K through the H-K mechanism either in an excited or ground state.

Finally, we look at the mechanisms of H 2 formation when the REMPI yield is substantially higher than the QMS-TPD. At steady state, the number of atoms on the surface is dN/dt = 0 = F -WN, where F is the flux, W = v × exp(-E des /kT)isthe desorption rate, v = 10 12 Hz is the atom-surface fundamental vibration frequency and E des is the desorption energy of an atom, taken to be 44 meV from fitting rate equations to H 2 formation data on amorphous silicates (Perets et al. 2007). We get N ∼ 1 × 10 7 cm -2 . The rms distance traveled by an atom is ( R 2 ) 1/2 = 4Dt, where D is the diffusion constant (va 2 /4 exp(-E diff /kT)), a is a typical atom's single jump distance, and E diff = 35 meV is the single atom diffusion energy on an amorphous silicate (Perets et al. 2007). This distance ( R 2 ) 1/2 ∼20 Å at T = 36 K is much smaller than the inter-particle distance, ∼(N) -1/2 = 3 × 10 4 Å. Therefore, the L-H mechanism can contribute only negligibly to H 2 formation at this temperature, and another mechanism has to be invoked. Because of the mass mismatch between the H atom and the solid, the atom can travel quite a distance before it becomes equilibrated with the surface. Theoretical calculations of the distance traveled by a hydrogen atom landing on the surface of amorphous ice, arguably much rougher than the current amorphous silicate film, via the hot-atom mechanism yields ∼10-160 Å depending on the surface temperature and kinetic energy of the particle.

DISCUSSION AND CONCLUSIONS

Although both rate equations derivations (Katz et al. 1999) based on experiments on silicates (Perets et al. 2007;Vidali et al. 2007) or theoretical work on graphite (Cuppen & Hornekaer 2008;Gavardi et al. 2009) show that the window of high efficiency for H 2 formation is extended and shifted to higher dust temperature with respect to crystalline surfaces, the upper range does not go much beyond 18 K. However, there are regions in the ISM where temperatures can be higher, for example, in PDRs. In this work, we find that the detection of rovibrationally hot D 2 formed on the surface of a silicate extends well beyond that range. While the formation of H 2 at low temperature (<30 K) is well explained by the L-H reaction (Morisset et al. 2004(Morisset et al. , 2005)), at higher temperature the residence time of H atoms is too short. At 36 K, the rms diffusion path is 2 nm, but the coverage is negligible since the residence time is ∼10 -7 s. Therefore, the formation of D 2 can only occur if D atoms travel a significant distance to find other atoms that just landed on the surface (the use of D atoms prevents the contribution of atoms from a hypothetically hydrogenated surface). We propose that the H-K mechanism provides an explanation for the formation of D 2 at high surface temperature. From Figure 4 we clearly see that the decay rate of the yield with surface temperature is different for the QMS versus the REMPI signals. Therefore, the new and innovative picture that emerges is that there is simultaneous formation of D 2 via the L-H and the H-K mechanisms. The former leads to formation of molecules that remain on the surface, indicating, as it has been amply observed in experiments (Roser et al. 2003;Hornekaer et al. 2003;Congiu et al. 2009), that the mechanism of transfer of the energy gained in the formation of the molecule is quickly taken by the substrate. The H-K mechanism leads to the formation of rovibrationally hot D 2 and the two coexist in the temperature range explored (5-70 K). In the preliminary analysis presented in the previous section, we find that more than one third of D 2 molecules are formed on a silicate surface at 12 K through the H-K mechanism either in an excited or ground state. However, while above T ∼ 18-28 K the L-H mechanism produces few molecules because of the short residence times, rovibrational hot molecules are still detected (see Figure 4). In our interpretation, it is the H-K mechanism that provides the fast motion that allows atoms to find each other, overcoming the constraint of the short residence time. Calculations should be done (as for H diffusion on amorphous ice (Masuda et al. 1998;Al-Halabi & Van Dishoeck 2007;Matar et al. 2008;Watanabe et al. 2010;Matar et al. 2010)) to check this. Given the recent discovery on how H atoms interact on the surface of graphite (Hornekaer et al. 2006), there might be other mechanisms at play. These results should be of major interest to modelers of the dynamics of infall of an interstellar cloud as well as to theoreticians working on models of the chemical evolution of a cloud, where the inclusion of gas-surface processes is recognized to be important (Herbst 2001;R öllig et al. 2007).

Finally, we comment on the non-detection of H 2 formation pumping in the ISM that our results should explain. A reason for that is that so far no place has been found where H 2 excitation due to formation pumping could be separated in the line of sight from others dominant processes, radiative and collisonal pumping. These processes are searched either from groundbased telescopes in the near infrared and particularly in the K-band (1.8-2.5 µm) or from space in the 5-28 µm range (for transitions between rotational levels in the lower v ′′ = 0 ground state). According to our findings, emission lines starting from v ′ = 4tolowerv ′′ will have to be searched in the X, J,orH (resp. 1/1.25/1.65 µm) infrared bands of atmospheric windows. In the case of diffuse or dark clouds, depending on the stage of star formation, one needs to look for places with temperature in the 5-60 K range and where silicates are still bare. Furthermore, they could be observed only in restricted and narrow regions in molecular clouds, most likely edge-on, because forbidden transitions require large column density in the line of sight. The advent of the new generation of giant telescopes (including adaptive optics) giving access to spatial resolution down to a few milliarcsecond in the spectroscopic mode will lead to the disentangling of this observational puzzle. This is not exclusive of multiwavelength studies with Herschel and ALMA for simultaneous H 2 observations and dust characterization as in the case of NGC6720 (van Hoof et al. 2010).

Funding was provided by ANR (Agence Nationale de la Recherche; contract 07-BLAN-0129), Conseil Régional d'Ile de France (SESAME contract I-07-597R), by Conseil Général du Val d'Oise and through the European Community FP7-ITN-Marie-Curie Programme (grant agreement no. 238258). G.V. is supported by NSF-Astronomy & Astrophysics Division, grant no. 0908108.

Figure 1 .

 1 Figure 1. Schematic of the "FORMOLISM" (FORmation of MOLecules in the InterStellar Medium) apparatus used for these experiments.

Figure 2 .

 2 Figure2. in red, shows the formation of D 2 in the v ′′ = 4, J ′′ = 2 excited state when irradiating the sample at 12 K with D atoms from time t = 0 to 600 s. Blue (scale offset): QMS D 2 signal; its rise at t = 0 is due to the increased gas load from the beam. Orange: sample temperature. During the temperature ramp of the TPD (beam valve closed), a sharp peak of D 2 molecules desorbing from the surface is seen whereas no signal appears on the REMPI (the second broad and smaller peak during the TPD is an artifact due the cryostat, see the text). When the temperature reaches 70 K, the irradiation resumes. The QMS signal goes to the same level as at t = 0 (which is due to the extra D 2 gas load from the beam) while the REMPI signal is almost negligible.

Figure 3 .

 3 Figure3. Same as Figure2with the sample at 36 K during irradiation. There is no TPD-QMS signal of D 2 from the surface during the temperature ramp (t > 360 s) as no deuterium stick on the surface at such temperature.

Figure 4 .

 4 Figure 4. Integrated intensities of the REMPI (D 2 ) signal during irradiation and of the QMS (D 2 ) signal of the TPD after irradiation vs. the sample temperature during irradiation. The HD signal is due to exchange reactions on the sample and/or walls of the apparatus.


Benoît Braïda
email: benoit.braida@upmc.fr

Julien Toulouse

Michel Caffarel

C J Umrigar

Quantum Monte Carlo with Jastrow-Valence-Bond wave functions

We consider the use in quantum Monte Carlo calculations of two types of valence bond wave functions based on strictly localized active orbitals, namely valence bond self-consistent-field (VB-SCF) and breathing-orbital valence bond (BOVB) wave functions. Complemented by a Jastrow factor, these Jastrow-Valence-Bond wave functions are tested by computing the equilibrium well depths of the four diatomic molecules C2, N2, O2, and F2 in both variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC). We show that it is possible to design compact wave functions based on chemical grounds that are capable of describing both static and dynamic electron correlation. These wave functions can be systematically improved by inclusion of valence bond structures corresponding to additional bonding patterns.

I. INTRODUCTION

Quantum Monte Carlo (QMC) methods (see e.g. Refs. [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF][START_REF]Quantum Monte Carlo Methods in Physics and Chemistry[END_REF][3] constitute an alternative to standard quantum chemistry approaches for accurate calculations of the electronic structure of atoms, molecules and solids. The most commonly-used approach consists of optimizing a flexible trial wave function in a variational Monte Carlo (VMC) calculation, and then using the resulting wave function in a more accurate fixed-node diffusion Monte Carlo (DMC) calculation.

For atoms and molecules, the most common form of trial wave function consists of a Jastrow factor, expected to describe a major part of the dynamic electron correlation effects, multiplied by a single Slater determinant of orbitals expanded in a localized one-particle basis. For systems featuring important static correlation effects, the single Slater determinant is usually replaced by a linear combination of several Slater determinants. This has been shown to lead to systematic improvements of both VMC and DMC total energies, provided that the wave function is properly reoptimized [4][5][6]. In particular, fullvalence complete active space (FVCAS) expansions, i.e. including all the Slater determinants that can be generated by distributing all the valence electrons in all the valence orbitals, were found to yield near chemical accuracy for bonding energies of first-row homonuclear diatomic molecules at the DMC level [6].

However, FVCAS expansions involve large numbers of Slater determinants that scale exponentially with system size, so this approach cannot be applied for large systems. Instead, more compact and yet systematic wave functions are desired for QMC calculations. Naively-truncated expansions of Slater determinants of delocalized molecular orbitals cannot be seen as a general solution since they of-ten suffer from lack of consistency for calculating energy differences between two systems or two states. In contrast, pairing-type wave functions should be an adequate, compact and systematic way of describing static correlation in QMC. Indeed, several forms of pairing wave functions have recently been used in QMC, namely antisymmetrized geminal power (AGP) [7][8][9], pfaffians [10,11] and generalized valence bond (GVB) [12], with promising results on atoms and small molecules.

Here, we explore the use in QMC of valence bond (VB) wave functions based on nonorthogonal active orbitals that are localized on a single atom (for a review of VB theory, see Refs. [START_REF] Mcweeny | Methods of Molecular Quantum Mechanics[END_REF][START_REF] Gallup | Valence Bond Methods: Theory and applications[END_REF][START_REF] Shaik | A Chemist's Guide to Valence Bond Theory[END_REF]. The merits of valence bond self-consistent field (VBSCF) wave functions [START_REF] Van Lenthe | [END_REF], are tested, along with more flexible breathing-orbital valence bond (BOVB) wave functions [17,18]. These approaches can describe static correlation effects with compact wave functions containing only a limited number of VB structures (i.e., linear combinations of Slater determinants made of spin-singlet pairs of orbitals) which are selected on chemical grounds. In comparison to the GVB method, which is usually considered in the perfect-pairing and strong-orthogonality approximations, the VBSCF and BOVB approaches are more flexible and provide a greater chemical interpretability thanks to the use of localized nonorthogonal orbitals. In contrast to traditional analytical calculations, nonorthogonal orbitals do not lead to extra computational cost in QMC. On the contrary, the use of localized nonorthogonal orbitals could potentially reduce the computational cost as compared with a delocalized molecular-orbital approach in three ways: a) localized orbitals are negligibly small and need not be computed beyond some cutoff distance, b) the number of determinants needed may be reduced by having nonorthogonal orbitals, and, c) in the absence of spatial symmetry, the number of orbital parameters to be optimized is smaller for localized orbitals.

BOVB wave functions improve upon VBSCF wave functions by allowing a different set of orbitals in each VB structure, effectively introducing dynamic correla-tion into the Slater determinants, but in a different way than through backflow coordinate transformations [19] or orbital-attached multi-Jastrow factors [20]. In VB practice, it is known that bonding energies at the BOVB level are significantly improved as compared with the VBSCF level, especially in cases where dynamic correlation is important. In QMC, this form of wave functions has already successfully been tested on acetylene dissociation energies [21]. However, it is not known whether the breathing-orbital relaxation yields significant improvement when the wave function is fully optimized in QMC, because it is conceivable that the Jastrow factor in QMC provides much of the same variational freedom.

The paper is organized as follows. In Section II, we offer a brief review of the VBSCF and BOVB approaches, and describe our corresponding Jastrow-Valence-Bond wave functions used in QMC. In Section III, we present results for the four diatomic molecules C 2 , N 2 , O 2 , and F 2 , which span different types of bonding (single, multiple, and three-electron bond) and electron correlation (strong dynamic correlation for O 2 and F 2 , and strong static correlation for C 2 and N 2 ). We discuss the type of electron correlation to include in the wave functions to obtain accurate equilibrium well depths, as well as the importance of orbital optimization in QMC, and the computational cost of our wave functions. Section VI summarizes our conclusions.

II. METHODOLOGY

A. Valence-Bond wave functions

The different valence bond methods can be classified into two main families, depending of the degree of localization of the active orbitals: semilocalized orbitals and strictly localized orbitals [START_REF] Shaik | A Chemist's Guide to Valence Bond Theory[END_REF]22]. In the first family, the optimized active orbitals are expanded on all the basis functions of the entire molecule, and are thus free to fully delocalize, although most of the time they appear to be fairly localized. The GVB method [START_REF] Bobrowicz | Methods of Electronic Structure Theory[END_REF]24] and the spin-coupled valence bond (SCVB) method [25] belong to this family. In the second family, each active orbital is expanded in basis functions centered on a single atom or the atoms of a fragment of the molecule (alternative definition of strictly localized orbitals exists, Ref. 26). The VBSCF [START_REF] Van Lenthe | [END_REF]27] and BOVB [17,18] methods belong to this family. The semilocalized approaches usually provide more compact wave functions, but the delocalization tails of the orbitals may hinder their chemical readability. The strictly localized approaches usually involve more VB structures, but permit the distinction between covalent and ionic contributions to the bonds, thus providing more insight into the nature of the chemical bonds. The interpretative power of the strictly localized VB approaches is hence greater, and allows, among other examples, the characterization of charge-shift bonding [28][29][30], three-electron bonds [31][32][33] or hypervalent bond-ing [33][34][35]. Further, most applications involving VB curve-crossing diagrams [START_REF] Shaik | A Chemist's Guide to Valence Bond Theory[END_REF]36] are within the strictly localized VB formalism.

Any VBSCF wave function can be written as a linear combination,

|Ψ VBSCF = I c I |Φ I , (1)
where |Φ I are spin-adapted VB structures, each one being determined by a choice of an orbital configuration (or orbital occupation) and a spin coupling of these orbitals. We consider VB structures of the form (disregarding normalization)

|Φ I = inactive p â † p↑ â † p↓ active pairs (ij) â † i↑ â † j↓ -â † i↓ â † j↑ × active unpaired q â † q↑ |vac , (2)
where â † pσ (σ =↑, ↓) is a spin-orbital creation operator and |vac is the vacuum state of second quantization. The VB structures are thus made of inactive (always closed-shell) orbitals p, spin-singlet pairs of active orbitals (ij), and possibly remaining unpaired spin-up active orbitals q. We use inactive orbitals that are either localized (expanded on the basis functions centered on a single atom), e.g. for core orbitals, or delocalized (expanded on all the basis functions of all the atoms), e.g. for bonds made of inactive orbitals that do not mix with the active orbitals. We use active orbitals that are always localized on a single atom, and they are typically identified with valence atomic (hybrid) orbitals. Note that Eq. ( 2) encompasses the case of spin-singlet pairing of an active orbital with itself, i.e. i = j giving simply â †

i↑ â † i↓ -â † i↓ â † i↑ = 2â † i↑ â † i↓ .
The spin-coupling scheme based on singlet pairing used in Eq. ( 2) is usually referred to as the Heitler-London-Slater-Pauling (HLSP) scheme. There exist other spin-coupling schemes, but the HLSP scheme has the advantage of providing a clear correspondence between each VB structure and a Lewis chemical structure, two singlet-paired active orbitals representing either a bond or a lone pair. In principle, considering all possible pairings exhausts, for a given orbital configuration, all the spin eigenstates (or spin couplings) of fixed quantum numbers S = N unpaired /2 and M S = +S (N unpaired is the number of spin-up unpaired electrons). In fact, considering all possible pairings leads to an overcomplete set of spin couplings, but they can be reduced to a complete basis of (non-redundant) spin couplings, also called Rumer basis [START_REF] Shaik | A Chemist's Guide to Valence Bond Theory[END_REF]37,[START_REF] Pauncz | Spin Eigenfunctions[END_REF]. In practice, for most practical applications, only a small number of chemically relevant VB structures are kept in the calculation.

In practical VBSCF calculations, each VB structure is expanded in 2 Npairs Slater determinants, |Φ I = µ d I,µ |D µ , where N pairs is the number of pairs of different active orbitals (i = j). The coefficients of the determinants d I,µ for a given VB structure are all equal in absolute value. The VB structure coefficients and orbital coefficients on the basis are then optimized using a direct generalization of the usual multi-configuration self-consistent-field (MCSCF) algorithm for arbitrary nonorthogonal orbitals [27], or algorithms more specific to VB theory [START_REF] Wu | [END_REF]40]. In fact, this VBSCF procedure is general enough to also permit optimization of GVB and SCVB wave functions if the active orbitals are allowed to delocalize. Note that GVB wave functions, in the most commonly used perfect-pairing approximation, are made of only one VB structure, while the SCVB wave functions include all possible spin couplings.

As an illustration, consider the single-bonded Li 2 molecule. A typical VBSCF wave function is

|Ψ VBSCF = c 1 |Φ cov + c 2 |Φ ionA + c 3 |Φ ionB , (3)
with a covalent VB structure obtained by spin pairing the 2s orbital of the first Li (A) with the 2s orbital of the second Li (B)

|Φ cov = â † 1sA↑ â † 1sA↓ â † 1sB ↑ â † 1sB ↓ × â † 2sA↑ â † 2sB ↓ -â † 2sA↓ â † 2sB ↑ |vac , (4)
and two ionic VB structures obtained by either spin pairing 2s A with itself,

|Φ ionA = â † 1sA↑ â † 1sA↓ â † 1sB ↑ â † 1sB ↓ â † 2sA↑ â † 2sA↓ |vac , (5) or
, symmetrically, by spin pairing 2s B with itself,

|Φ ionB = â † 1sA↑ â † 1sA↓ â † 1sB ↑ â † 1sB ↓ â † 2sB ↑ â † 2sB ↓ |vac , (6)
the core inactive orbitals 1s A and 1s B remaining always doubly occupied. These three VB structures correspond respectively to three Lewis chemical structures, Li-Li, -Li Li + , and + Li Li -, respectively. Obviously, in this case, spatial symmetry implies that c 2 = c 3 and the two ionic VB structures can thus be combined in a single "symmetry-adapted structure", which we refer to as a configuration state function (CSF). In this simple case, the VBSCF wave function is essentially a complete active space (CAS) wave function with the two valence electrons distributed in the valence bonding σ g and antibonding σ u molecular orbitals, i.e. CAS(2,2), and would be strictly identical to it in a minimal basis set. In less trivial cases, the selected VB structures correspond only to a subspace of a CAS. However, if structure selection is correctly made, the VBSCF wave function can lead to observable predictions in close agreement to a CAS wave function. Indeed, in addition to the most chemically significant structures, the CAS wave function includes a lot of improbable covalent and (multi-)ionic structures that are of little importance. Thus, a VB wave function can represent an intelligent way of truncating a CAS expansion based on chemical grounds.

We adopt the usual definitions of electron correlations used in quantum chemistry. Static correlation corresponds to the correlation included in a FVCAS wave function as compared with a Hartree-Fock reference, and dynamic correlation corresponds to the remaining correlation resulting from all the outer-valence excitations in a full configuration interaction expansion. The VBSCF method handles static correlation only. Dynamic correlation can be added in the usual way by performing configuration interaction on top of VBSCF [41,42], referred to as VBCI, or through a perturbative treatment [42,43]. The BOVB method is an alternative way of incorporating some dynamic correlation, keeping the wave function as compact as the VBSCF one. This method consists in allowing a different set of (inactive and active) orbitals in each VB structure. The orbitals can thus adapt to the different charge distributions in the VB structures. For example, for the Li 2 wave function of Eqs. ( 3)-( 6), in BOVB the 2s orbitals can become more diffuse in the ionic structures than in the covalent structure. This extra degree of freedom is seen as introducing dynamic correlation, since it roughly corresponds to adding excitations to atomic orbitals of higher principal quantum numbers. Just as in the VBSCF method, inactive valence orbitals can be chosen to be either localized or delocalized in the BOVB method. Finally, in the ionic structures, a doublyoccupied active orbital (the 2s orbital for Li 2 ) is sometimes split into two independent orbitals that are spinsinglet coupled, which introduces additional dynamic correlation. The most flexible BOVB method is thus obtained by both splitting of active orbitals in the ionic structures and delocalization of inactive orbitals, which is referred to as split-delocalized-BOVB (SD-BOVB) in the VB literature.

B. Jastrow-Valence-Bond wave functions

We discuss now the parametrization of the two types of Jastrow-Valence-Bond wave functions, J-VBSCF and J-BOVB, that we use in QMC. The J-VBSCF wave function is parameterized as (see Ref. 6)

|Ψ J-VBSCF = Ĵ e κ I c I |Φ I , (7)
where Ĵ is a Jastrow operator, e κ is an orbital "rotation" operator, and |Φ I are CSFs, i.e. VB structures or linear combination of a few VB structures to accommodate spatial symmetry. Each CSF is a short linear combination of products of spin-up and spin-down Slater determinants, composed generally of nonorthogonal orbitals expanded on Slater basis functions. The orbital operator κ is κ = kl κ kl Êkl where κ kl are the orbital "rotation" parameters and Êkl is the singlet excitation operator from orbital l to orbital k, Êkl = â † k↑ bl↑ + â † k↓ bl↓ , written with dual biorthogonal orbital creation and annihilation operators â † kσ and blσ (see, e.g., Ref. 44). The parameters to optimize are the parameters in the Jastrow factor, the CSF coefficients, the orbital rotation parameters and the exponents of the basis functions.

The J-BOVB wave function is similar except that each CSF is allowed to have a different set of orbitals, which can be written as

|Ψ J-BOVB = Ĵ I c I e κI |Φ I , ( 8
)
where κI is the orbital operator for the set of orbitals in the CSF |Φ I . In practice, Eq. ( 8) can be recast in the form of Eq. ( 7) by defining the orbital operator κ to operate on all sets of orbitals: κ = I κI . We use a Jastrow factor consisting of the exponential of the sum of electron-nucleus, electron-electron, and electron-electron-nucleus terms, written as systematic polynomial and Padé expansions [START_REF] Umrigar | [END_REF] (see also Refs. 46 and 47).

C. Computational details

We start by generating a standard restricted Hartree-Fock (RHF) wave function with the program Gaussian [START_REF] Frisch | Gaussian 03, Revision C.02[END_REF] and a VBSCF or BOVB wave function with the program XMVB [START_REF] Song | [END_REF]. We use the core-valence triplezeta quality Slater basis (CVB1) of Ema et al. [50]. For the XMVB calculations, each Slater function is expanded into 15 (core basis functions), 9 (s and p valence basis functions), or 6 (d polarization functions) Gaussian functions. We use here 5 d spherical functions while 6 d Cartesian functions were used in Ref. 6. The standard RHF, VBSCF or BOVB wave function is then multiplied by our Jastrow factor, and QMC calculations are performed with the program CHAMP [START_REF] Champ | a quantum Monte Carlo program written[END_REF] using the true Slater basis set rather than its Gaussian expansion.

The Jastrow, CSF, orbital and basis exponent parameters are simultaneously optimized by minimizing the energy plus a small percentage of the energy variance (1%) with the linear optimization method [4][5][6] in VMC, using an accelerated Metropolis algorithm [START_REF] Umrigar | [END_REF][START_REF] Umrigar | Quantum Monte Carlo Methods in Physics and Chemistry[END_REF]. Once the trial wave function has been optimized, we perform a DMC calculation within the short-time and fixed-node (FN) approximations (see, e.g., Refs. [START_REF] Grimm | [END_REF][55][56][57][58]. We use an imaginary time step of τ = 0.01 hartree -1 in an efficient DMC algorithm with very small time-step errors [59].

For orbital optimization, the orbitals are partitioned into three sets: inactive (doubly occupied in all determinants), active (occupied in some determinants and unoccupied in others), and virtual (unoccupied in all determinants). All inactive (including 1s core) and active orbitals are optimized. The non-redundant excitations to consider are inactive → active, inactive → virtual, active → virtual and active → active. If the action of the excitation Êkl on the wave function is not zero but the reverse excitation Êlk is zero, then the "orthogonality" condition κ lk = -κ kl is imposed. For some active-active excitations, both direct and reverse excitations ( Êkl and Êlk ) may be allowed, and thus it makes sense for localized orbitals to remove the "orthogonality" constraint by treating κ kl and κ lk as independent parameters. This results in only a very few (if any at all) additional orbital parameters for the wave functions considered here. We note that, when considering active-active excitations, redundancies between two orbital wave-function derivatives or between an orbital wave-function derivative and a CSF frequently occur, and must be detected and eliminated. Localized orbitals do not have the point group symmetry of the system, so the number of orbital excitations cannot be reduced based on the non mixing of irreducible representations, as usually done. However, the number of orbital excitations is greatly reduced by forbidding mixing between orbitals of different localization classes, i.e. expanded on different subsets of basis functions. For BOVB wave functions, one has a different set of (occupied and unoccupied) orbitals for each CSF [60] and only orbitals that belong to the same set are allowed to mix.

III. RESULTS

The Jastrow-Valence-Bond wave functions are tested on the ground states of the molecules C

2 ( 1 Σ + g ), N 2 ( 1 Σ + g ), O 2 ( 3 Σ - g )
, and F 2 ( 1 Σ + g ) at their experimental bond lengths of 1.2425 Å, 1.0977 Å, 1.2075 Å, and 1.4119 Å, respectively. The core inactive 1s orbitals are always taken as localized. For O 2 and F 2 , we use delocalized inactive π orbitals both at the VBSCF and BOVB levels. This reduces the Pauli repulsion between the otherwise localized lone pairs. The precise composition of all wave functions is explained below. In addition to total energies, we present molecular well depths all computed with respect to the separated atoms calculated at the FVCAS level. For the N 2 , O 2 , and F 2 molecules, both VBSCF and BOVB wave functions in the limit of a complete basis of structures (Rumer basis) dissociate into the FVCAS atomic limit, which simply corresponds for these atoms to single-determinant wave functions. For C 2 , VBSCF in the limit of a complete basis of structures also dissociates into the FVCAS atomic limit, which corresponds for the carbon atom in its triplet ground state to a CAS(2,2) wave function, since the carbon lone pair can doubly occupy either a 2s or a 2p orbital. With this definition, the well depths actually allow a comparison of the absolute quality of the different molecular wave functions at equilibrium since the atomic reference is the same.

A. Single bonding-pattern wave functions

We first consider wave functions describing what we call a single bonding pattern. A bonding pattern is a compact Lewis-style picture of the molecular electronic structure, but which is made of several (covalent and ionic) VB structures. The dominant bonding patterns of C 2 , N 2 , O 2 , and F 2 are depicted in Fig. 1. The conventions of the drawings are the following. Each single line connecting two atoms represents a σ or π two-electron bond. In our wave functions, each two-electron bond is generally composed of one covalent and two ionic components. Hence, F 2 (single bond) is described by a wave function containing 3 VB structures, whereas N 2 (triple bond) leads to a wave function with 3 3 = 27 VB structures. Each pair of dots on a single atom represents a lone pair, which is described by two opposite-spin electrons occupying the same orbital.

The single dots connected by a curve which appears for C 2 represent a (formally covalent) singlet coupling between two electrons occupying two different orbitals, one located on each carbon in that case. The two corresponding coupled orbitals are sp hybrids pointing in opposite direction, so that they overlap and thus interact very weakly. Hence, this corresponds to a singlet biradical state, the coupling between the two single electrons being only formal to ensure the proper spin symmetry. Very recent VB calculations by Su et al. [62] have shown this biradical triple-bond situation to be the dominant bonding pattern for C 2 at equilibrium distance. In order to keep the number of VB structures low, only the single-ionic structures, as well as the double-ionic structures with opposite charges on each carbon, are included in our wave function. This leads to a wave function with 21 VB structures.

The O 2 molecule is considered in its spin-triplet (paramagnetic) ground state. The picture in Fig. 1 corresponds to the textbook molecular-orbital diagram for this molecule, displaying a σ two-electron bond (single line connecting the two oxygens), and two π three-electron half-bonds, each one represented by a three-dot symbol. In our wave functions, the σ bond is composed of one co-valent and two ionic components, and each three-electron bond is composed of the two resonating structures typical of this kind of bond [31,32]. This leads to a wave function with 12 VB structures.

In each wave function, structures equivalent by symmetry are further grouped together into CSFs. For each dimer, the VBSCF and BOVB wave functions have exactly the same number of CSFs, the only difference being that different sets of orbitals are used for the different CSFs of the BOVB wave functions.

A rough analysis of the relative contributions of the different types of electron correlation to the equilibrium well depths for the four molecules is given in Table I. This analysis leads us to consider separately the case of the O 2 and F 2 molecules which have strong dynamic correlation, and the case of the C 2 and N 2 molecules which have strong static correlation.

Systems with strong dynamic correlation: O2 and F2

The total energies and well depths of O 2 and F 2 computed with standard VB methods and QMC methods are displayed in Table II. We start by analyzing the standard VB results. The difference between the RHF and VBSCF energies can be taken as the definition of the "left-right" part of static correlation. In VB theory, the left-right correlation is viewed as the energy gained by restoration of the correct balance between covalent and ionic contributions to the bond, starting from the incorrect 50% covalent -50% ionic RHF description. The contributions of left-right correlation to the well depths appear to be very large here, 59.3 and 45.1 kcal/mol for O 2 and F 2 , respectively. A significant part of dynamic correlation is retrieved by going from VBSCF to BOVB, improving well depths by 33.2 and 18.9 kcal/mol for O 2 and F 2 , respectively. Another 3.9 kcal/mol improvement is obtained for F 2 by splitting the doubly-occupied active orbital in the ionic structures, which brings radial correlation to the active electron pair. The agreement between the BOVB and exact well depths is quite reasonable considering the compactness of these wave functions: only 12 VB structures (combined into 5 CSFs) for O 2 and 3 VB structures (combined into 2 CSFs) for F 2 .

We now discuss the VMC results. The Jastrow factor keeps electrons apart and describes part of dynamic and static correlation. It is a particularly efficient way of including dynamic correlation, that otherwise requires in configuration interaction calculations a very large number of highly excited determinants. It also includes some static left-right correlation since it greatly reduces the probability of ionic configurations [6]. As expected, a sharp decrease of the total energies is observed in VMC, indicating that a large amount of the total correlation energy is recovered by the Jastrow factor. However, for well depths, and most chemical properties, we are instead interested in a balanced description of the differential correlation energy in the molecule and the separated TABLE II. Total energies (in hartree) and well depths (in kcal/mol) for the ground states of O2 and F2 at their experimental bond lengths calculated in RHF, VBSCF and BOVB, and in VMC and DMC using corresponding Jastrow-Valence-Bond wave functions with all the parameters (Jastrow, CSF coefficients, orbitals and basis exponents) fully optimized in VMC. For comparison, DMC results for FVCAS wave functions from a previous study are also shown. The J-VBSCF VMC well depths remain in relatively poor agreement with the exact values, underestimated by 13.9 and 14.3 kcal/mol for O 2 and F 2 , respectively. It is expected that these errors are mainly due to some missing dynamic correlation, rather than static correlation corresponding to some missing VB structures. This point will be checked in Section III.B. Additional correlation can be introduced in VMC by replacing the VBSCF determinant expansion by the more flexible BOVB form. By using a different set of orbitals in each VB structure, the BOVB wave function allows the orbitals to be more diffuse on the negatively charged atom, and more compact on the positively charged atom of the ionic structures, compared to the neutral covalent structure. This "breathing-orbital effect" can be seen as a manifestation of dynamic correlation. It is clear that such breathing-orbital effects are partly redundant with the correlations introduced by the Jastrow factor. In practice, it is observed that with the limited form of the Jastrow factor employed, the J-BOVB wave functions improve a bit upon J-VBSCF wave functions: the VMC well depths are bettered by 5.1 and 2.5 kcal/mol for O 2 and F 2 , respectively. As expected, these improvements are much smaller than those obtained when going from standard VBSCF to BOVB, showing that a large part of dynamic correlation is already described by the Jastrow factor. The agreement between the VMC J-BOVB well depths and the exact values is still not good, with an underestimation of 8.8 kcal/mol for O 2 , and 10.7 kcal/mol (8.4 kcal/mol with BOVB-split) for F 2 , respectively. Even more disappointing, the VMC J-BOVB well depth is actually worse than the standard BOVB one for O 2 , and the same is true for BOVB-split for F 2 . This indicates that our Jastrow factor describes correlation effects in the atom more accurately than in the (less symmetric) diatomic molecule. By incorporating more correlation in the separated atoms than in the molecule, the Jastrow factor thus makes the BOVB well depths less accurate. In fact, it is significant that for F 2 splitting the ionic electron pairs improves the VMC well depth by 2.3 kcal/mol, indicating that this radial correlation effect in the molecule was not well described by the Jastrow factor. The more flexible multi-Jastrow approach of Ref. 20 could solve this problem by better handling the delicate balance between atomic and molecular correlation effects [63].

Let us now consider the fixed-node DMC results. The total energies and the well depths reflect the quality of the nodes of the different wave functions. Both total energies and well depths are much improved compared to the VMC results. J-RHF wave functions give DMC well depths of moderate accuracy, underestimated by 6.2 and 9.2 kcal/mol for O 2 and F 2 , respectively. J-VBSCF wave functions improves DMC well depths but by rather limited amounts, 1.8 and 4.2 kcal/mol for O 2 and F 2 , respectively, showing that most of left-right correlation was already recovered with the J-RHF nodes. Using J-BOVB wave functions allow further gains in accuracy, leading to results close to the ones obtained with much less compact J-FVCAS wave functions. The best DMC J-BOVB calculations underestimate the well depths by only 2.4 and 1.1 kcal/mol for O 2 and F 2 , respectively.

Systems with strong static correlation: C2 and N2

Let us now consider C 2 and N 2 , two molecules having strong static correlation even at equilibrium distance. RHF is expected to be particularly poor for these molecules, because of the triple bond nature of these systems. In a RHF wave function, each two-electron bond is constrained to be 50% covalent -50% ionic, and in addition all (multi-)ionic forms are constrained to have the same weights, which is unphysical. In particular, bi-ionic structures with two positive charges on one atom and two negative charges on the other will have the same weight as neutral bi-ionic structures (one positive and one negative charge on each atom). Similarly, tri-ionic structures with three (positive or negative) charges of the same sign on each atom will be forced to have the same weight as the mono-charged-atom tri-ionic structures. Note that GVB wave functions also have these constraints on the weights of multi-ionic structures, this form of VB wave function being only able to take into account intrapair left-right correlation (covalent vs. ionic balance within only one two-electron bond) but not interpair correlation. On the other hand, a VBSCF wave function can include both intrapair and interpair left-right correlation. More importantly, for C 2 , the singlet biradical character is expected to generate additional static correlation beyond left-right correlation, which will require extra bonding patterns.

The total energies and well depths of these molecules are displayed in Table III. As expected, RHF strongly underestimates the well depths (note that the negative well depth obtained for C 2 is in fact due to the choice of the FVCAS reference for the separate atoms). VBSCF appears to retrieve a large part of the missing correlation, as it improves the well depths by 120.6 and 78.8 kcal/mol for C 2 and N 2 , respectively, as compared with RHF values. The VBSCF well depths are however still underestimated by more than 30 kcal/mol for both molecules. It was impossible to carry out standard BOVB calculations to assess dynamic correlation effects at the VB level, as the number of VB structures (21 for C 2 and 27 for N 2 ) is too large to allow such calculations. Last, we checked that GVB (not shown in Table III) greatly underestimates the well depths as expected, giving 84.0 and 164.9 kcal/mol for C 2 and N 2 , respectively, values which are about 30 kcal/mol worse than those from VBSCF.

As for O 2 and F 2 , addition of a Jastrow factor in VMC considerably improves the results. This improvement is dramatic for single-determinant wave functions, the J-RHF well depths being bettered by 110.7 and 88.4 kcal/mol for C 2 and N 2 , respectively, as compared with the RHF ones. The improvement is much less important but still significant for the VBSCF wave functions, with J-VBSCF well depths improved by 21.6 and 22.6 kcal/mol for C 2 and N 2 , respectively, compared to the VBSCF ones. As for O 2 and F 2 , these findings show that the Jastrow factor can retrieve a significant amount of left-right correlation. The accuracy of the VMC J-VBSCF calculations is similar to what has been found for O 2 and F 2 , with well depths underestimated by 15.1 and 11.3 kcal/mol for C 2 and N 2 , respectively.

The DMC results follow the trends found in VMC, but the differences between C 2 and N 2 show up more clearly. The N 2 well depth is underestimated by 7.3 kcal/mol with a J-RHF wave function, which is actually in the range of what was observed for O 2 and F 2 . A J-VBSCF wave function leads to a DMC well depth close to chemical accuracy, indicating that the quality of the J-VBSCF nodes are good enough get good energy differences. The situation is quite different for C 2 . The J-RHF DMC well depth is underestimated by 26.7 kcal/mol, far from the accuracy reached for the other molecules. Using a J-VBSCF wave function greatly improves the DMC well depth, but with an error of 8.6 kcal/mol, it is still far from chemical accuracy. For O 2 and F 2 , a further improvement is possible by using J-BOVB wave functions, but for C 2 even if it were possible to use a J-BOVB wave function a significant improvement is not expected. Because of the multiradical character of C 2 , we expect the missing correlation to still be of static nature, rather than dynamic. The route for improvement is to go beyond a single bonding-pattern description.

B. Multiple bonding-pattern wave functions

A straightforward way of improving our J-VBSCF wave functions, particularly for the C 2 molecule, is to increase the number of VB structures included in the wave function. Adding VB structures should be done in a systematic but selective way, as the total number of non-redundant VB structures (for all orbital configurations and spin couplings) can be very large, e.g. up to 1764 structures for C 2 (for the formula giving the number of VB structures in a complete (Rumer) basis, see Ref. [START_REF] Weyl | Theory of Groups and Quantum Mechanics[END_REF]. A route to extend the VB description beyond a single bonding pattern is to systematically derive possible alternative bonding patterns, and classify them in terms of stability following Lewis-type rules. In particular, the most bonded patterns will be preferred over the less bonded ones, a σ bond will be preferred over a π one, and a three-electron (or one-electron) bond will be counted as half a bond. Following these rules of thumb, TABLE III. Total energies (in hartree) and well depths (in kcal/mol) for the ground states of C2 and N2 at their experimental bond lengths calculated in RHF and VBSCF, and in VMC and DMC using corresponding Jastrow-Valence-Bond wave functions with all the parameters (Jastrow, CSF coefficients, orbitals and basis exponents) fully optimized in VMC. For comparison, DMC results for FVCAS wave functions from a previous study are also shown. a description extended to multiple bonding patterns is proposed in Fig. 2 for the C 2 , N 2 , O 2 , and F 2 molecules. The main bonding pattern for C 2 is of acetylenic type, displaying one σ and two π bonds, as well as two singletcoupled (biradical) electrons. It is natural to postulate for this molecule a secondary bonding pattern of ethylenic type, with one σ and one π bond, thus leaving one pair of singlet-coupled electrons on each atom. Another doubly bonded pattern could be postulated for C 2 , with two π bonds. However, by applying our simple rules, it is expected to be less stable than the ethylenic pattern, as σ bonds are usually more stable than π ones. In a recent VB study, this π-doubly bonded pattern has indeed been found to be of very little importance for the ground state at equilibrium distance [62], so it will not be considered here. Considering a C 2 description with two bonding patterns (acetylenic and ethylenic), and including the two covalent and all mono-ionic and neutral bi-ionic structures, leads to a set of 35 VB structures, which can be combined into only 12 CSFs when all the structures that must have equal weights by symmetry are combined. All calculated well depths are gathered in Table IV. The first three lines just repeat, for ease of comparison, the well depths obtained using single bonding-pattern (J-)VBSCF wave functions, previously shown in Tables II and III. The next three lines display the well depths obtained using multiple bonding-patterns (J-)VBSCF wave functions, and the last three lines some reference data. At the standard VBSCF level already, the well depth of C 2 is significantly improved, by 15.4 kcal/mol, when going from the single bonding-pattern to the multiple bondingpattern wave function. Interestingly, a large part of this gain is retained at the VMC J-VBSCF level, with a 9 kcal/mol improvement. This suggests that the Jastrow function is not able to significantly recover the part of static correlation that is included when going from single bonding pattern to multiple bonding patterns. This observation partly transfers to the DMC results as well, with a significant 4.9 kcal/mol improvement of the well depth. In conclusion, the multiple bonding-pattern J-VBSCF wave function for C 2 gives both VMC and DMC well depths in reasonable agreement with estimated exact values (errors of 6.1 and 3.7 kcal/mol, respectively), while still being a very compact wave function (only 12 CSFs compared to 165 CSFs for the FVCAS wave function). The remaining error could be progressively eliminated by adding VB structures corresponding to additional bonding patterns, starting with the pattern with two π bonds and no σ bond, and then including patterns with only one (σ or π) bond. It is expected that a J-VBSCF wave function including all the main VB structures belonging to these bonding patterns would give almost exact dissociation energies, as has recently been shown to be the case at the VBCI level [62], while still having a significantly more compact wave function than the FVCAS one.

Following the same reasoning as for C 2 , a secondary bonding pattern displaying one σ and one π bond can be postulated for the N 2 molecule, corresponding to picture b in Fig. 2. This bonding pattern could also be obtained by starting from the common Lewis structure of spinsinglet O 2 , and doubly ionizing one of the two lone pairs. For N 2 , a π-doubly bonded pattern is likely to be of even less importance than for C 2 , as this would result in an N atom with a completely filled σ space, i.e. two (highly repulsive) σ lone pairs facing each other, so we rule it out. Last, singly bonded patterns with either a σ or a π bond, corresponding to picture c, are also included, by analogy with the secondary bonding pattern of C 2 . For all these extra patterns, the covalent and all the neutral and mono-charged-atom ionic VB structures are included. This gives 27 extra VB structures, and so a total of 54 structures (with the main bonding pattern) which can be combined into 24 CSFs. This multiple bondingpattern wave function improves the VBSCF well depth of N 2 by 15.4 kcal/mol in comparison to the single bondingpattern calculation, while still being a substantially more compact wave function than the FVCAS one (24 vs. 107 CSFs). The improvements in VMC and DMC are less TABLE IV. Well depths (in kcal/mol) for the ground states of C2, N2, O2 and F2 at their experimental bond lengths calculated in VBSCF, VMC J-VBSCF and DMC J-VBSCF using single and multiple bonding-pattern wave functions. All the parameters of the J-VBSCF wave functions (Jastrow, CSF coefficients, orbitals and basis exponents) have been fully optimized in VMC. For comparison, results for FVCAS wave functions from a previous study are also shown. pronounced, though not negligible, with a 2-3 kcal/mol gain. The DMC J-VBSCF well depth is very close to the estimated exact value, indicating the good quality of the nodes of the multiple bonding-pattern wave function.

For O 2 , the main bonding pattern corresponds to a double-bond situation, with a σ bond and two π half bonds. There are not many possibilities for deriving alternative patterns that are compatible with the spatial and spin-triplet symmetry of the ground state, and in particular it is impossible to preserve the σ bond. Pat-tern b keeps the two three-electron bonds, while pattern c has one three-electron and one one-electron π half-bonds and a completely filled σ on each atom. There are a total number of 12 VB structures emerging from patterns b and c, and when added to the previous set we arrive at a total number of 24 structures (10 CSFs).

Obviously, even if these extra patterns can be associated with a formal bond order of one, they are not expected to be of much importance in the ground state. This expectation is already confirmed at the standard VBSCF level, as the well depth is improved by only 2.4 kcal/mol, much less than what was observed for C 2 and N 2 . Quite logically, the gains are even smaller in VMC and DMC, with about 1 kcal/mol improvement only. It thus seems that for O 2 the single bonding-pattern wave function already contains most of the static correlation, and that further gain in accuracy requires instead improvement of the description of dynamic correlation, as previously found when going from J-VBSCF to J-BOVB wave functions.

For F 2 , alternatives to the usual σ bond situation are limited too. A π bond situation can be postulated, however it comes with a (highly repulsive) completely filled σ space on each atom. A situation with a σ three-electron bond can be envisaged (picture b) as it corresponds to the (quite stable) ground state of the difluorine anion, but it comes here with a singly-ionized σ lone pair. Considering all these bonding patterns generates a set of 13 VB structures (5 CSFs). As expected, the multiple bondingpattern wave function gives very small improvements. As for O 2 , a much more significant improvement was previously obtained for F 2 by using a BOVB wave function.

C. Effect of orbital optimization

Although for some systems reoptimizing the orbitals in QMC in the presence of the Jastrow factor leads to only small improvements, it is known that it can be very important for other systems, e.g. for calculating certain excitation energies [START_REF] Schautz | [END_REF][66][67]. To examine the effect of orbital optimization, we show in Table V the VMC and DMC well depths of the four molecules for J-RHF and single bonding-pattern J-VBSCF wave functions with or without VMC optimization of the orbital coefficients and the basis exponents. In each case, of course, the same level of optimization has been used for the atoms. At the VMC level, optimization of the orbital coefficients and the basis exponents improves the well depths by the order of 10 kcal/mol for J-RHF and J-VBSCF wave functions. In DMC, the gain from optimization of the orbital coefficients and the basis exponents is smaller but still significant, ranging from about 2 to 7 kcal/mol. This shows that orbital optimization in QMC has a significant impact of the nodes of the wave functions. These improvements in VMC and DMC are primarily due to the optimization of the orbital coefficients, the gains from reoptimization of the basis exponents being generally smaller, TABLE V. Well depths (in kcal/mol) for the ground states of C2, N2, O2 and F2 at their experimental bond lengths calculated in VMC and DMC using single bonding-pattern J-VBSCF wave functions with partial or complete VMC optimization of the parameters: Jastrow (J), CSF coefficients (c), orbital coefficients (o) and basis exponents (e). except in the case of O 2 .

D. Computational cost

Finally, we discuss the computational cost of our Jastrow-Valence-Bond wave functions. In QMC calculations, using multideterminant wave functions, the major bottleneck is the evaluation of Slater determinants along with their first-and second-order derivatives with respect to electron coordinates, so that the computational cost is usually directly determined by the number of these determinants. In addition, orbital optimization in VMC using a super-configuration-interaction-type algorithm [4,6] requires evaluation of numerous singly excited determinants, which is often the bottleneck of the optimization calculations. It is thus desirable to have compact wave functions with small numbers of determinants. In practice, each determinant is actually broken into the product of a spin-up and a spin-down determinant, and only the list of unique spin-up and spin-down determinants is computed. In Table VI, we have reported the numbers of CSFs, numbers of unique spin-up + spin-down determinants, and the numbers of unique spin-up + spin-down singly excited determinants for the VBSCF and BOVB wave functions employed in this work. It is seen that there are about 4-7 times less determinants to evaluate for the single bonding-pattern VBSCF wave functions, in comparison to the FVCAS wave functions. For orbital optimization, these single bonding-pattern VBSCF wave functions also generate about 3-5 times less singly excited determinants. On the other hand, single bonding-pattern BOVB wave functions do not appear as computationally advantageous in comparison to FVCAS wave functions. For O 2 , the BOVB wave function has only a few less determinants, and for F 2 the BOVB wave function in fact has more determinants than the FVCAS wave function. Multiple bonding-pattern VBSCF wave functions give less determinants than FVCAS wave functions for C 2 , N 2 , and O 2 , but the reduction is rather modest with at most about a factor of 2 less determinants for C 2 . For F 2 , the multiple bonding-pattern VBSCF wave function actually generates more singly excited determinants than the FVCAS wave function. This happens because for these small F 2 wave functions, it is more advantageous to restrict the number of single excitations based on the spatial symmetry of the molecular orbitals, as done in the FVCAS calculations, rather than based on a localization criterion of the VBSCF orbitals.

IV. CONCLUSIONS

We have performed VMC and DMC calculations using compact Jastrow-Valence-Bond wave functions based on strictly localized active orbitals. The J-VBSCF wave functions are based on a chemically intuitive description of bonding, including either a single bonding pattern (Fig. 1) or multiple bonding patterns (Fig. 2). The J-BOVB wave functions introduce additional correlation through orbital relaxation.

For N 2 , O 2 , and F 2 , the single bonding-pattern J-VBSCF wave functions yield DMC equilibrium well depths with errors ranging from about 2 to 5 kcal/mol. This is quite a reasonable accuracy considering the compactness of the wave functions. For C 2 , a system with very strong static correlation, a two-bonding-pattern J-VBSCF wave function is necessary to reach a similar accuracy. Optimization of all the wave-function parameters in VMC, and in particular the orbitals, is important to reach this accuracy. As regards computational cost, these VBSCF wave functions contain significantly less Slater determinants than FVCAS expansions, and are therefore more efficient to evaluate. The J-BOVB wave functions are much more costly to evaluate than the J-VBSCF wave functions, and give only moderate improvement of well depths.

The VB ansatz provides a general way to design compact wave functions that capture only the essential physics of the chemical bond, by including only the most chemically relevant structures. These wave functions could be further improved, by including more and more structures belonging to secondary bonding patterns, depending on the accuracy needed.

first-row diatomic molecules this approach is of limited interest since it is easily possible to use FVCAS expansions. However, for molecules requiring larger active spaces, very compact valence bond wave functions can still be designed, and this becomes a great advantage over molecular orbital based multi-determinant approaches. Also, for large systems with regions of interest localized on a few atoms, which is usually the case in chemical reactions, strictly localized orbitals enable both a reduction of the active space and a reduction of the orbital optimization cost, as compared to delocalized molecular orbital approaches.

FIG. 1 .

 1 FIG. 1. Dominant bonding patterns for the C2, N2, O2, and F2 molecules, defining our single bonding-pattern VB wave functions.

FIG. 2 .

 2 FIG.2. Most important resonant bonding patterns for the C2, N2, O2, and F2 molecules, defining our multiple bondingpattern VB wave functions.

TABLE I .

 I Relative contributions of the different types of electron correlation to the equilibrium well depths for C2, N2, O2, and F2. Starting from the RHF well depth, the "left-right" static correlation is defined as the contribution recovered by a single-bonding pattern VBSCF calculation, the other parts of static correlation as the additional contributions recovered by a FVCAS MCSCF calculation, and the dynamic correlation as the remaining contribution necessary to reach the exact well depth.

		Static correlation Dynamic correlation
		left-right other
	C2	77%	18%	5%
	N2	70%	14%	16%
	O2	60%	8%	32%
	F2	60%	8%	32%

TABLE VI .

 VI Numbers of CSFs (NCSF), numbers of unique spin-up + spin-down determinants (N det ) in ground-state wave functions and numbers of unique spin-up + spin-down singly excited determinants (N exc det ) for orbital optimization for the single bonding-pattern (SBP) and multiple bonding-pattern (MBP) VBSCF and BOVB wave functions employed in this work. Values for FVCAS wave functions are also shown for comparison.

			C2		N2		O2		F2
		NCSF N det N exc det NCSF N det N exc det NCSF N det N exc det NCSF N det N exc det
	SBP VBSCF	8	32 1488	10	16 1280	5	10 848	2	4 232
	SBP BOVB					5	32 2144	2	20 1428
	MBP VBSCF	12	64 2992	24	72 5152	10	28 2282	5	16 1572
	FVCAS	165 140 6120	107 112 5864	30	52 2548	8	16 984

ACKNOWLEDGMENTS

BB thanks Philippe C. Hiberty for helpful discussions and Arne Lüchow for comments which helped improve the manuscript. CJU acknowledges support from NSF grant number CHE-1004603.