Christophe Marcenat

Thierry Klein

David Leboeuf

Alexandre Jaoui

Gabriel Seyfarth

Jozef Kačmarčík

Yoshimitsu Kohama

Hervé Cercellier

Hervé Aubin

Kamran Behnia

Benoît Fauqué
email: benoit.fauque@espci.fr

Wide critical fluctuations of the field-induced phase transition in graphite

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

In the immediate vicinity of the critical temperature (Tc) of a phase transition, there are fluctuations of the order parameter, which reside beyond the mean-field approximation. Such critical fluctuations usually occur in a very narrow temperature window in contrast to Gaussian fluctuations. Here, we report on a study of specific heat in graphite subject to high magnetic field when all carriers are confined in the lowest Landau levels. The observation of a BCS-like specific heat jump in both temperature and field sweeps establishes that the phase transition discovered decades ago in graphite is of the second-order. The jump is preceded by a steady field-induced enhancement of the electronic specific heat. A modest (20 percent) reduction in the amplitude of the magnetic field (from 33 T to 27 T) leads to a threefold decrease of Tc and a drastic widening of the specific heat anomaly, which acquires a tail spreading to two times Tc. We argue that the steady departure from the mean-field BCS behavior is the consequence of an exceptionally large Ginzburg number in this dilute metal, which grows steadily as the field lowers. Our fit of the critical fluctuations indicates that they belong to the 3DXY universality class, similar to the case of 4 He superfluid transition.

A phase transition is accompanied by sharp discontinuities of the thermodynamic properties. Quantifying entropy by measuring the specific heat across the transition pins down the order of the transition and informs on the underlying microscopic interaction. Of particular interest is the critical regime of the phase transition, which allows the identification the universality class of the transition and provides direct information on the order parameter [1]. Critical fluctuations are important when their amplitude is comparable with the amplitude of the jump of the specific heat, ∆C, which occurs roughly when the reduced temperature τ = T -Tc Tc is smaller that of others τ G [2], the Ginzburg criterion with:

τ G = α 2 ( k B ∆Cξ 3 m ) 2 (1)
Here, ξ m is the correlation length (averaged in presence of anisotropy) and α = 1 4 √ 2π is a numerical factor. In most cases, τ G 1, the critical fluctuations are located in the extreme vicinity of the transition, and therefore, hardly observable. One notorious exceptions is the 4 He superfluid transition for which the shape of the transition is determined by the critical fluctuation [3]. Near a quantum critical point, thermal fluctuations are replaced by quantum mechanical zero-point fluctuations, which can produce new quantum phases [START_REF] Sachdev | Quantum Phase Transitions[END_REF].

Here, we report on the electronic specific heat (C el ) of graphite, using the state of the art of calorimeters, when all the carriers are confined in the lowest Landau levels (LLL), the so-called quantum limit, which can be easily achieved in this dilute metal. We find that this regime is marked by a steady field-induced enhancement of C el , signaling the enhancement of electron-electron correlations. Deep in this regime, we detect a jump in C el , unambiguously establishing a second-order phase transition induced by the magnetic field. As the magnetic field is reduced, the anomaly shifts to lower temperature and widens. It evolves from a BCS mean-field type transition at 33 T to a cross-over regime below 25 T. At the lowest critical temperature (T c ∼ 1K), fluctuations can be observed up to two times T c . We identify them as critical fluctuations associated with an exceptionally large Ginzburg number due to a low energy condensation. Comparison with a number of other phase transitions shows what distinguishes this phase transition and its criticality. In graphite, the difference in the heat capacity between the normal and ordered phase within a coherence volume rapidly drops with decreasing magnetic arXiv:2011.02435v1 [cond-mat.str-el] 4 Nov 2020 field and becomes exceptionally low. With an electronic density as low as n=4×10 18 cm -3 and light in-plane mass carrier (m a,b =0.05m 0 ), a magnetic field of 7.5T, labelled B QL , oriented along the caxis of graphite is large enough to confine the carriers into the lowest LL (n=0). In the early 80s Tanuma and co-workers [START_REF] Tanuma | Springer Series in Solid State Sciences[END_REF] discovered the onset of an electronic phase transition at B=25T and T=1.3K. Since then, extensive electrical [START_REF] Yaguchi | [END_REF][7][8][9][10][11], thermo-electrical [12] and ultrasound measurements [13] at high magnetic field have established that graphite hosts a succession of, at least, two field induced phases [7,11] arising from electron-hole instabilities. Depending on the nesting vector considered and the strength of the electron-electron interaction, various types of charge [10,14], spin [15] density waves or an excitonic insulating phase [8,9,11,16] have been proposed. To the notable exception of ultrasound measurements [13] these studies have employed transport probes. Due to the low electronic density of graphite, thermodynamic studies are challenging and, nevertheless, crucial. Fig. 1a) shows the field dependence of γ= C el T at T=1.6K up to 35T of graphite. The specific heat setup is described in the supplement [17]. At zero field, γ is as small as 20 ± 3 µJ.mol -1 .K -2 in good agreement with the value expected from the Slonczewski Weiss Mc-Clure (SWM) band model [18,19]. This is several orders of magnitude smaller than in metals due to its low density and lightness of carrier. Sweeping the magnetic field, we found that γ peaks at each LL depopulation. Above B QL , γ increases linearly up to 28T where it presents yet another peak, which evolves with temperature. The temperature and magnetic field dependence of the 28T peak can be tracked by field sweeps (at different temperatures), as shown in the inset of Fig. 1 and by temperature sweeps (at different magnetic fields), as shown on Fig. 2 a). where T p is the temperature of the peak position. At the lowest field the tail of the transition extended up to two times Tc. c) Specific heat anomaly in a BCS transition. The amplitude of the jump at Tc is such that ∆C el γ el Tc =1.43 d) Specific heat anomaly in a BEC transition (in black) and the singularity caused by 3DXY critical fluctuation regime (in green) in the case of the superfluid transition in helium.

Let us begin by commenting the evolution of γ with magnetic field between 10 T and 35 T. Its magnitude in the normal state (labelled γ N, ) can be deduced from temperature sweeps, shown in Fig. 2a), and is represented by open black circles in Fig. 1. γ N increases linearly with magnetic field. The threefold enhancement between 10 T and 35 T is larger than the enhancement of γ reported in Sr 3 Ru 2 O 7 across its quantum critical point [20]. This remarkable enhancement is driven by the change of the density of states (DOS) induced by the magnetic field. When B>B QL the DOS of the LLL is the product of the in-plane degeneracy (which scales linearly with the field) and the one dimensional DOS along the field direction. As long as the Fermi energy (E F ) is far from the bottom of the LLL γ N,el scales linearly with B and with ∂γ ∂B ∝ (

m z E F
) where m z is the mass along the magnetic field and E F (see the supplement [17]). The small E F 40meV and the large c-axis mass m c 10-20m 0 of graphite set the large slope, which we find to be ∂γ N,el ∂B =2.2 ± 0.5µJ.K -2 .mol -1 .T -1 ). Decades ago, starting from the SWM model, Jay-Gerin computed ∂γ ∂B =1.8µJ.K -2 .mol -1 .T -1 [21] in quantitative agreement with our result. At the highest magnetic field (B=33.5T), the specific heat anomaly occurs at T c 3K (see Fig. 2a)). The transition shifts to lower T c with decreasing magnetic field, in contrast to a superconducting transition. By reducing the magnetic field by twenty percent, i.e. by passing from 33.5T to 27.7T, T c decreases by a factor of 3. Fig. 3a) shows how T c derived from the peak in C el T compares with those deduced from anomalies in the in-plane (R xx ) and the out-of-plane resistance (R zz ) measurements. Each symbol represents a T c (B) (or B c (T)) anomaly. The ver-tical axis represents ln(T ) and the horizontal axis B -1 . Thus, the BCS-like expression T c =T exp -B Bc [14] becomes a straight line, which is the behavior of T c (B) according to R xx measurements. In this formula, T and B are phenomenological temperature and field scales and the underlying assumption is that the DOS linearly increases with magnetic field [14]. The onset of ordering has different manifestations in R xx (which jumps by 30% but does not diverge) and in R zz [7,11], which shows an activation behavior and a detectable energy gap (see [17] for complementary measurements). As seen in Fig. 3a), the peak in C el T tracks the onset of the activation energy in R zz and below 3 K, they both occur below the R xx anomaly and the BCS line. Fig. 3b), shows the evolution of the ratio of the activation gap, ∆ c , to the critical temperature with magnetic field, which is close to what is expected in the BCS picture and it weakens with decreasing field. The same figure shows that the normalized specific heat jump ∆C el γTc presents a similar evolution towards weak coupling as the field decreases.

The steady evolution towards weak coupling with decreasing magnetic field is accompanied by a drastic change in the shape of the specific heat anomaly (see Fig. 2a andb). At B=33.5T, the amplitude of the jump is ∆C el γ N Tc =1.1, smaller than what is expected for a meanfield weak-coupling case and has a small tail caused by fluctuations above the critical temperature. With decreasing T c the transition widens, ∆C el γTc decreases and the tail extends to higher temperatures. At B c =27.7T, the tail of the transition extends up to twice T c . The anomalies at two fields differing merely by a factor of 1.2 are compared in Fig. 2. At B=33.5T, the anomaly looks similar to a mean-field BCS transition (sketched in Fig. 2c)), but at B c =27.7T, it acquires a cusp shape and the fluctuating contribution weighs as much as the mean-field jump. This is to be compared and contrasted with the case of the superfluid transition in 4 He (sketched in Fig. 2d)). In the latter case, there is no specific heat jump [22] and the 3DXY critical fluctuations induce the non-analytical behavior of the specific heat across the 'λ'transition [3]. Thus, what distinguishes the case of the field induce state of graphite from a typical superconductor and the archetypal superfluid is the simultaneous presence of a BCS jump and strong fluctuations and the contrasting evolution of these two components of the specific heat anomaly in a narrow field range. . Here T p c is the temperature at which specific heat peaks. ∆C el quantifies the excess in specific heat relative to a mean-field transition. The log-log plot clearly shows that it does not evolve as power law and points therefore to a non-Gaussian origin (for which ∆C Gauss. el = 1 ξ D 0 τ 2-D/2 with D is the dimension and ξ 0 is the BCS coherence length [START_REF] Larkin | Theory of Fluctuations in Superconductors[END_REF]). On the other hand, as shown in the inset of the same figure, a simple empirical law of the form : ∆C el = C 0 exp(-τ τ0 ) captures the evolution of the data with B and τ . The field dependence of the two parameters δC 0 and τ 0 is shown in Fig. 4b). With decreasing magnetic field τ 0 increases and its extrapolation implies τ -1 0 = 0 when B<27 T, where no specific heat anomaly is detected down to 0.6K ( see Fig. 2a). As seen in Fig. 4e), it is possible to collapse all the curves on top of each other by plotting ∆C el /C 0 vs. =T p c (1 + 0.5τ 0 )). The success of this simple scaling procedure has implications for the origin of the broadening of the transition caused by a small variation in the amplitude of the magnetic field. It is unlikely that disorder plays a major role. In this layered material, however, the in-plane and out-of-plane length scales differ by several orders of magnitude. The mean free path within the graphene planes, ⊥ e is very long, in the range of 5µm to 50µm, three orders of magnitude longer than 1

k -1 F,⊥
=7nm, and the magnetic length, B = ( eB )=5nm at 25 T which is a plausible candidate to represent the in-plane coherence length (ξ ⊥ ) beyond the quantum limit. On the other hand, the mean free path along the c-axis, // e is two or three orders of magnitude shorter than ⊥ e . However, it remains longer than the c-axis inter-electron distance, d e-e,//

1 k F,//
1nm and the c-axis coherence length (ξ // ). The latter can be estimated using the BCS coherence length ξ // =ξ 0 = 2 k F,// πm z ∆c . Here ∆ c is the c-axis gap (deduced from the activated behavior of R zz and shown in the inset of Fig. 3b)) and k F,// = π 4a (where a is the interlayer distance). As seen in Fig. ??d) ξ // (33T ) 5nm and ξ // (25T ) 15nm. The latter number may approach the estimated // e below 25T. Therefore the role played by stacking disorder across the planes in causing the broadening cannot be rule out at the lowest magnetic field.

The most plausible source of the observed broadening are critical fluctuations. This can be seen by quantifying the Ginzburg criterion and the expected width of the thermal window for critical fluctuations. Plugging the measured ∆C el and the estimated ξ //,⊥ (ξ m = (ξ 2 // ξ ⊥ )

3 ) in Eq.1, we deduced G . As seen in Fig. 4d), this estimation of τ G closely matches τ 0 deduced from our fits to the data. At the highest field, B=33T τ G 0.05, which is already non negligible. As the field decreases to 27.7T, τ G becomes as large as 0.45, which is exceptionally large in comparison to any other known second-order phase transition. This is seen in Table I which compares our case with a few other systems. The list consists of a conventional superconductor (Sn), two superconductors with higher T c s and shorter coherence lengths (MgB 2 and YBa 2 CuO 7 ), and a Charge Density Wave (CDW) solid (K 0.3 MoO 3 ). One can see that it is the small magnitude of ∆C el , which distinguishes graphite, a system in which 10000 atoms share a single electron and hole. Even at 33 T, the amplitude of the specific heat jump is many orders of magnitude smaller than others. This ∆C el 0.16mJK -1 mol -1 (combined with a molar volume of V m =5.27 cm 3 .mol -1 ) implies that the difference in average specific heat of the normal and ordered phases within a coherence volume becomes comparable to the Boltzmann constant (∆Cξ 3 m 0.3k B ), making this competition critically fragile. With decreasing magnetic field, fluctuations grow and the transition widens because of further decrease in ∆C el . So far we have discussed the temperature dependence of ∆C el through a phenomenological two-parameter exponential fit. As shown in Fig. 4e), on the whole temperature range, it is equally possible to fit ∆C el with a simplified version of the asymptotic form of 3DXY universality class expression with three parameters : A 0 (1 + C 0 |τ 0.5 | + D 0 τ ) (see the supplement [17]). This universality class provides a natural explanation for the saturation of ∆C el at low τ due to its almost vanishing critical exponent (α ≈-0.01 [START_REF] Guillou | [END_REF]). Therefore, we conclude that the observed broadening of the transition is caused by critical fluctuations belonging to the 3DXY model. Let us note that Gaussian fluctuations expected outside the critical window, i.e. for τ > τ G would bring a correction to the mean field behavior. Such a correction is not detected in our data, presumably because of the width and predominance of the critical fluctuations.

Let us also note that below a critical temperature, T p,0 c 1K, and a critical field B p,0 c 25T the entrance in the ordered phase is no more accompanied by a well defined jump in C el but it is replaced by a cross-over. In this case, like in the case of the λ-transition, critical fluctuations dominate the transition. However, in contrast to the λ-transition [3], they are not restricted to the extreme vicinity of the transition.

The present result helps to understand the evolution of the observed anomalies in the Nernst effect [12,29] (a measure of the entropy per carrier [30]) and the ultrasound measurements [13] caused by the transition. At low temperature, one expects that ordering induces a smooth variation in entropy (see [17]) and therefore a rounded drop in the Nernst response [12,29]. As the temperature increases the mean field component of the transition strengthens and the Nernst anomaly becomes a clear kink. Furthermore, the origin of the relatively large jump in the sound velocity becomes clear. It was argued [13] that it can be caused by either a strong anisotropy in the strain dependence of T c or a large jump in the specific heat (for example due to a putative lattice deformation accompanying the transition). According to the present study, ∆C ≈ γT c and the transition is purely electronic. Therefore, we can safely conclude that the jump in the sound velocity is caused by a strong anisotropic strain dependence of T c .

Future studies of the specific heat at higher magnetic field would shed light to the BCS-BEC crossover as one approaches the maximum transition temperature around 47T, where the degeneracy temperature and critical temperature become close to each other [29]. Specific heat studies on other dilutes metals pushed to extreme quantum limit and hosting field-induced state (such as bismuth [31,32], InAs [START_REF] Jaoui | Giant seebeck effect across the field-induced metal-insulator transition of InAs[END_REF], TaAs [START_REF] Ramshaw | [END_REF] or ZrTe 5 [35]) could bring interesting insights.

In summary, we measured the specific heat of graphite in high magnetic field and detected a second-order phase transition jump at high field. The specific heat anomaly drastically evolves in a narrow field window as a consequence of the change in the balance between critical fluctuations and a mean-field jump. The field-induced phase transition in graphite emerges from this study as possessing an exceptionally wide critical window compared to any other electronic phase transition.

FIG. 1 .

 1 FIG.1. Field dependence of γ el = C el T at T=1.6K. The vertical black line indicates the quantum limit regime. When B>BQL=7.5T, holes and electrons are confined to their n=0 Landau levels (with two spin-polarized sub-levels for each). Open circles represents γ in the normal state (labelled γN ), deduced from temperature sweeps shown in Fig.2. The inset shows field sweeps at different temperatures.

FIG. 2 .

 2 FIG. 2. Temperature dependence of γ: a) γ = C el T function of temperature for different magnetic fields. b) C el /T γ N as

FIG. 3

 3 FIG. 3. a) Phase diagram of graphite showing the evolution of critical temperature (Tc) as a function of the inverse of magnetic field, (Tc,B -1 c ), according to different experimental probes. Anomalies in Rxx (black open square points) and Rzz (black open diamonds points) (from [13]) are compared with peaks in C el p c in purple circles). The black line represents Tc=T exp -B Bc and reasonably describes the evolution of Rxx anomalies. b) Field dependence of ∆C el γ N,el Tc and ∆c k B Tc . ∆c is the activation gap deduced from the c-axis resistance measurements (see [17] for the raw datas). The field dependence of ∆c is shown in the inset.

Fig. 4a )

 4a Fig.4a) shows ∆C el =C el -γT as a function of the reduced temperature τ = T -T p c T p c

FIG. 4 . 1 3

 41 FIG. 4. Fluctuations regime of the electronic specific heat: a) ∆C el = C el -γN T as function of τ = T -Tc Tc for τ > 0. The dot lines correspond to a phenomenological fit ∆C el = C0 exp(-τ τ 0 ). Inset shows ln(∆C el ) as function of τ . b) Field dependence of C0 and τ -1 0 . The red dot line is linear fit of τ -1 0 . c) Field dependence of ξ ⊥ = B , ξ // = v F,// π∆c where v F,// is the Fermi velocity along the field, ∆c is the gap deduced from the c-axis resistance (shown on Fig.3b)) and ξm=(ξ 2 ⊥ ξ // ) 1 3 . d) Temperature dependence of τ0 compare with the Ginzburg criterion τG. Inset shows a plot τ -1 0 vs T p c with a linear fit (blue line). e) Normalised plot ∆C el as function of (T -T ren c ) T ren c τ 0 where T ren c =T p c (1 + 0.5τ0). The black dot lines is fit with the 3DXY model (see the text and [17]).

TABLE I .

 I Comparison of τG = 1 32π 2 ( k B Vm is the amplitude of the jump of specific heat at the transition in J.K -1 .m -3 , Vm is the molar volume, ξ0 the BCS coherence length. For anisotropic compounds the average value of ξ0, ξm, is used.

			0 ∆Cξ 3	) 2 for different
	superconductors (SC) (Sn, MgB2 and YBa2CuO7), the CDW
	system (K0.3MoO3) and in graphite. ∆C=∆C el × Sample ∆C(J.K -1 .m -3 ) ξ0( Å)	τG
	Sn (SC)[1]	800	2300	10 -14
	MgB2 (SC)[24]	3360	26	10 -3
	K0.3MoO3 (CDW)[25, 26]	35800	6.2(ξm) 0.01
	YBa2CuO7 (SC) [27]	4993	6.7(ξm) 0.2
	graphite (33.2T)	30.1	47(ξm) 0.05
	graphite (27.7 T)	3.7	68(ξm) 0.45

We thank Zengwei Zhu for useful discussions. We acknowledge the support of the LNCMI-CNRS, member of the European Magnetic Field Laboratory (EMFL). This work was supported by JEIP-Collège de France and by the Agence Nationale de la Recherche (ANR-18-CE92-0020-01; ANR-19-CE30-0014-04). J.K. was supported by the Slovak grants No. APVV-16-0372 and VEGA 2/0058/20.


Steve A Ndengué

Yohann Scribano

David M Benoit

E A Milne

Fabien Gatti

Richard Dawes

Intermolecular rovibrational bound states of H 2 O-H 2 dimer from a MultiConfiguration Time Dependent Hartree approach

Keywords: Quantum Dynamics, Bound states, Time-dependent, Van der Waals, MCTDH

We compute the rovibrational eigenstates of the H 2 O-H 2 Van der Waals complex using the accurate rigid-rotor potential energy surface of Valiron et al. [J. Chem. Phys. 129, 134306 (2008)] with the MultiConfiguration Time Dependent Hartree (MCTDH) method. The J = 0 -2 rovibrational bound states calculations are done with the Block Improved Relaxation procedure of MCTDH and the subsequent assignment of the states is achieved by inspection of the wavefunctions' properties. The results of this work are found to be in close agreement with previous time independent calculations reported for this complex and therefore supports the use of the MCTDH approach for the rovibrational spectroscopic study of such weakly bound complexes.

Introduction

Simulating material properties from first principles at the molecular level is a very challenging goal for all technological materials development or for a more fundamental perspective. This may require an accurate ab initio potential energy surface (PES) which should be accurate enough to be used for a large range of applications starting from gas-phase to the study of the bulk-phase properties. For many years, the field of van der Waals molecules, or weakly bound complexes (implying non-covalent bonding between atoms), has emerged and their molecular spectroscopic properties (like tunneling splittings in water clusters) have been used to test or/and develop accurate manybody PESs. Usually those molecular cluster PESs are expanded as a sum of low-order many body terms, as has been reported for water clusters [START_REF] Bukowski | Predictions of the properties of water from first principles[END_REF][START_REF] Huang | New ab initio potential energy surface and the vibration-rotation-tunneling levels of (H 2 O) 2 and (D 2 O) 2[END_REF][START_REF] Shank | Accurate ab initio and hybrid potential energy surfaces, intramolecular vibrational energies, and classical ir spectrum of the water dimer[END_REF][START_REF] Wang | Flexible, ab initio potential, and dipole moment surfaces for water. I. Tests and applications for clusters up to the 22-mer[END_REF][START_REF] Medders | On the representation of many-body interactions in water[END_REF][START_REF] Cisneros | Modeling Molecular Interactions in Water: From Pairwise to Many-Body Potential Energy Functions[END_REF][START_REF] Paesani | Getting the right answers for the right reasons: toward redictive molecular simulations of water with many-body potentials[END_REF]. It was shown that including up to three-body terms is enough to correctly reproduce the infrared and Raman spectro-copy of liquid water [START_REF] Medders | Infrared and Raman Spectroscopy of Liquid Water through "First-Principles" Many-Body Molecular Dynamics[END_REF][START_REF] Medders | Many-body convergence of the electrostatic properties of water[END_REF]. The accuracy of the ab initio two-body term, the dimer potential, may be evaluated by a direct comparison between observed experimental spectrum and theoretical spectrum obtained by solution of the Schrödinger equation as was done recently in the study of the water dimer [START_REF] Wang | Using monomer vibrational wavefunctions to compute numerically exact (12D) rovibrational levels of water dimer[END_REF][START_REF] Babin | Development of a first principles water potenital with flexible monomers: Dimer potential energy surface[END_REF].

The H 2 O-H 2 van der Waals complex is an example of a molecular dimer with relevance in domains beyond chemical physics. It is indeed an important system within the astrophysical context since H 2 O is the third most abundant polyatomic system after CO and H 2 . Thus the complex involves two of the most abundant molecules in the Interstellar Medium (ISM), which are also important on Earth. The knowledge of an accurate PES is thus required in order to determine the efficiency of the molecular (H 2 O) cooling through its inelastic collisions with molecular hydrogen [START_REF] Philips | Collision dynamics for an asymmetric top rotor and a linear rotor: Coupled channel formalism and application to H 2 O-H 2[END_REF][START_REF] Dubernet | Collisional excitation rates of H 2 O with H 2 . I. Pure rotational excitation rates with para-H 2 at very low temperature[END_REF][START_REF] Faure | Quasi-classical rate coefficient calculations for the rotational (de)excitation of H 2 O by H 2[END_REF][START_REF] Dubernet | Rotational excitation of ortho-H 2 O by para-H 2 (j 2 =0, 2, 4, 6, 8) at high temperature[END_REF][START_REF] Scribano | Communication: Rotational excitation of interstellar heavy water by hydrogen molecules[END_REF][START_REF] Scribano | Rotational excitation of H 2 O by para-H 2 from an adiabatically reduced dimensional potential[END_REF]. This van der Waals complex is also of interest for gas storage materials such as hydro-gen clathrate hydrates [START_REF] Mao | Clathrate hydrates under pressure[END_REF][START_REF] Schüth | Hydrogen and hydrates[END_REF][START_REF] Homayoon | Fulldimensional, high-level ab initio potential energy surfaces for H 2 (H 2 O) and H 2 (H 2 O) 2 with application to hydrogen clathrate hydrates[END_REF] due to their potential as efficient and environmentally friendly materials for hydrogen storage. Indeed, most of the theoretical studies of one or more H 2 molecules encapsulated in clathrate hydrate make the assumption that the interaction between the guest molecule and the host ice-water structure can be described as a sum of pairwise water-hydrogen dimer interactions [START_REF] Xu | Hydrogen molecule in the small dodecahedral cage of a clathrate hydrate: Quantum five-dimensional calculations of the coupled translation-rotation eigenstates[END_REF][START_REF] Powers | The effect of the condensed-phase environment on the vibrational frequency shift of a hydrogen molecule inside clathrate hydrates[END_REF].

The quantum dynamics of H 2 O-H 2 is complicated by the inherently large amplitude motions which characterize weakly-bounded complexes with light-atoms and their associated large zero-point-energies (ZPE). All intermolecular motions are delocalized over the various minima of the PES. The previous theoretical spectroscopy studies have used a time-independent approach to solve the Schrödinger nuclear equation and compute the rovibrational bound states within the rigid rotor model. Those previous spectroscopic calculations have been performed by Wang and Carrington [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF] using a combination of Discrete Variable Representation (DVR) and Finite Basis Representation (FBR) with Lanczos diagonalization as well as by Avoird [START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF] who used a Close Coupling (CC) methodology combined with a Davidson algorithm to determine the rovibrational states. Those two time independent calculations have also used the accurate five dimensional ab initio potential of Valiron et al. [START_REF] Valiron | R12-calibrated H 2 O-H 2 interaction: Full dimensional and vibrationally averaged potential energy surfaces[END_REF] which was obtained from the original nine-dimensional potential energy surface by averaging over the intramolecular vibrational ground state probability densities of water and hydrogen.

In this work, we describe the numerical solution of the rovibrational Schrödinger equation for this complex using an efficient time-dependent method: the MCTDH approach. This approach was recently used by Ndengué et al. [START_REF] Ndengué | Atom-Triatom Rigid Rotor Inelastic Scattering with the MultiConfiguration Time Dependent Hartree approach[END_REF] in an inelastic scattering study of the Ar-H 2 O atom-triatom complex and has demonstrated remarkable accuracy in comparison to CC calculations on this same system. This work on the water-hydrogen complex is an immediate continuation of the previous work and represents the first application of the MCTDH methodology for the rovibrational spectroscopy of a diatomic-triatomic system.

This article is organized as follows: first we will discuss the computational procedure, introducing the MCTDH method, the Hamiltonian (kinetic operator and the potential energy surfaces) and how the rovibrational states calculations are done. Next we will present and discuss our results and conclude with perspectives.

Computational procedure

Rigid Body diffusion monte carlo calculations

In order to benchmark the ground state calculation, we have performed the ground state calculation of the para-para H 2 O-H 2 complex using Diffusion Monte Carlo (DMC). Indeed, this is a practical method to solve accurately the many-body vibrational Schrödinger equation and is particularly suited to the determination of the vibrational ground state. DMC replaces the resolution of the time-independant Schrödinger equation by the simulation of a diffusion process in imaginary time τ = it on a given multidimensional PES. In particular, DMC avoids the convergence problems that can be encountered in discretevariable representation or in spectral approaches due to the finite basis set representation, and thus is often used as a reference method to provide benchmark values. We use the rigid-body version of DMC (RB-DMC), initially proposed by Buch [START_REF] Buch | Treatment of rigid bodies by diffusion Monte Carlo: Application to the para-H 2 • • • H 2 O and ortho-H 2 • • • H 2 O clusters[END_REF], and used by other authors [START_REF] Benoit | Quantum simulation of phenol-water clusters[END_REF][START_REF] Viel | Importance sampling in rigid body diffusion monte carlo[END_REF] to compute the vibrational ground state of weakly bounded clusters. For these molecular systems, decoupling the fast intramolecular modes from the slow intermolecular modes is a very good approximation and we can consider each rotor as a rigid body since this approximation was already validated in the case of water dimers [START_REF] Kelly | Water dimer vibration-rotation tunnelling levels from vibrationally averaged monomer wavefunctions[END_REF]. Our RB-DMC simulations are performed with the Xdmc code developed by Benoit [START_REF] Benoit | Xdmc, Quantum Diffusion Monte Carlo code[END_REF] (see also Ref. [START_REF] Benoit | Quaternion formulation of diffusion quantum monte carlo for the rotation of rigid molecules in clusters[END_REF] for implementation details). In this study we use 500 replicas, a stabilization period of 11 × 100 cycles with ∆τ = 30 a.u. and an averaging phase of 2000 × 100 cycles with ∆τ = 5 a.u. We get a RB-DMC vibrational ground state of -33.58 ± 0.16 cm -1 in very good agreement with the previous values of Wang and Avoird. The water-hydrogen system is indeed strongly anharmonic and has a very weakly bound vibrational ground state.

The MCTDH method

The MultiConfiguration Time Dependent Hartree [START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Beck | The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets[END_REF] method is a time-dependent method in which each degree of freedom is associated with a small number of orbitals or single particle functions (SPFs) which, through their time dependence, allow efficient description of the molecular dynamical processes. The total MCTDH wave function is expanded in Hartree products, that is, products of single-particle functions

Ψ(Q 1 , . . . , Q f , t) = n1 j1=1 • • • n f j f =1 A j1•••j f (t) f κ=1 φ (κ) jκ (Q κ , t) = J A J Φ J , (1)
where f is the number of degree of freedom of the system, Q 1 , . . . , Q f are the nuclear coordinates, A j1•••j f denotes the MCTDH expansion coefficients, and φ

(κ)
jκ (Q κ , t) are the n κ SPFs associated with each degree of freedom κ. The subsequent equation of motion for the coefficients and single particle functions are derived after replacing the wave function ansatz into the time dependent Schrödinger equation. To solve the equations of motion, the κ SPF are represented on a primitive basis or discrete variable representation (DVR) grid of N κ points, where in general the n κ of Equation ( 1) is such that n κ < N κ . Thus, the MCTDH method propagates a wave function on a small time-dependent, variationally optimized basis set of single-particle functions, which in turn are defined on a time-independent primitive basis set.

ϕ (κ) jκ (Q κ , t) = Nκ iκ=1 c (κ) iκjκ (t)χ (κ) iκ (Q κ ) (2)

The Kinetic Energy Operator

For the MCTDH algorithm to be efficient, the Hamiltonian operator must be written as a sum of products (SOP) of single-particle operators. The Kinetic Energy Operator (KEO) has always the required form when using polyspherical coordinates, such as the Jacobi coordinates used in this work where we followed the approach of Ref.35 the subsystem KEO derivation presented by Gatti and Iung [START_REF] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF] which was used recently to describe an asymmetric rotor -atom collision [START_REF] Ndengué | Atom-Triatom Rigid Rotor Inelastic Scattering with the MultiConfiguration Time Dependent Hartree approach[END_REF]. As was the case with previous MCTDH calculations involving fragments, we do not work in the Body-Fixed frame but in the E 2 frame [START_REF] Gatti | Exact and constrained kinetic energy operators for polyatomic molecules: The polyspherical approach[END_REF] which is obtained by rotation of the two first Euler angles of the SF frame. A choice of the vectors for this calculation is displayed in Figure 1. The KEO in the E 2 frame is generally expressed as

2 T = - 1 µ R ∂ 2 ∂R 2 + 2 TA + 2 TB + 1 µ R R 2 J † J + ( L A + L B ) 2 -2( L A + L B ) J E2 , (3)
where µ R is the reduced mass of the H 2 O-H 2 cluster, the A and B subscripts refer respectively to the H 2 O and H 2 fragments. J, L A and L B are respectively the total angular momentum of the system and the angular momenta of the fragments A (H 2 O) and B (H 2 ). The rigid rotor Hamiltonian of the H 2 O molecule is expressed as [START_REF] Leforestier | Grid method for the Wigner functions. Application to the van der Waals system Ar-H 2 O[END_REF]

TA = A 2 L 2 A,+ + L 2 A,-+ L A,+ L A,-+ L A,-L A,+ - C 2 L 2 A,+ + L 2 A,--L A,+ L A,--L A,-L A,+ +BL 2 z BF A , (4)
where the rotational constants have values A = 27.8572 cm -1 , B = 14.5145 cm -1 and C = 9.2799 cm -1 . The rigid rotor kinetic energy of the H 2 fragment is written simply as TB = B H2 L 2 B where the rotational constant B H2 = 59.2434 cm -1 . The final form of the KEO translated and implemented in the MCTDH code is then

2 T = - 1 µ R ∂ 2 ∂R 2 + 2 TA + 2 TB + 1 µ R R 2 J(J + 1) + L 2 A + L 2 B -2L 2 A,z -2L 2 B,z + 1 µ R R 2 (L A,+ L B,-+ L A,-L B,+ -2L A,z L B,z ) + 1 µ R R 2 (C + (J, K)(L A,+ + L B,+ )) + 1 µ R R 2 (C -(J, K)(L A,-+ L B,-)) , (5)
with

L A(B),± = L A(B),x ± iL A(B),y , (6)
and C ± (J, K) = (J(J + 1) -K(K ± 1). ( 7)

The Potential Energy Surface

The potential energy operator just as the KEO needs to be expressed in the product form. For systems of low dimensionality which are not already in the product form, there exists an efficient fitting procedure (Potfit) [START_REF] Jäckle | Product representation of potential energy surfaces[END_REF][START_REF] Jäckle | Product representation of potential energy surfaces. II[END_REF] implemented in the MCTDH package [START_REF] Worth | The MCTDH Package, Version 8.2[END_REF] to obtain the desired representation. For high dimensionality potentials MultiGrid [START_REF] Peláez | The multigrid POTFIT (MGPF) method: Grid representations of potentials for quantum dynamics of large systems[END_REF], Multi-Layer [START_REF] Otto | Multi-layer Potfit: An accurate potential representation for efficient high-dimensional quantum dynamics[END_REF] and Monte Carlo [START_REF] Schröder | Transforming high-dimensional potential energy surfaces into sum-of-products form using Monte Carlo methods[END_REF] versions exist to transform a general potential to the product form. However those implementations often lead to a large number of potential terms which will slow down the computations, for example in scattering calculations where several long propagations are necessary to obtain the cross-sections. In this work, we are able to go around this difficulty since the potential [START_REF] Valiron | R12-calibrated H 2 O-H 2 interaction: Full dimensional and vibrationally averaged potential energy surfaces[END_REF] is already provided in a multipolar form and we only had to make a transformation to express it in the coordinates representation used for the MCTDH calculations.

The surface from Ref.25 is expressed as

V (R, θ, ϕ, θ , ϕ ) = i v i (R) ti (θ, ϕ, θ , ϕ ) (8)
where i = {l A , m A , l B , l} and tl A ,m A ,l B ,l is expressed as

ti (θ, θ , ϕ, ϕ ) = 1 2π α l A ,m A l A r1=-l A l B r2=-l B β i,r1 l A l B l r 1 r 2 r × P r2 l B (cosθ )P r l (cosθ)cos(r 2 ϕ + rϕ) (9)
with

α l A ,m A = √ 2l A + 1 2 + 2δ m A ,0 , (10)
and

β i,r1 = 1 2π δ m A ,r1 + (-1) l A +m A +l B +l δ -m A ,r1 . ( 11
)
The maximum order of the terms involved in the expansion are l A = 11, l B = 6, resulting in 149 v i (R) terms in equation ( 8) with m A (≥ 0), l A and l A + l B + l always even because of symmetry considerations. The BF of the original PES is defined such that its origin is the center of mass of the H 2 O molecule, the z -axis is its C 2 axis with the positive z in the direction of the O atom and the xz plane is the molecular plane. The angles (θ, ϕ) and (θ , ϕ ) represent respectively the collisional direction and the H 2 fragment orientation in the Body-Fixed (BF) frame. In the MCTDH implementation, the dynamics in the E 2 frame is described by 6 coordinates, R the fragments' separation and 5 angles: α A , β A , γ A , θ B and φ B . The origin of the E 2 frame is the center of mass of H 2 O as was the case with the BF frame. The z -axis is the axis of R, the connecting vector between the centers of mass of H 2 O and H 2 . The 3 Euler angles α A , β A , γ A give the orientation of of H 2 O in our frame while the other 2 spherical angles θ B and φ B give the orientation of the H 2 fragment. We point out that the BF Frame and the E 2 Frame mentioned here correspond respectively to the Molecule Fixed (MF) Frame and the Dimer Fixed (DF) Frame described in the work of Wang and Carrington [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF]. A cut of the PES expressed in the set of coordinates used by Valiron et al is shown in Figure 2.

As Van der Avoird and Nesbitt [START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF] noticed, there are two ways to transform the PES from the MF Frame set of coordinates to DF Frame set of coordinates. In the first approach the coordinates are related according to θ = β A , ϕ = π -γ A and θ , ϕ can be expressed in terms of α A , β A , γ A , θ B , φ B with the use of the inverse Euler rotation matrix; then the PES in the appropriate frame can be numerically generated. This procedure, used by Wang and Carrington, [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF] is not convenient for this work as it would require to build a new SOP expansion of the PES after the surface is numerically transformed. The second approach, [START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF] which is the one we selected, relates the coefficients of the multipolar expansion in both frames and thus allows one to generate the PES exactly in the E 2 Frame directly in the SOP form. The correctness of the This second approach saves a significant amount of time and allows a lot of flexibility in the calculations.

Rovibrational states calculations with MCTDH

The rovibrational bound states of the H 2 O-H 2 complex are obtained with the block improved relaxation method [START_REF] Doriol | Computation of vibrational energy levels and eigenstates of fluoroform using the multiconfiguration time-dependent hartree method[END_REF][START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF], implemented in the MCTDH package [START_REF] Worth | The MCTDH Package, Version 8.2[END_REF]. The block improved relaxation is derived from the improved relaxation method [START_REF] Meyer | Calculation and selective population of vibrational levels with the Multiconfiguration Time-Dependent Hartree (MCTDH) algorithm[END_REF], an MCSCF approach where the SPFs are optimized by relaxation [START_REF] Kosloff | A direct relaxation method for calculating eigenfunctions and eigenvalues of the Schrödinger equation on a grid[END_REF] (propagation in imaginary time) and the coefficients vector (A-vector) is determined by diagonalization of the Hamiltonian matrix evaluated in the set of SPFs. The working equations of the improved relaxation and block improved relaxation have already been presented elsewhere [START_REF] Meyer | Quantum molecular dynamics: propagating wavepackets and density operators using the multiconfiguration time-dependent Hartree method[END_REF][START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF] and thus aren't given here.

The primitive basis, its range and the number of the SPFs used for the calculations of the rovibrational states are summarized in Table 1. We see that not more than 48,000 Single Particle Functions are used to run the calculations. As the calculation method is the Block Improved Relaxation, we have to perform several of those to obtain the complete set of expected rovibrational levels. Each calculation for a block of 6 functions with 48,000 SPFs and using 32 processors on a Linux based cluster takes approximately 48 hours. While the time for each calculation is not very long, the blocks are typically done sequentially, waiting for each block to complete (or be near completed) before starting the next block of higher states: this procedure adds to the overall duration of the calculation which therefore could even take longer than traditional time-independent calculations, depending on the number of desired states. However, quite often in this type of application, one is interested only in a limited number of low lying states making this procedure advantageous. A primitive basis composed of Fast Fourier Transform (FFT) functions for the intermolecular distance R was coupled with a Wigner-DVR [START_REF] Leforestier | Grid method for the Wigner functions. Application to the van der Waals system Ar-H 2 O[END_REF] basis where α A , γ A are replaced by their momentum representation k α , k γ and a two dimensional extended Legendre-K DVR [START_REF] Sukiasyan | Reaction Cross Sections for the H+D 2 (ν 0 =1) → HD+D and D+H 2 (ν 0 =1) → DH+H Systems. A Multiconfiguration Time-Dependent Hartree (MCTDH) Wave Packet Propagation Study[END_REF] (replacing φ B by k φ ) to describe respectively the orientation of the H 2 O and the H 2 fragments in the E 2 frame. The Wigner and Legendre-K DVRs switch between their respective FBR (Wigner and Spherical harmonics) and DVR to calculate the action of the derivatives of operators (and the singularities) and to evaluate the potential respectively. Leforestier [START_REF] Leforestier | Grid method for the Wigner functions. Application to the van der Waals system Ar-H 2 O[END_REF] defined the unitary transformation that goes from the 3D DVR to the 3D FBR: this transformaton is what is implemented in the MCTDH code for the Wigner-DVR, and a 2D version for the Legendre-K DVR. For testing purposes, we also ran (J = 0) calculations with the angular primitive basis composed of the "usual" Wigner-DVR basis and a two dimensional Legendre DVR. The (Wigner,K,K) and (KLeg,K) DVRs are then replaced by their counterparts the (Wigner,Exp,Exp) and (PLeg,Exp) DVRs, where KLeg and PLeg are respectively the extended Legendre DVR and two dimensional Legendre DVR. Results of the two sets of calculations were confirmed to be identical, with the second set running faster than the first one. However that approach is not yet implemented in the MCTDH package for J > 0: thus, in the rest of this work we used only the first method for our calculations.

After transforming the PES from Equation [START_REF] Medders | Infrared and Raman Spectroscopy of Liquid Water through "First-Principles" Many-Body Molecular Dynamics[END_REF] to the E 2 frame we can write

V (R, β A , γ A , α A , θ B , φ B ) = r β ,rγ rα,r θ Ṽr β ,rγ rα,r θ (R)frβ,rγ rα,r θ (ω A , ω B ) (12)
where

frβ,rγ rα,r θ (ω A , ω B ) = D (r β ) rαrγ (α A , β A , γ A ) C r θ ,-rα (θ B , φ B ) (13)
with D and C being respectively the Wigner D-matrix and the Racah normalized spherical harmonics. The action of the potential on the wavefunction can then be obtained using k α , k γ and k φ to express the angles α A , γ A and φ B , where in the following we drop the indices A and B for simplicity.

V Ψ(R, β, k γ , k α , θ, k φ ) = r β ,rγ rα,r θ Vr β ,rγ rα,r θ (R, β, k γ , k α , θ, k φ ) ×Ψ(R, β, k γ -r γ , k α -r α , θ, k φ + r α ). ( 14
)
The primitive basis set was selected from previous work [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF] done on this system for which a convergence of the results to 0.002 cm -1 or less for the low lying states was re- ported. Also, the number of SPFs was increased in the calculations from a relatively small number for the lower levels to significantly higher values for higher excited states. This however presents a challenge since currently symmetry is not implemented for the Wigner primitive basis in the MCTDH package. The calculations for para/ortho H 2 O and H 2 were done without differentiating these states and the results we obtained were general and did not separate bound from resonance states for the various symmetries. For this work we will only address bound states of this system as reliable data are available for comparison. The calculations of resonances will be described and discussed in future work where we will take advantage of a recently presented procedure to obtain resonances with MCTDH [START_REF] Ndengue | Resonances of HCO Computed Using an Approach Based on the Multiconfiguration Time-Dependent Hartree Method[END_REF].

Results and discussion

Rovibrational states and symmetry

The rovibrational states of the H 2 O-H 2 cluster for J = 0, 1, 2 were computed with the method described above. For this work the assignment was straightforward as two sets of data [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF][START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF] (obtained with the same PES) which agreed very well with our calculations are already available in the literature and thus served as reference for the assignment. Comparisons of our calculations with previous results are displayed in Tables 3,4, 5 and 6. It is worth pointing out that the difference between the two sets of previous results probably arise from the difference in the rotational constants values used by each group of authors; in our work we used the values used by Wang and Carrington [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF]. Note also that their computed pH 2 O-pH 2 ground state also agrees well with our RB-DMC results.

As described previously, the assignment of states resulting from MCTDH calculations (at least in its current implementation) while not trivial, is about as difficult as it is with other methods for rigid rotor Van der Waals systems. Though the calculations that were performed here did not include the symmetries of the dimer, it is possible [START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF] or Wang and Carrington [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF], the primitive basis was constrained such that K, the projection of the total angular momentum writes K = m A + m B . Then, for example, K = 0 requires m A = -m B and all the states that will result from their rovibrational calculations will present those characteristics. While this type of constraint is difficult to implement within the MCTDH program, it is still possible to obtain from the results of our calculations the rovibrational states that fullfill those conditions for a specific J and K and discard the other dummy states (with J = 0 and K = m A + m B = ±1 for example) which may also appear in the calculations. We show in Table 2 the convergence of the first few states for J = 0, an output from the MCTDH calculations. We use for these calculations a primitive basis which is similar in its range and size to the results reported by Wang and Carrington [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF], but where we varied the number of SPFs. We can see from the Table the convergence of the calculation, but also the dummy (or spurious) states in italic, as mentioned earlier.

The results of the rovibrational calculations from the Heidelberg MCTDH program provide relevant information that can be used to assign the various states. First, it is possible to remove spurious or non-physical rovibrational states (as we mentioned before) with m A + m B > J from the calculations by individually analysing each state, doing for example a short time single-state calculation. The value of m A + m B allows also one to determine the character of the rovibrational state (Σ, Π, ∆, ...). It is worth pointing out that for our calculations we didn't notice any spurious states but for another system (H 2 O-HCN) currently being investigated this check is relevant. Secondly, in a block calculation where no use of symmetry Table 3: Rovibrational energy levels (cm -1 ) of pH 2 O-pH 2 for J = 0, 1, 2. The parity (P ) e and o stand for even and odd respectively. The character and dominant j kakc are expressed as Λ(j kakc ) where Λ can be Σ, Π, ∆ for K = 0, 1, 2 (see text for details). Comparisons are made with other calculations reported for the same PES: Wang & Carrington [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF] and Van der Avoird & Nesbitt [START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF]. Energies are in cm -1 . P Λ(j kakc ) Ref. [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF] Ref. [START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF] was made such as in our calculations, the character ortho or para of the rovibrational states could be determined from the expectation value of the grid population of the θ DOF. In other words the average value j of the rovibrational state where j is the angular momentum quantum number of the H 2 fragment. As it is well known for H 2 , para and ortho states have respectively an even and odd j progression and therefore the value of j provides the necessary information to differentiate the states. Alternatively, using the Extended Legendre DVR filtered with even or odd functions could also allow one to differentiate pH 2 from oH 2 states. Third, the j of the dominant j kakc use in the assignment can be determined from the expectation value of the grid population of the β DOF. Finally, as an alternative to all the previous methods, one could use projection methods like the one described by Wang and Carrington [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF] to assign the states.

Conclusion

We report calculations of rovibrational energies on the H 2 O-H 2 system in the rigid rotor approximation, obtained with the MCTDH algorithm on the Valiron PES [START_REF] Valiron | R12-calibrated H 2 O-H 2 interaction: Full dimensional and vibrationally averaged potential energy surfaces[END_REF]. Our results are found to be in very close agreement with the values previously reported. The calculations performed here with the MCTDH program are relatively fast, computationally cheaper and relatively easier to implement than P Λ(j kakc ) Ref. [START_REF] Wang | Theoretical study of the rovibrational spectrum of H 2 O-H 2[END_REF] Ref. [START_REF] Van Der Avoird | Rovibrational states of the H 2 O-H 2 complex: An ab initio calculation[END_REF] standard calculations and therefore appear as a bright prospect for the investigation of similar or other even more complex systems. Additionally, the extension to full dimensionality (9D) calculations on this system can be easily implemented and performed at a lower cost than other traditional methods, even though it may become more difficult to implement the symmetry in higher dimensionality. This work opens the path to a systematic study of various similar systems [START_REF] Kalugina | Interaction of H 2 O with CO: potential energy surface, bound states and scattering calculations[END_REF] and is an important step towards our study of H 2 O+H 2 collisions which will be presented in a future paper.

Figure 1 :

 1 Figure 1: Definition of the Jacobi vectors for the H 2 O-H 2 dimer. G j is the center of mass of the jth monomer. R 1 joins the center of mass of G 1 and G 2 and defines the z -axis. In the current calculations the R (j) i vectors have fixed lengths.

Figure 2 :

 2 Figure 2: 2D cut of the PES in the set of coordinates used by Valiron et al[25]. The figure shows the global and secondary minima of the PES.

  then be verified by comparison of the PES values with the ones generated from the first approach.

Table 1 :

 1 Parameters of the primitive basis used for the rovibrational calculations of H 2 O-H 2 . FFT stands for the Fast Fourier Transform. Wigner stands for the Wigner DVR. KLeg is the extended Legendre DVR. The units for distance and angle are bohrs and radians respectively.

	Coords. Primitive Number of	Range	SPF
		Basis	Points		basis
	R	FFT	96	4.0-26.0 10-20
	β	Wigner	9	0-π	20-80
	γ	K	17	-8,8
	α	K	17	-8,8
	θ	KLeg	9	0-π	20-30
	φ	K	17	-8,8

Table 2 :

 2 Convergence of the rovibrational energies (cm -1 ) of H 2 O-H 2 for J = 0. The results are displayed for various SPF basis and reported as n 1 /n 2 /n 3 where the n i are the number of SPF functions for R, (β,γ,α), (θ,φ). DVR and KLeg-DVR and thus run 4 separate type of calculations to treat the specific pH 2 O-pH 2 , pH 2 O-oH 2 , oH 2 O-pH 2 and oH 2 O-oH 2 cases presented in previous work. Thus in our computations, amongst the results, we expected to obtain the real bound states for some cases along with resonances from other cases as we do successive block improved relaxation calculations. In previous work by van der Avoird and Nesbitt

	10/40/20 20/60/25 20/70/30 10/80/30 20/90/40
	-33.551	-33.562	-33.562	-33.563	-33.563
	-12.826	-12.845	-12.845	-12.845	-12.846
	-7.412	-7.418	-7.421	-7.421	-7.421
	-7.412	-7.418	-7.421	-7.421	-7.421
	-2.057	-2.064	-2.064	-2.064	-2.064
	1.418	1.411	1.411	1.411	1.411
	within the MCTDH code to implement and use symmetry
	of the Wigner-

Table 4 :

 4 Same as Table 3 for oH 2 O-pH 2 .

					This Work
			J = 0
	e	Σ(0 00 )	-33.576	-33.567	-33.563
	e	Σ(0 00 )	-2.064	-2.062	-2.064
			J = 1
	o	Σ(0 00 )	-32.194	-32.184	-32.181
	o	Σ(0 00 )	-1.400	-1.398	-1.398
			J = 2
	e	Σ(0 00 )	-29.443	-29.433	-29.430
	e	Σ(0 00 )	-0.163	-0.163	-0.156
	P Λ(j kakc ) Ref. [23] Ref. [24] This Work
			J = 0
	e	Σ(1 01 )	-12.850	-12.834	-12.846
	o	Σ(1 10 )	15.618	15.657	15.638
	e	Σ(1 01 )	21.939	21.946	21.944
			J = 1
	o	Σ(1 01 )	-12.071	-12.055	-12.065
	e	Π(1 01 )	-6.077	-6.062	-6.058
	o	Π(1 01 )	-5.452	-5.437	-5.435
	e	Π(1 10 )	6.951	6.992	6.963
	o	Π(1 10 )	7.314	7.355	7.325
	e	Σ(1 10 )	17.300	17.339	17.320
	o	Σ(1 01 )	21.885	21.891	21.897
	e	Π(1 01 )	22.961	22.967	22.963
	o	Π(1 01 )	23.614	23.619	23.613
			J = 2
	e	Σ(1 01 )	-10.261	-10.245	-10.253
	o	Π(1 01 )	-3.367	-3.352	-3.349
	e	Π(1 01 )	-1.759	-1.744	-1.743
	o	Π(1 10 )	9.073	9.114	9.085
	e	Π(1 10 )	10.094	10.135	10.105
	o	Σ(1 10 )	20.579	20.618	20.598
	e	(1 01 )	22.524	22.529	22.525

Table 5 :

 5 Same as Table 3 for pH 2 O-oH 2 .

	P Λ(j kakc ) Ref. [23] Ref. [24] This Work
			J = 0
	e	Σ(0 00 )	64.879	65.083	64.892
	e	Σ(0 00 )	109.911	110.115	109.913
	e	Σ(1 11 )	115.120	115.335	115.075
	o	Σ(1 11 )	117.733	117.956	117.958
			J = 1
	o	Σ(0 00 )	66.294	66.498	66.307
	o	Π(0 00 )	80.176	80.378	80.233
	e	Π(0 00 )	80.217	80.419	80.274
	o	Σ(0 00 )	110.897	111.100	110.898
	e	Π(1 11 )	111.370	111.585	111.378
	o	Π(1 11 )	111.442	111.656	111.439
	o	Σ(1 11 )	116.408	116.623	116.448
	e	Π(0 00 )	117.361	117.561	117.393
	o	Π(0 00 )	117.609	117.807	117.616
			J = 2
	e	Σ(0 00 )	69.114	69.318	69.127
	e	Π(0 00 )	82.910	83.112	82.966
	o	Π(0 00 )	83.031	83.232	83.086
	e	Σ(0 00 )	112.766	112.966	112.770
	o	Π(1 11 )	113.602	113.816	113.610
	e	Π(1 11 )	113.761	113.974	113.756
	o	Π(0 00 )	118.451	118.654	118.566

Table 6 :

 6 Same as Table 3 for oH 2 O-oH 2 .

Acknowledgment

RD is supported by the US National Science Foundation (No. CHE-1566246). S.A.N thanks Hans-Dieter Meyer for various discussions related to the MCTDH implementation.

Highlights

 Rovibrational states of the H 2 O-H 2 complex were computed using the MCTDH method.  Precise agreement between these calculations and previous time-independent calculations was obtained.

 This confirmation of the accuracy and efficiency of the MCTDH method in this context opens up many possible applications.