Sotirios Fragkos

Roberto Sant

Carlos Alvarez

Alexei Bosak

Polychronis Tsipas

Dimitra Tsoutsou

Hanako Okuno

Gilles Renaud

Athanasios Dimoulas

Room temperature commensurate charge density wave in epitaxial strained TiTe2 multilayer films

The group IVB 2D transition metal dichalcogenides (TMD) are considered to be stable in the high symmetry trigonal octahedral (1T) structure due to the lack of unpaired delectrons on the metal site. It is found here that multilayer epitaxial TiTe2 is an exception adopting a commensurate 2×2×2 charge density wave (CDW) structure at room temperature with an ABA type of stacking as evidenced by direct lattice imaging and reciprocal space mapping. The CDW is stabilized by highly anisotropic strain imposed by the substrate with an out-off-plane compression which reduces the van der Waals gap increasing the interlayer coupling. A weaker 2×2 CDW is also confirmed at RT for epitaxial monolayer TiTe2. The addition of epitaxial strained TiTe2 to the family of CDW materials is expected to deepen our understanding on the CDW formation mechanisms in other TMD materials, still under debate and enable real world applications that take advantage of a CDW ground state at room temperature.

The group IVB 2D TMDs 1,2,3 with the chemical formula MX2 (M=Ti, Hf, Zr; X=S, Se, Te) are generally stable adopting a high symmetry trigonal octahedral (1T) structure 4,5 since they do not possess unpaired d electrons 6 at the metal site that could otherwise be a source of instability. A notable exception is 1T-TiSe2 which, below 205 K makes the transition to a lower symmetry commensurate CDW ground state associated with a 2×2×2 periodic lattice distortion (PLD). Despite a large number of studies 2,3 over the years since the first discovery 7 and a couple of comprehensive reviews 8,9 the actual mechanism for PLD/CDW formation is still under debate.

The most recent experimental [START_REF] Chen | Charge density wave transition in single-layer titanium diselenide[END_REF][START_REF] Chen | Dimensional Effects on the Charge Density Waves in Ultrathin Films of TiSe2[END_REF][START_REF] Sugawara | Unconventional Charge-Density-Wave Transition in Monolayer 1T-TiSe2[END_REF] and theoretical [START_REF] Guster | First principles analysis of the CDW instability of single-layer 1T-TiSe2 and its evolution with charge carrier density[END_REF] works focus on the thickness dependence of CDW and the possible unconventional behavior in the ultimate 2D limit of a single layer TiSe2.

On the other hand, the other Ti dichalcogenides namely TiS2 and TiTe2 did not show any clear evidence until very recently when a CDW state was reported only for 1 ML-thin TiTe2 at temperatures lower than 92 K [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] . It is surprising that the CDW in TiTe2 was found to be totally suppressed for films thicker than 1 ML [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] , unlike the case of other TMDs where 1 ML and bulklike films both make the transition to a CDW at nearly the same temperature.

The interest about TiTe2 is continuously increasing in view of theoretical predictions [START_REF] Zhang | Series of topological phase transitions in TiTe2 under strain[END_REF] and more recent experimental evidence [START_REF] Zhang | Pressure-induced topological phase transitions and structural transition in 1T-TiTe2 single crystal[END_REF] about pressure induced topological phase transitions in TiTe2. The possibility to also manipulate superconductivity by external pressure as predicted [START_REF] Xiao | Manipulating superconductivity of 1T-TiTe2 by high pressure[END_REF] and more recently evidenced [START_REF] Dutta | Pressure-induced superconductivity in semimetallic 1T-TiTe2 and its persistence upon decompression[END_REF] in bulk TiTe2 creates the prospect to explore the emergence of topological superconductivity in this material. In the latter work [START_REF] Dutta | Pressure-induced superconductivity in semimetallic 1T-TiTe2 and its persistence upon decompression[END_REF] it has been shown that under non-hydrostatic pressure, a CDW-like state with estimated transition temperature above room temperature (RT) appears in bulk TiTe2 at around 0.5-1.8 GPa. These results call for a reexamination of the possibility to obtain a CDW in multilayer TiTe2 and indeed at RT with good potential for real world applications utilizing the properties of the CDW state [START_REF] Liu | A chargedensity-wave oscillator based on an integrated tantalum disulfide-boron nitride-graphene device operating at room temperature[END_REF][START_REF] Vaskivskyi | Fast electronic resistance switching involving hidden charge density wave states[END_REF][START_REF] Hollander | Electrically Driven Reversible Insulator-Metal Phase Transition in 1T-TaS2[END_REF][START_REF] Khan | Anomalous electron transport in back-gated field-effect transistors with TiTe2 semimetal thin-film channels[END_REF] .

In this paper it is shown that multilayer films (50 ML ~ 32 nm), as well as single layer TiTe2 epitaxially grown on InAs(111)/Si(111) substrates by molecular beam epitaxy exhibit, in ambient pressure conditions, a CDW distortion at room temperature which is sustained up to higher temperatures, at least 400 o C as evidenced by reflection high energy electron diffraction (RHEED) (Fig. S1). The results are explained in terms of anisotropic strain imposed by the substrate.

Stability of 1T TiTe2: First principles calculations

1T-TiTe2 belongs to P3 ̅ m1 space group with trigonal symmetry assuming an octrahedral structure. The stability of 1 ML and bulk 1T TiTe2 structure is studied by DFT and results on phonon dispersion are presented in Fig. 1a and S2. Employing Perdew-Burke-Ernzerhof (PBE) functionals, the 1 ML 1T structure is predicted to remain stable since only a small phonon softening at M occurs (Fig. S2) similar to previous report [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] . However, if Heyd-Scuseria-Ernzerhof (HSE) hybrid functionals are used (see Methods) a large Kohn anomaly is predicted in one of the acoustic branches near M (Fig. 1a), implying a large 1T structural instability which could induce a transition to a lattice with lower symmetry. This is further supported by DFT/HSE energy minimization calculations in a 2×2 supercell predicting that the structure relaxes to a periodically distorted lattice as shown in Fig. 1b which is more stable than the 1T structure by 86 meV per unit cell. In the distorted lattice, one out of four Ti atoms in the 2×2 unit cell (TiI in Fig. 1b) remains fixed at the same position as in undistorted 1T structure forming six equal bonds (2.767 Å) with the top and bottom Te atoms. On the other hand, the other three Ti atoms (TiII, in Fig. 1b) in the 2×2 cell are displaced from the ideal 1T lattice positions by 0.021a (a is the inplane lattice constant); Te atoms are also slightly displaced by 0.002a. Due to the TiII atom displacements, two short and two long bonds with lengths 2.701 Å and 2.840 Å, respectively, are formed resulting in the configuration of Fig. 1b which is similar to the 1 ML TiSe2 2×2 commensurate CDW structure as previously reported 7,[START_REF] Guster | First principles analysis of the CDW instability of single-layer 1T-TiSe2 and its evolution with charge carrier density[END_REF] . In accordance with phonon calculations, when PBE functionals are used, energy minimization results in the 1T TiTe2 as the most stable structure. The structure is predicted to be essentially planar (confirmed by electron microscopy, Fig. 3 below) with no significant buckling and with all atoms lying on the lattice planes as in the undistorted 1T TiTe2.

The distorted lattice calculated by HSE hybrid functionals is found to be more stable by 90 meV per unit cell than the 1T, also in the case where enlarged lattice parameters are considered as a result of tensile strain from the substrate as evidenced by GIXD measurements (see in the text below). The theoretical atomic fractional coordinates for the relaxed and strained 1 ML TiTe2 are given in tables S1 and S2, respectively.

CDW order at the surface probed by in-situ scanning tunneling microscopy (STM)

The surface structure probed by atomic resolution STM at RT and their corresponding fast Fourier transform (FFT) patterns are presented in Fig. 2. The 1 ML TiTe2 STM image and 2×2 FFT (Fig. 2a,b) are very similar to data previously reported for relaxed 1 ML TiTe2 at 4.2 K and attributed to CDW superstructure [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] . The main difference is that in the present work the Moiré pattern attributed to the interaction with the InAs substrate is also visible. As expected, for the thicker 50 ML film, the Moiré pattern is absent (Fig. 2c,d) since the top layer is not in direct contact with the substrate. Bright intensity modulation in a hexagonal formation (white arrows) in the real image indicates the doubling of periodicity (Fig. 2c) which however cannot be attributed to height differences between top Te atoms given that the Te plane is essentially flat as already discussed above and evidenced by STEM (see Fig. 3 below). Instead, this STM intensity modulation could be explained as charge density modulation associated with a transition to the low symmetry structure as depicted in Fig. 1b. It should be noted that the 2×2 superstructure signal in the FFT of the 50 ML (Fig. 2d) is significantly stronger compared to 1 ML film (Fig. 2b) which is relatively weak and diffuse, in agreement with the weak (and barely visible) superstructure streaks in RHEED (Fig S1). This indicates a much stronger CDW formation in the multilayer TiTe2 compared to 1 ML.

Direct observation of CDW structure in multilayer TiTe2 by high resolution scanning tunneling spectroscopy (HR-STEM).

STEM and nanobeam precession electron diffraction (N-PED) measurements are performed on the 50 ML-thick film and the main results are presented in Fig. 3 and in Fig. S3. The highresolution STEM image along the [21 ̅ 1 ̅ 0] zone axis in Fig. 3a shows a flat structure with great similarity to 1T TiTe2. However, a closer look reveals that the intensity of the Ti columns in Fig. 3a,c,d is modulated between sharp/intense (green arrows) and weak/diffuse (red arrows) essentially doubling the periodicity parallel to the layer. The same pattern appears along the vertical direction (c-axis) where Ti column brightness is modulated between high and low values in two consecutive TiTe2 layers, thus doubling the periodicity along c axis too, resulting in an overall 2×2×2 periodicity in the film. This periodicity is also evident in the FFT (Fig. 3b) and the N-PED pattern (Fig. 3e) since weak superstructure spots, indicated by vertical blue thick arrows appear in between the main ones. The real image directly points to a CDW ABA type of stacking as shown in Fig. 3f where the aligned Ti atom columns A produce the sharp and bright Ti atom fringes (Fig. 3c,d) while the slightly misaligned atom columns B produce the weaker and diffuse Ti atom fringes (Fig. 3c,d Considering also the data in Fig. S4, it can be said that, more generally, the half integer l RSMs are dominated by the CDW superstructure peaks and the integer l RSMs are dominated by the 1T TiTe2 Bragg peaks. A diffuse elongated scattering pattern more visible in the half integer l-planes (Fig. 4 andS4) is attributed to thermal phonon diffuse scattering originating from Si substrate 23. The reciprocal space images of the h0l and hhl planes are shown in Figs 4c and4d, respectively. 1T TiTe2 Bragg peaks have integer hl values, while CDW superstructure peaks have half-integer hl values, clearly defining the ×2 pattern along the in-and-out-of plane directions. The finite scattering intensity in the vertical direction (along ΓA) which is more pronounced in the h0l plane, is attributed to Si peak contributions (green arrows) and to finite sample thickness resulting in finite intensity in the half-integer l RSMs (Fig. 4a andS4b,d) at the Γ points where it should normally be absent. The h0l pattern is compatible with the FFT and N-PED patterns observed in STEM (Fig. 3b,e).

Selected experimental l-scans are presented in Fig. 4e showing 1T TiTe2 (left panel) peaking at integer l and CDW superstructure peaks (right panel) at half integer l. Order of magnitude lower intensity peaks at integer l (shown by the vertical arrows in right panel of Fig. 4e), correspond to the M point weak superstructure spots as discussed above. In Fig. 4f the lscans are simulated [START_REF] Drnec | Integration techniques for surface X-ray diffraction data obtained with a two-dimensional detector[END_REF][START_REF] Vlieg | ROD: a program for surface X-ray crystallography[END_REF] (for details see SI) using the CDW model with the ABA type of stacking (tables S1 and S2). A qualitative agreement with the experimental data of Fig. 4e can be seen, since the main features are reproduced at the correct hkl positions within an order of magnitude agreement in intensity. The out of plane ×2 periodicity manifested by the strongest superstructure peaks at half-integer hkl is a direct consequence of the breaking of translation symmetry along caxis in the ABA stacking configuration considered here. Markedly, AAA CDW stacking model fails to reproduce the experimentally observed ×2 periodicity (data not shown).

The 2×2×2 CDW reconstruction is a characteristic of our epitaxial TiTe2 films since bulk TiTe2 shows a stable 1T structure where superstructure peaks at the L points of half integer l RSMs are totally absent as can be seen from Fig. S5. Also by a direct comparison, it can be inferred that our epitaxial films exhibiting sharp diffraction spots are of superior quality compared to bulk (Fig. S5). Our films with low in-plane mosaic spread (see Fig. S6 and related discussion) are very well aligned with the substrate following the in-plane epitaxial relationship:

InAs [112 ̅ ](111) //TiTe2 [101 ̅ 0](0001), despite the large lattice mismatch (~ -12%). The 1×1 TiTe2 in-plane a and out-of-plane c lattice parameters are accurately measured to be a=3.842 (±0.005) Å in 1ML TiTe2 and a=3.862(±0.003) Å, c =6.349 (±0.008) Å in 50 ML TiTe2. These values differ substantially from the experimentally determined free-standing bulk values a= 3.777 Å, c = 6.498 Å [START_REF] Dutta | Pressure-induced superconductivity in semimetallic 1T-TiTe2 and its persistence upon decompression[END_REF][START_REF] Arnaud | Etude comparative des composés TiX2 (X = S, Se, Te). Structures de TiTe2 et TiSeTe[END_REF] . For the 50 ML film, the in-plane parameter a is enlarged by 2.2% while the out of plane parameter c is reduced by 2.3% indicative of epitaxial tensile strain in the film imposed by the substrate. Similarly, an epitaxial tensile strain of 1.7% can be inferred for the 1 ML TiTe2 film.

Modifications in the electronic band structure imaged by in situ angle resolved photoemission spectroscopy (ARPES)

The electronic band structure of 1 ML (Fig. 5a) and 50 ML TiTe2 (Fig. 5b) are imaged by in-situ ARPES. They reveal a semimetalic material with two Te 5p valence bands crossing EF near Γ, which overlap with Ti 3d conduction bands at M(L) [START_REF] Rossnagel | Three-dimensional Fermi surface determination by angle-resolved photoelectron spectroscopy[END_REF] by 0.2-0.3 eV, as expected from DFT calculations [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] and experimental data [START_REF] Rossnagel | Three-dimensional Fermi surface determination by angle-resolved photoelectron spectroscopy[END_REF][START_REF] De Boer | Titanium ditelluride: Band structure, photoemission, and electrical and magnetic properties[END_REF][START_REF] Claessen | Complete band-structure determination of the quasi-two-dimensional Fermi-liquid reference compound TiTe2[END_REF] of undistorted 1T TiTe2, suggesting that the PLD/CDW formation in our epitaxial films does not have a large influence on the electronic band structure.

In contrast to common wisdom that a CDW should open a gap at the Fermi level (Ef), no such gap is observed here. This is not in contradiction though with other CDW cases, as for example in 1ML TiTe2 [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] where valence and conduction bands cross at Ef, also in TiSe2 [START_REF] Chen | Charge density wave transition in single-layer titanium diselenide[END_REF][START_REF] Chen | Dimensional Effects on the Charge Density Waves in Ultrathin Films of TiSe2[END_REF][START_REF] Sugawara | Unconventional Charge-Density-Wave Transition in Monolayer 1T-TiSe2[END_REF] where despite evidence for a small gap, Ef often lies within one of the bands having finite density of states. It is worth noticing that the multilayer TiTe2 (Fig. 5b) with a stronger CDW structure compared to 1 ML as evidenced from Fig. 2, shows a lowering of one of the bands near Γ, about 0.25 eV below Ef (blue marks in Fig. 5b). This indicates a reduction of the density of states near Ef corresponding to a partial gap opening (or a pseudogap) which could result in electronic energy gain, counterbalancing the losses associated with the PLD/CDW lattice distortions.

Notably, there is a weak and diffuse (blur) intensity near Ef distributed along the ΓM directions (E-k and contour plots in Figs. 5 andS7), which could be considered as a dispersion-less band (DB), not predicted by theory. Similar dispersion-less intensity (blur) indicating localized conduction electrons has been observed previously in TiTe2 [START_REF] De Boer | Titanium ditelluride: Band structure, photoemission, and electrical and magnetic properties[END_REF] and has been attributed to polaron effects associated with CDW order although it was not reproduced in later works [START_REF] Claessen | Complete band-structure determination of the quasi-two-dimensional Fermi-liquid reference compound TiTe2[END_REF] .

It should be noted that overlapping of conduction and valence band was expected at M as a result of periodicity doubling and associated zone folding as previously reported [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] at T <92K in synchrotron ARPES experiments. However such overlapping is not observed here possibly due to small spectral weights of the Te-p valence bands at M yielding signals beyond the sensitivity and resolution of our lab-based room temperature ARPES.

Discussion and conclusions

After several decades of research on CDW in 2D TMDs and the debate about the associated mechanisms, the common wisdom is that the softening of an acoustic phonon 8,9,[START_REF] Yoshida | Electron Lattice Interaction and Lattice Instability of 1T-TiSe2[END_REF] and the electron-phonon interaction are involved in the formation of a PLD and an associated CDW. In our multilayer TiTe2 the stronger intensity superstructure spots are obtained at the L points of the Brillouin zone (Fig. 4a) implying that softening of an acoustic phonon with wavevector 𝑞 ⃗ =ΓL ⃗⃗⃗⃗⃗ is at the heart of the commensurate 2×2×2 PLD/CDW formation in multilayer TiTe2 at RT, similar to the TiSe2 case 7 , while hints of localized conduction electrons at the Fermi surface from ARPES measurements (Fig. 5) indicate polaronic effects [START_REF] Rossnagel | Three-dimensional Fermi surface determination by angle-resolved photoelectron spectroscopy[END_REF] likely associated with the CDW.

Softening of an acoustic phonon with 𝑞 ⃗ =ΓM ⃗⃗⃗⃗⃗⃗⃗ as also predicted by theory (Fig. 1a) and evidenced by the very weak superstructure peaks at M (Fig. 4b,e) in GIXD could be responsible for the CDW observed at RT in 1 ML TiTe2 (Fig. 2a,b).

In free-standing TiTe2 multilayer films with relaxed lattice parameters [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] and a sizeable v.d. Waals gap, the weakly coupled Te-Ti-Te (tri)layers do not favor the propagation of the CDW-driving 𝑞 ⃗ =ΓL ⃗⃗⃗⃗⃗ phonons through the layers, probably explaining the absence of CDW in multilayer films [START_REF] Chen | Emergence of charge density waves and a pseudogap in single-layer TiTe2[END_REF] . In our epitaxial films though the situation is different, because a significant interaction with the substrate is established producing an anisotropic strain such that the in-plane parameter is enlarged and the out of plane parameter along the c-axis is compressed by 2.3% as accurately measured by synchrotron GIXD, yielding an estimated stress ~0.6-0.8 GPa if a measured [START_REF] Rajaji | Structural, vibrational, and electrical properties of 1T-TiTe2 under hydrostatic pressure: Experiments and theory[END_REF] bulk modulus B ~ 0.28-0.40 GPa is assumed. The compression results in an appreciable reduction of van der Waals gap from ~ 3.2 Å 18 to the value of 3.1 Å that we measured by STEM bringing the Te-Ti-Te trilayers closer between each other thus enhancing interlayer Te-Te interactions, signifying a transition from a 2D to a 3D-like TiTe2 crystal. We anticipate that this transition to a 3D structure facilitates the propagation of the important phonon modes along ΓL which are required to propagate across tightly coupled TiTe2 layers in thick films or bulk material to provide the source of instability that drives the CDW distortion at RT which is sustained up to 400 o C.

The analysis based on the anisotropic strain correlates well with predictions [START_REF] Fu | Modulation of doping and biaxial strain on the transition temperature of the charge density wave transition in 1T-TiSe2[END_REF] in TiSe2 stating that the CDW critical temperature is enhanced under biaxial tensile (stretching) strain (compression along the c-axis). Our analysis correlates also nicely with the recent experimental results [START_REF] Dutta | Pressure-induced superconductivity in semimetallic 1T-TiTe2 and its persistence upon decompression[END_REF] reporting the emergence of CDW in a TiTe2 multilayer under 1.8 GPa non-hydrostatic pressure (effectively uniaxial compression along the c-axis as in our films) also providing hints of a re-entrant 2×2 CDW with an estimated critical temperature above RT. The latter appears during the decompression phase around 0.5 GPa, which is similar in magnitude to the epitaxial stress exerted on our films (0.6-0.8 GPa, see above).

Methods:

Growth The In-terminated-InAs(111)/Si(111) wafers prepared by MOCVD are used as the substrate materials. The substrate are chemically treated in a solution of HCl(5N) in isopropyl alcohol and subsequent dipping in isopropyl alcohol, followed by an annealing step at 400 C in UHV until a 2×2 superstructure is obtained by RHEED, indicating a clean, oxide-free surface.

The molecular beam epitaxy (MBE) of thin TiTe2 films was performed in an ultra-high vacuum (UHV) chamber at 400 C. High purity metal Ti was evaporated from an e-gun evaporator whereas Te was evaporated from effusion cells. TiTe2 epilayers were grown under Te rich conditions with a growth rate of 1 ML/30sec.

Scanning Tunneling Microscopy

The surface morphology of the TiTe2 epilayers on InAs(111) is examined by in-situ RT Omicron UHV-STM. The scanning conditions are It =500 nA, Vb =+0.05 V, loop gain 1.7. Vb in our system is applied on the sample so that positive Vb means empty state STM images.

Electron Microscopy. STEM and N-PED experiments were performed using a probe-corrected FEI Titan Themis S/TEM operated at 200 keV. High angle annular dark field STEM images were collected using a convergence semi-angle of 18 mrad with the annular detector collecting scattering above 65 mrad. N-PED was collected using the FEI Epsilon acquisition software, and more details about the method can be found by these works [START_REF] Béché | Improved precision in strain measurement using nanobeam electron diffraction[END_REF][START_REF] Cooper | Combining 2 nm spatial resolution and 0.02% precision for deformation mapping of semiconductor specimens in a transmission electron microscope by precession electron diffraction[END_REF] . The N-PED data was acquired using a convergence semi-angle of 1.5 mrad, a precession angle of 0.3 degrees, and precession frequency of 200 Hz. Electron diffraction patterns were simulated using the JEMS software package. The specimen was prepared by focused ion beam lift-out protocol using a FEI dualbeam Strata 400S.

Synchrotron GIXD Measurements: Diffraction measurements were performed at the European Synchrotron Radiation Facility by means of the UHV-MBE CVD diffractometer installed at the BM32 CRG/IF, optimized for GIXD, and ID28 diffraction side station. The experimental setup energy and incident angle were set at 11 keV (1.13 Å) and 0.2°, the latter set slightly below the critical angle value for total reflection in order to enhance the 2D film signal while minimizing the background. During the measurements the sample was maintained under UHV (base pressure: 3x10 -3 mbar) in the UHV chamber that is coupled to the diffractometer. ID28 data were collected with the fixed incidence angle 1 deg. at 17.8 keV (0.697 Å) on the sample rotating around its normal; scattering intensity was collected by PILATUS3X 1M detector in shutterless mode with angular step of 0.25 deg. CrysAlis package 35 was employed for refinement of experiment geometry, high qualiry reciprocal space reconstructions were produced with locally developed software.

Photoelectron spectroscopy. ARPES measurements were carried out with a PHOIBOS 100 (SPECS) hemispherical analyzer, made of μ-metal at a pass energy of 15 eV, using a 2D CCD detector and a He plasma source with He I excitation energy at 21.22 eV. The energy resolution of the system was better than 40 meV, although the measurement resolution is limited by thermal excitation processes to around 100 meV.

First Principles Calculations. Density Functional Theory calculations were performed using the Vienna Ab initio Simulation Package (VASP) [START_REF] Kresse | Ab initio molecular dynamics for open-shell transition metals[END_REF][START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF] and projector-augmented waves [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] . The Heyd-Scuseria-Ernzerhof (HSE06) hybrid functionals [START_REF] Heyd | Hybrid functionals based on a screened Coulomb potential[END_REF] were used as the exchange correlation potential to determine the atomic displacements in the (2×2) CDW TiTe2 structure, total energy and phonon dispersion calculations. Atomic positions were optimized until the residual forces were less than 0.005 eV•Å -1 . Also, a vacuum 15 Å was used to avoid interaction between the periodically repeated layers. The plane wave kinetic energy cutoff was set at 500 eV, employing a 16×16×1 (8×8×1) k-mesh for the normal (CDW) phase. Perdew-Burke-Ernzerhof functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] was also used for the phonon dispersion calculation of the (1×1) 1T-TiTe2 structure. Phonon calculations were performed using the Phonopy package [START_REF] Togo | First principles phonon calculations in materials science[END_REF] .

Phonon dispersion using PBE functionals

The phonon dispersion 2 is also calculated using PBE functionals 3 to determine the stability of the (1×1) 1T TiTe2 monolayer (Fig. S2) assuming enlarged lattice parameter a as determined by XRD (see main text). In contrast to the phonon calculations using HSE hybrid functionals, no imaginary frequencies are predicted, indicating that the structure should remain stable since only a small phonon softening at M occurs. These results are in line with previously published data 4 for the 1 ML TiTe2 with relaxed lattice parameter a, using PBE functionals. This indicates that taking into account electron-correlation effects through HSE functionals is essential for the prediction of a Kohn anomaly at M leading to structural instabilities. The l=0 and l=2 RSMs, similar to the l=1 RSM (Fig. 4b) are dominated by the 1T TiTe2 Bragg peaks at the Γ points of the BZ having all indices hkl integer. In the l=2 RSM very weak superstructure spots at M points (marked by yellow arrows) are seen similar to the l=1 RSM in Fig. 4b, however their appearance is not confirmed for l=0 since the scattered beam is below the sample horizon and thus attenuated by few orders of magnitude. The half integer l RSMs (1.5 and 2.5) are dominated by the CDW superstructure spots, located at the L points of the BZ and having at least one in-plane index (k or h) half integer, similar to l=0.5 RSM in Fig. 4a (main text).

Simulations of Crystal Trancation Rods (l-scans)

The structure factor amplitude along crystal truncation rods (CTRs) is calculated as a function of the reciprocal space coordinate l for (h,k) integer (left) and (h,k) half-integer (right) out of plane directions for the model described in, with Debye-Waller factors set equal to zero. When both h and k are integer, only Bragg peaks at integer l coordinates are visible, whereas superstructure peaks are not observed in between. For (h,k) half-integer coordinates instead, peaks occur at half-integer values of l as consequence of the ABA stacking of layers) with intensity of at least one order of magnitude lower than Bragg peaks in the left figure. Moreover, along the same reconstructed directions, additional weaker peaks appear also at integer l, in agreement with experimental observations at M points (see main text). The calculations have been performed with PyRod, software developed at BM32 (ESRF) and based on the well known ROD software 4,5 for surface diffraction data analysis and processing.

  Figure 4a,b show in-plane reciprocal space maps (RSM) for the 50 ML InAs/TiTe2 sample, for

FiguresFig. 1 .

 1 Figures

Fig. 2 .

 2 Fig. 2. Surface ordering of epitaxial TiTe2 by in-situ STM at room temperature. Images are taken at

Fig. 3 :

 3 Fig. 3: HR STEM and diffraction pattern simulation. (a) STEM image of TiTe2 film along the [21 ̅ 1 ̅ 0]

Fig. 4

 4 Fig. 4 Reciprocal space mapping of CDW 50 ML TiTe2. (a) RSM for l=0.5 (a) and l= 1 (b). The

Fig 5 :

 5 Fig 5: Electronic band structure imaging of CDW TiTe2 by by in-situ ARPES. Energy dispersion E-k

Fig. S2 .

 S2 Fig. S2. Phonon frequency dispersion of 1ML 1T TiTe2 along high symmetry points in

Fig. S3 :

 S3 Fig. S3: STEM measurements. (a) STEM image of TiTe2 film on InAs (111) substrate. (b)

7 .

 7 Diffraction from bulk 1T TiTe2 crystals 200-300 μm-thick bulk TiTe2 crystals fabricated by chemical vapor transport by HQ Graphene are investigated by GIXD (Fig. S6) to compare with our epitaxial films. All integer l RSMs show the 1T TiTe2 bragg peaks at the center (Γ point) of the BZ as expected. In contrast to the epitaxial films (Figs 4 and S4), the superstructure spots expected at the L points of the BZ are absent from all half-integrer l RSMs (Fig S6) in the bulk TiTe2. This means that bulk TiTe2 is stable in the 1T structure and that the CDW is characteristic of epitaxial films due to the influence from the substrate.

Fig. S5 .

 S5 Fig. S5. Reciprocal Space Mapping of bulk 1T TiTe2 for different l values. The yellow

Acknowledgements

This work was financially supported by the LANEF Chair of Excellence program of U. Grenoble Alpes and CEA (A.D.), the Greek State Scholarships Foundation (IKY) Program for the strengthening of postdoctoral research (P.T.), the project MIS 5002567, co-financed by Greece and the EU (European Regional Development Fund) (D.T), the French state fundsANR-10-LABX-51-01 (Labex LANEF du Programmed'Investissements d'Avenir) and Equipex ANR-11-EQPX-0010 (G.R.). The authors thank Dr. P. Rodiere of CNRS/Neel Institute for very useful discussions on soft phonon mediated X-Ray scattering and Dr. T. Baron of Université Grenoble Alpes, CNRS, CEA/Leti Minatec, LTM for providing us with the InAs/Si substrates.

Author contributions

Sotirios Fragkos performed DFT calculations and was involved in ARPES, Roberto Sant, Alexei Bosak, and Gilles Renaud performed the GIXD, Carlos Alvarez and Hanako Okuno performed STEM, P. Tsipas performed the MBE growth, RHEED and STM, D. Tsoutsou performed ARPES and A. Dimoulas was involved in DFT, MBE growth, STM and ARPES

Competing Interests

The authors declare no competing interests.

Materials & Correspondence

Should be addressed to Athanasios Dimoulas

SUPPLEMENTARY INFORMATION 1. In-situ reflection high energy electron diffraction

The RHEED patterns recorded along the InAs [110] and [112 ̅ ] azimuths for the MBE-grown 1 ML and 50 ML TiTe2 as well as for the bare InAs substrate are presented in Figure S1. After insitu substrate preparation, a clean, oxygen and carbon-free In-terminated InAs(111) surface was obtained, as evidenced from the recorded 2×2 reconstruction, which is attributed to In vacancies 1 . The 2×2 reconstruction is an indication that the InAs surface is clean and ready for the TiTe2 growth. The streaky RHEED patterns of both 1 ML and 50 ML TiTe2 films indicate very good surface atomic ordering of the epilayers and rotational commensurability. Additional satellite peaks are observed in the vicinity of the main (1×1) streaks of the 1 ML TiTe2 better visible in off-angle patterns (not shown here), which originate from Moiré interference fringes as described in more detail in the main text in the section discussing STM data (Fig. 2a andb). In the 1ML TiTe2 RHEED spectrum, a very faint and broad ×2 reconstruction is also observed mainly in the [112 ̅ ]InAs azimuth which is associated to the similarly weak and diffuse STM/ FFT pattern (Fig. 2b).

The same ×2 pattern is clearly visible in both azimuths in the thick (50 ML) film which correlates nicely with the more intense and sharp 2×2 superstructure in the STM/FFT of the 50 ML TiTe2 in Fig. 2d S2. DFT calculated atomic positions of (2×2) CDW-TiTe2 1ML unit cell expressed in fractions of the experimental lattice parameters. The latter are deduced from synchrotron GIXD in the present work and are found to be enlarged compared to free-standing values determined from DFT, due to the influence of the substrate.

RSM for l= 0 and higher diffraction orders of multilayer TiTe2 CDW structure

The RSM for l=0.5 and l=1 have already been described in detail in Fig. 4 a, b in the main text.

Here we complementary give the l=0 (a), and higher order diffractions l=1.5 (b), l=2 (c) and l=2.5 (d). symmetry in the first BZ (unprimed) and in neighboring BZs (primed and double-primed).

GIXD-Rocking curves and mosaicity

Furthermore, low in-plane mosaic spread of 1.


Fanny Bocquentin

Hamoudi Khalaily

Nicolas Samuelian

Omry Barzilai

Gaëlle Le Dosseur

Kolska Liora

Aline Horwitz

Emery-Barbier

Renewed excavation on the PPNB site of Beisamoun, Hula Basin

à 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.

Introduction

The Hula Basin, located in the northern part of the Jordan Valley, enjoys a moderate Mediterranean climate. This region has played an important role in the prehistoric occupation of the central Levant since its rich ecosystem offered people the possibility of long-term occupation of one area. As part of the central Jordan Valley, it would have served as a conduit for East-West and North-South migration/exchange of populations and/or materials. Thus, the Hula Basin is particularly suitable for exploring environmental exploitation systems and cultural as well as biological transmission over time as well as between groups inhabiting the same region. The site of Beisamoun in the Hula Basin is a key PPNB location that can provide important information on these issues.

The major occupation phase at Beisamoun dates to the Pre-Pottery Neolithic B but Pottery Neolithic and Bronze Age deposits are also present, as attested to by archaeological finds and recent rescue excavations conducted by one of us (H. Kh.).

Since the early 1950s the area of the PPNB site was extensively investigated as fish-breeding ponds were dug into the underlying Neolithic settlement and partly destroyed the archaeological layers. The subsequent drainage of these ponds enabled A. Assaf (Curator of the Mayan Baruch Museum of Prehistory of the Hula Basin) to collect a large PPNB assemblage and to demar-cate the extent of the site to an area of 10 ha. In 1969, a French team supervised by M. Lechevallier (CNRS) began mapping the structures discernible at the bottom of the different ponds, and this was followed in 1972 by a 70 m2 salvage excavation (Lechevallier 1978). This area yielded a rectangular structure with a well preserved plaster floor. It was associated with several graves, numerous artefacts and, among other unexpected findings, two plaster skulls. No radiometric dating was undertaken.

During the summer of 2007, a two-week long exploratory field season was held from 24 August to 7 September in order to test a large area at the north-western part of the site (at the location of the modern fish pond number 11). A surface of 3000 m 2 was mechanically cleaned. Within this area, one sector of 76 m 2 was excavated from 0.20 m to 0.50 m in depth (Sector E) and three long trenches were opened (Trenches 1,2,3).

Sector E Stratigraphy

Four levels were recognized in Sector E:

• Level 0 (from 10 to 30 cm deep) comprises a black, modern sediment very rich in archaeological finds dating to the PPNB. This level is superficial in the southwestern part of the sector but deepens towards the northeastern side.

• Level I is a brown-greyish sediment, rich in finds and contains the uppermost stone structures.

Renewed Excavation of the PPNB Site of Beisamoun, Hula Basin

Fanny Bocquentin1,2 , Hamudi Khalaily3 , Nicolas Samuelian4 , Omry Barzilai5 , Gaelle Le Dosseur 1,2 , Liora Kolska Horwitz 5 and Aline Emery-Barbier 1

• Level Ia is a fine-grained brown-yellowish sediment, relatively poor in finds but has yielded the second level of structures. It is underlain by the earlier occupation layer,

• Level Ib, which was found at the end of the excavation season, under a floor level associated with Layer Ia that covered a large area of the sector. This is a finegrained sediment similar to that found in Layer Ia but very compact and richer in charcoal and ochre specks. Six charcoal samples from Layers I, Ia and Ib have been submitted for dating.

Within these levels, nine loci were recognised and partly excavated (Fig. 1). At least two phases of construction were recognized, and both probably relate to the very end of the PPNB.

From the earlier occupation layer found so far, we have exposed five loci. Three of them are segments of walls:

• Locus 203 is a wall comprising two rows of stones oriented NNW-SSE. This wall was exposed for one meter and then underlies Locus 200.

• Locus 202 is oriented east-west and is made up of stones and mud bricks. It is preserved for 2.4 m in length.

• Locus 207 is a wall of two-courses oriented perpendicularly to the previous one.

• Locus 205 is a sort of platform made of a 20 cm thick, hard calcified white material showing a stair pattern on its west side and a curved pattern on the north-east edge (Figs. 1 and2). A bloc-sample of this locus was taken and is being examined by F. Berna (Boston University). It may represent the remains of a poorly preserved plaster floor. Its northern limit is strictly parallel to wall 202, probably part of the same structure.

• Lastly, Locus 208 is a stony area strewn with numerous items discovered in a flat position, on or between the stones. They are the testimony to an unplastered floor level.

Four loci are part of a second phase of occupation.

• Locus 200 is a wall of two-ranges, oriented NW-SE, and currently exposed only for 5.30 m (Figs. 1 and3).

• Locus 201 is an oval structure made of two or three rows of stones adjoining wall 200.

• Locus 204 comprises two rows of stones which are visible on the surface of the excavated area.

• Lastly Locus 207 is an oval pit that penetrates Layers I and Ia, destroying part of floor 208.

Both phases of construction are directly on top of each other, suggesting that there was no major interruption in occupation between them. We have investigated the later period less than the earlier one. This may appear as a paradox but, in fact, the later structures are preserved only in the SW corner of Sector E, an area enlarged only during the second week of the excavation.

Tools

Preliminary observations of the flint assemblage from the 2007 season at Beisamoun (more than 7000 artefacts) indicate that the flint industry is domestic and most of the knapping activities took place on-site, which is typical of large, permanent sites. This conclusion is based on the variety of the lithic material (cores, waste, debris and tools) and of the chipped stone technologies (blade, bladelet, flake and bifacial tools). Despite the stratigraphic distinctions, the flint assemblage seems to be homogenous in raw material acquisition and mode of production. Technologically, the high frequency of bidirectional blade technology that is characteristic of the PPNB (e.g. [START_REF] Abbès | Les outillages néolithiques en Syrie du Nord : méthode de débitage et gestion laminaire durant le PPNB[END_REF]Khalaily 2006) is noteworthy. The assemblage also includes a classic naviform blade technology made on high quality purple flint. As for the tools, two types of sickle blades are present: sickle blades that display fine denticulation, mainly on the ventral face and high luster on their working edges, together with deep denticulated sickle blades that have been shaped on bidirectional products; a few unpolished axes and long Amuq points that have been shaped by pressure retouching were also recovered.

Based on these techno-typological observations, it is possible to attribute the assemblage to the later stage of the PPNB period. This attribution is also based on the appearance of bladelet technology and standardized tools such as deep denticulated sickle blades and pressure retouched Amuq points that are characteristic of the Early PN Yarmukian culture [START_REF] Rollefson | The Neolithic Period[END_REF][START_REF] Gopher | Trends in Sickle Blade Production in the Neolithic of the Hula Valley, Israel[END_REF].

Knapping industries are also represented by finegrained basalt cores and flakes as well as obsidian pieces. Fourteen of them were collected during this season (three blade fragments, four bladelet fragments and seven debris). All belong, a priori, to the same raw material, a homogeneous, grey translucent obsidian with small black inclusions.

Fragments of grinding stones, numerous mullers, a few bowls, and discoidal items made out of limestone or vesicular basalt, were also found.

Twenty-three artefacts manufactured from animal elements (bone and antler) were found. Most of them derive from Level Ia. Two main tool categories were recognised in the sample -pointed tools and cutting ones. Most are made of bone. Metapodials were selected for awls and ribs for flat knives. Obvious technical reasons have guided these choices. Several techniques have been recognized -percussion, sawing, shaving, grindingwhile at least two methods of debitage were practiced: fracturing and bipartition. These technical choices are usual in PPNB contexts, and some are inherited from the Natufian (Le [START_REF] Dosseur | La fabrication et l'entretien des poinçons en os au Levant Sud : investissement technique et règles de production du Natoufien au Néolithique Précéramique B récent[END_REF]. A bigger sample is needed to discuss the possible unique features of the Beisamoun bone tool kit.

Fauna

The new excavations at Beisamoun have yielded a small (ca. 400 bones from Layers I and Ia) but interesting sample of well preserved animal bones, the majority of which comprise remains of what have tentatively been identified as wild species: Bos cf. primigenius, Sus scrofa fer., as well as Gazella gazella. The caprine (Ovis/Capra) remains comprise a minority of the identified finds (Table 1). The majority of caprine bones are undetermined as to species, but goats, and no sheep, have been positively identified to date. At least one slightly twisted goat horn core of an animal, tentatively identified as Capra cf. hircus, was found in Layer I, but the status of the caprines from the other layers cannot be absolutely determined at this stage of the analysis.

The low frequency of caprine remains found in this sample is unusual for PPNB assemblages from the south- ern Levant, since by the MPPNB caprines comprise 40-50% or even more in other Jordan Valley assemblages [START_REF] Horwitz | Counting Cattle: Trends in Neolithic Bos Frequencies from the Southern Levant[END_REF]. The high cattle frequencies in the new Beisamoun sample closely resembles the PPNC assemblages from Hagoshrim and Tel Te'o in the Jordan Valley, where cattle frequencies comprise 30% or more of the identified fauna relative to frequencies of caprines and gazelles, and the LPPNB sites in this region with frequencies of circa 25% cattle [START_REF] Horwitz | Counting Cattle: Trends in Neolithic Bos Frequencies from the Southern Levant[END_REF]. These researchers have suggested that such high proportions indicate that by the PPNC in the southern Levant, cattle were domesticated or in the course of this process. A cursory examination indicates that in size and robusticity the Beisamoun cattle are similar to wild forms, but this is a feature that needs to be examined in greater detail.

The new Beisamoun assemblage differs in the range of taxa represented and their relative frequencies compared to the samples previously recovered from this locality and analyzed by [START_REF] Davis | Etude de la faune[END_REF] and [START_REF] Ducos | La faune de Beisamoun dans les collections du Musée Préhistorique de la Vallée du Houleh[END_REF] (Table 1). These differences may stem from the smaller size of the new sample, differences in past human activities undertaken in this particular part of the site, or even chronological differences between the assemblages.

Botanical Remains

Soil samples from Layers I, Ia and Ib were checked for botanical remains (Fig. 4). Within Layers I and Ia, paly-nomorphs reflect a humid and poorly oxygenized milieu of fossilization of the different identified micro-remains. This might be a result of recent botanical pollution during the time of the fish-breeding ponds. In contrast, Layer Ib shows a spectrum typical of dry conditions, with a high proportion of Chenopodiacae in the pollen spectrum. However, very few taxa were identified in this layer and this result may be due to differential preservation of the fossils rather than a testimony to a dry climate. These first botanical samples taken at Beisamoun cannot as yet be interpreted in terms of palaeoenvironment, but this issue will be pursued as micro-remains preservation improves as the excavation of the archaeological layers deepens.

Trenches

Besides the exploration of Sector E, we also opened three long trenches (0.6 m wide and 10-15 m long) using a mechanical excavator. Two of them are oriented northsouth (Trenches 1 and 2) and the third one is oriented east-west (Trench 3). They have provided a long archaeological sequence with at least three different stratigraphic units preserved from 0.20 m to 1.20 m depth. This occupation sequence seems mainly to be associated with the Late PPNB. Like Sector E, the sub-units I, Ia, and Ib are also present in the trenches, except in the east side of Trench 3 where the stratigraphic sequence differs. Four thick walls were found in Trench 2 parallel to each other, oriented east to west. In Trench 3 one wall and two very well-preserved plaster floors were discovered.

At the base of all trenches, a dense black organic layer was found, probably the deposit of the now extinct Hula Lake.

Conclusion

Recent excavations in the Levant have reshaped our knowledge of the Late PPNB, in particular the megasite phenomenon. However, in the western part of the Jordan Valley and in the Jordan Valley itself, the end of the PPNB is still unclear. The results of the first season of renewed excavations at Beisamoun demonstrate the preservation of a promising and undisturbed part of a LPPNB village. This is an excellent opportunity to clarify the issue of the end of the PPNB in this central area of the Levant.

Fig. 1

 1 Fig.1 General map of the 2007 season at Beisamoun, Sector E (drawing: Marjolaine Barazani, CRFJ).

Fig. 2

 2 Fig. 2 View of Loci 202 and 205 from the east, Sector E.

Fig. 3

 3 Fig. 3 View of Loci 200 and 201 from the south, Sector E.

Fig. 4

 4 Fig. 4 Percentage of botanical micro-remains in three soil samples from Sector E, Layers I, Ia and Ib.

UMR 7041 CNRS, Nanterre <fanny.bocquentin@mae.u-paris10.fr>

Centre de Recherche Français de Jérusalem

Israel Antiquities Authority <hamudi@israntique.org.il>

Institut National de Recherches Archéologiques Préventives, France

Hebrew University of Jerusalem

Acknowledgements. The new excavations at Beisamoun were funded by the Irene Levi Sala Care Archaeological Foundation.


Grégoire Misguich

Frédéric Mila

Quantum dimer model on the triangular lattice: Semiclassical and variational approaches to vison dispersion and condensation

Keywords: number͑s͒: 05.50.ϩq, 71.10.Ϫw, 75.10.Jm

After reviewing the concept of vison excitations in Z 2 dimer liquids, we study the liquid-crystal transition of the quantum dimer model on the triangular lattice by means of a semiclassical spin-wave approximation to the dispersion of visons in the context of a "soft-dimer" version of the model. This approach captures some important qualitative features of the transition: continuous nature of the transition, linear dispersion at the critical point, and ͱ 12ϫ ͱ 12 symmetry-breaking pattern. In a second part, we present a variational calculation of the vison dispersion relation at the Rokhsar-Kivelson ͑RK͒ point, which reproduces the qualitative shape of the dispersion relation and the order of magnitude of the gap. This approach provides a simple but reliable approximation of the vison wave functions at the RK point.

I. INTRODUCTION

Since they have been shown to possess resonating valence bond ͑RVB͒ phases on the triangular, 1 kagome, 2 and other ͑nonbipartite͒ lattices, 3 quantum dimer models ͑QDMs͒ have been one of the main paradigms in the field of quantum spin liquids. These models, where the Hilbert space is spanned by hard-core dimer coverings of the lattice, are expected to capture the phenomenology of quantum antiferromagnets where the wave function is dominated by short-range valence bond configurations. On the triangular lattice, the simplest QDM is defined by the Hamiltonian:

H = -t r + H.c. +V r + ͑1͒
where the sum runs over all plaquettes ͑rhombi͒ including the three possible orientations. The kinetic term, controlled by t, flips the two dimers on every flippable plaquette, i.e., on every plaquette with two parallel dimers, while the potential term controlled by the interaction V describes a repulsion ͑V Ͼ 0͒ or an attraction ͑V Ͻ 0͒ between nearest-neighbor dimers.

The RVB phase is now relatively well understood. The first result goes back to Rokhsar and Kivelson, 4 who showed that, for V / t =1 ͑the Rokhsar-Kivelson or RK point͒, the ground state is the sum of all configurations with equal amplitudes:

͉RK͘ = 1 ͱ N ͚ c ͉c͘. ͑2͒
Since then, dimer-dimer correlations have been shown to be short ranged in a range of parameters below V / t =1, 1,5 and the excitation spectrum to be gapped. 1,[6][7][8] Besides, it has topologically degenerate ground states on non-simplyconnected clusters. 1,5,7 This degeneracy is not related to a standard symmetry breaking. Indeed, there is no local order parameter, 5,9 but only Z 2 topological order. 10 In QDM, the nature of the phase transition from a liquid to a solid is a long standing problem. It goes back to Jalabert and Sachdev 11 ͑see also Ref. 12͒, who studied a threedimensional frustrated Ising model related to the square lattice QDM.

In a more general context, Senthil and Fisher 13,14 showed that models of Mott insulators can be cast in the form of a Z 2 gauge theory. In this language, the transitions from a fractionalized insulator to a conventionally ordered insulator appears to be a condensation of Z 2 vortices ͑dubbed visons͒. Using a duality relation, 15 they showed that such transitions correspond to an ordering transition in a frustrated Ising model in transverse field. As we will see, this applies to the present QDM.

Building on a mapping between QDM's at V = 0 and Ising models in a transverse field, Moessner et al. 16,17 have developed a Landau-Ginzburg approach and suggested that the transition could be continuous and in the three-dimensional O͑4͒ universality class.

Numerical evidence in favor of this scenario has been obtained with Green's function quantum Monte Carlo by Ralko et al., who have shown that, at the transition point, the static form factor of the crystal decreases to zero on the crystal side of the transition, while the dimer gap also decreases to zero on the liquid side. 18 More recently, the vison spectrum has also been numerically determined using Green's function quantum Monte Carlo, with the conclusion that a soft mode, indeed, develops at the transition. 19 In spite of these results, a simple picture for the wave functions of the visons and the evolution of their spectrum is still missing. One difficulty is that, like any vortex, these excitations cannot be created by operators which are local ͑in the dimer variables͒. To attack this problem, we follow two strategies. First of all, we look at the problem in the context of a Z 2 gauge theory on the triangular lattice. As usual, this theory can be mapped onto a dual Ising model. 15 The duality transforms the nonlocal excitations of the gauge theory into local excitations, which we study using a semiclassical approach. This simple 1 / S expansion already captures most aspects of the confinement-deconfinement transition ͑soft mode and condensation͒ of the Z 2 gauge theory. Second, building on the explicit form of the excitations of the Z 2 gauge theory, we construct single vison wave functions for the QDM and study the properties of the system at the RK point in the variational Hilbert space spanned by these states. These "variational visons"-living on the triangular plaquettes and experiencing a flux emanating from each site of the lattice-turn out to be linearly independent. The associated dispersion is in good qualitative agreement with that obtained in Monte Carlo simulations, and the value of the gap at the RK point ͑⌬ var = 0.119͒, obtained with a single variational parameter, has the correct order of magnitude ͑Monte Carlo simulations 6 give 0.089͒. Improving quantitatively further these variational results would require more adjustable parameters to account in more details for the local ͑dimer-dimer, etc.͒ correlations in the vicinity of the core of the vortex, a task which has not been carried out here.

II. VISONS IN Z 2 GAUGE THEORY

The connection between QDMs and Z 2 gauge theories has already been discussed from different perspectives ͑see, in particular, Refs. 2, 3, 12, and 20͒, and is rooted in the existence, in both families of models, of Ising-like degrees of freedom subjected to local constraints ͑hard-core constraints for the dimers, Gauss law for the gauge theory͒.

In this section, we review two known mappings: ͑i͒ from a Z 2 lattice gauge theory on the hexagonal lattice to the triangular lattice QDM ͑valid in a particular limit, Sec. II A͒ and ͑ii͒ from the Z 2 theory to its dual frustrated Ising model ͑Sec. II B͒. The latter Ising model is then studied using a semiclassical approximation in Secs. II B 2 and II B 3. Since the Z 2 gauge theory-QDM mapping is formally justified only in the confined phase of the gauge theory, the relevance of the results obtained in the RVB ͑deconfined͒ phase of the QDM is not guaranteed a priori. Some reasons why this is expected to be the case are discussed in the next section, when we build on the results obtained for the excitations of the Z 2 gauge theory in its deconfined phase to construct elementary excitations in the RVB phase of the QDM.

A. Z 2 gauge theory

To write down a Z 2 gauge theory analogous to a QDM, the starting point is to define Pauli matrices ជ l on the bonds l of the triangular lattice such that l x = -1 if bond l is occupied by a dimer, and +1 if it is empty, while l z changes the state of bond l. Since each kinetic term of the QDM flips the dimers around a rhombus, it corresponds to a product of four l z around this rhombus. This would, however, not keep the structure of the simple Z 2 gauge theories, in which the kinetic term acts on an elementary plaquette of the lattice, a crucial ingredient to get a simple Ising model by duality. 15 An alternative is to consider the more standard Z 2 gauge theory which, in its Hamiltonian formulation, is defined by

H = H J + H ⌫ = -J ͚ l l x -⌫ ͚ i ͟ l͑i͒ l͑i͒ z , ͑3͒
where i runs over the sites of the dual honeycomb lattice, and l͑i͒ are the three bonds forming the triangular plaquette around site i ͑see Fig. 1͒. The hallmark of this model is to have local conserved quantities. Indeed,

ͫ H, ͟ l͓a͔ l͓a͔ x ͬ = 0,

͑4͒

where a is a site of the triangular lattice, and the product over l͓a͔ runs over the six links emanating from a ͑see Fig. 1͒. This allows one to define different sectors according to whether ͟ l͓a͔ ͓a͔ x is equal to +1 or -1. Since in the QDM the number of dimers emanating from a given site is exactly equal to 1, it is clearly better to consider the sector where

͟ l͓a͔ l͓a͔ x = -1 ͑5͒
for all a since this forces the number of dimers emanating from a site to be odd ͑defining an odd Ising gauge theory in the terminology of Ref. 20͒. Then, the true constraint is recovered in the limit ⌫ / J → 0 if J Ͼ 0 since, in the ground state, the number of dimers is then minimal. A bona fide QDM is then recovered if H ⌫ is treated with degenerate perturbation theory. Since H ⌫ changes the number of dimers, its effect vanishes to first order. To second order, however, one recovers exactly the QDM of Eq. ͑1͒ with V = 0 and t = ⌫ 2 / J. So, if J / ⌫ӷ1, the Z 2 gauge theory maps onto the QDM at V / t = 0. At finite J / ⌫, the model can be viewed as a "soft-dimer" model, where one, three, or five dimers may touch a given site.

As we shall see below, the limit J / ⌫ӷ1 lies deep inside the confined phase of the gauge theory, and from previous work on the QDM, it is known that at V / t = 0, the model is in a valence bond crystal phase. It would, of course, be very interesting to connect the two models away from this limit, when the gauge theory is in its deconfined phase and the QDM is in the dimer liquid phase. An interesting step in this a 1)

2) l[a] 3) l<i l(i) i i FIG.
1. ͑Color online͒ Some useful definitions of sets of bonds: ͑1͒ l͑i͒ is the set of the three bonds forming the edges of the plaquette i. ͑2͒ l͓a͔ is the set of the six ͑fat͒ bonds emanating from the triangular site a. ͑3͒ l Ͻ i denotes the bonds forming a "zigzag" string extending ͑to the left͒ from the triangular plaquette i to an edge of the lattice.

direction is provided by higher order perturbation theory. Indeed, to fourth order in ⌫ / J, a repulsion between dimers V = ⌫ 4 / 2J 3 is generated. However, other terms of the same order, involving dimer shifts along loops of length 6, are also generated ͑see Appendix A͒. So, a rigorous mapping does not extend beyond the V / t = 0 point.

B. Dual Ising model

As usual, this Z 2 gauge theory is best analyzed by mapping it onto a dual Ising model in a transverse field.

Model

This can be achieved by introducing spin- 1 2 operators ជ i on the dual honeycomb lattice:

i x = ͟ l͑i͒ l͑i͒ z , i z = ͟ lϽi l x , ͑6͒
where l Ͻ i represents all bonds cutting a straight path ͑say, horizontal, see Fig. 2 or 1͒ starting at i ͑we implicitly assume a finite lattice with open boundary conditions͒. Combined with the constraint, this definition implies that

i z j z = M ij l x ͑7͒
for two neighboring sites i , j separated by the bond l. M ij = Ϯ 1 is such that each hexagon has exactly one M ij =-1 bond ͑see Fig. 3͒. In terms of these spin operators, the Hamiltonian is the fully frustrated Ising model 21 ͑FFIM͒ discussed by Moessner and Sondhi: 17,26

H = H J + H ⌫ = -J ͚ ͗i,j͘ M ij i z j z -⌫ ͚ i i x . ͑8͒
Note that choosing other paths to define i z leads to other signs for M ij , which, however, are always such that an odd number of minus signs appear around each hexagon, leading to the same Ising model up to a gauge transformation. The present choice leads to the smallest unit cell ͑four honeycomb sites͒.

In the following, we study the phase diagram and the excitations of this model within a semiclassical ͑large S͒ approximation.

Classical phase diagram

First, we determine the classical ground state of the model as a function of ⌫ / J. To this end, we replace the spin-1 2 operators by classical three-component vectors of unit length. Since the y component does not appear in the Hamiltonian, it is clear that in the ground state, the magnetization must lie in the x-z plane. To find the lowest energy solution, we use the following parametrization:

i z → ͑i͒, ͑9͒ i x → ͱ 1 -͑i͒ 2 ,

͑10͒

with ͉͑i͉͒ ഛ 1. The corresponding energy is

E = - J 2 ͚ ij M ij ͑i͒͑j͒ -⌫ ͚ i ͱ 1 -͑i͒ 2 .

͑11͒

Initially defined for nearest neighbors ͓Eq. ͑7͔͒, the coefficients M ij have been upgraded to a matrix M by setting all other elements to 0. This matrix describes the motion of a particle on a honeycomb lattice with four sites per unit cell and a flux per hexagonal plaquette. It reduces to a 4 ϫ 4 matrix after Fourier transformation. The eigenvalues associated with the momentum

k = ͑k x , k y ͒ are 17 m k 1,2,3,4 = Ϯ ͱ 3 Ϯ ͱ 2͓3 + cos͑2k x ͒ -cos͑k x + k y ͒ + cos͑k x -k y ͔͒.
Likewise, we denote by the column vector of components ͕ i ͖.

⌫ / J Ͼ ͱ 6. For ⌫ӷJ, we expect the ͑i͒'s to be small ͑or zero͒. We can, therefore, expand Eq. ͑11͒ to quadratic order:

E = - J 2 t M + 1 2 ⌫ t . ͑12͒
The largest eigenvalue of M being ͱ 6, we find that, as long as ⌫ / J Ͼ ͱ 6, the energy is minimized by = 0 and all spins point in the x direction.

0 Ͻ⌫/ J ഛ ͱ 6. At ⌫ / J = ͱ 6, all the real eigenvectors of M for the eigenvalue ͱ 6 ͑satisfying ͉ i ͉ ഛ 1͒ minimize Eq. ͑12͒. Let us choose four complex vectors which form a basis of the subspace associated with the eigenvalue ͱ 6: ͑Ref. 17͒

v 1 ͑x,y͒ = v 3 * = ΄ e 5i/12 /F e -i/6 /F 1 e -i/12 ΅ expͩ i 6 x + i 2 yͪ , v 2 ͑x,y͒ = v 4 * = ΄ e i/12 e -5i/6 1/F e -5i/12 /F ΅ expͩi 5 6 x + i 2 yͪ , F = 2 sin͑5/12͒, ͑13͒
where the four entries of the vector refer to the four sites in the unit cell ͑numbered as in Fig. 3͒ andx , y are the ͑Bravais͒ coordinates of the unit cell. These eigenvectors correspond to the four points labeled B in the rectangular Brillouin zone of Fig. 6.

The most general real eigenvector can be parametrized by three angles ␣ 1 , ␣ 2 , and ␤ and a normalization factor :

͑i͒ = ͓cos͑␤͒R"v 1 ͑i͒e i␣ 1 … + sin͑␤͒R"v 2 ͑i͒e i␣ 2 …͔. ͑14͒
To find the ground state, the energy has to be minimized as a function of these three parameters. The analysis is made easier when one realizes that the second and fourth moments of the spin deviations are independent of the three angles. Indeed,

1 4 4 , ͑16͒

N ͚ i ͑i͒ 2 = r 2 2 , ͑15͒ 1
N ͚ i ͑i͒ 4 = 2r
with

r 2 = 1 4 ͩ1+ 1 F 2 ͪ. ͑17͒
Replacing by Eq. ͑14͒ in the expression for the energy ͓Eq. ͑11͔͒, we get up to order 4 :

E N = -⌫ + 1 2 ͑⌫ -J ͱ 6͒r 2 2 + 1 4 ⌫r 4 4 + O͑ 6 ͒. ͑18͒
Minimizing with respect to 2 gives

2 = J ͱ 6 -⌫ f 2 ⌫ . ͑19͒
For ⌫ϾJ ͱ 6, the energy is minimized for = 0. For ⌫ഛJ ͱ 6, we have to expand Eq. ͑11͒ to the sixth order to find the angles ␣ 1 , ␣ 2 , and ␤ which minimize the classical energy. Indeed, Eq. ͑18͒ shows that, up to order 4 , the energy is independent of the three angles. At this order, the energy is constant and minimum on a three-dimensional sphere, a con-sequence of the O͑4͒ symmetry discovered by Moessner and Sondhi. 17 One, therefore, has to go to the next order to find the actual minima. At sixth order in , the energy is minimized when 1

N ͚ i ͑i͒ 6 ͑20͒
is minimum. A numerical investigation shows that the solutions ͑for ⌫ / J close to but below ͱ 6͒ are 48-fold degenerate and can be deduced from each other by symmetry operations ͑12 translations and 4 point-group operations͒. This confirms the result obtained previously on symmetry grounds. 16,17 Motivated by the FFIM-dimer model correspondence, we are interested in the average "dimer density:"

d ij = 1 2 ͑1 -M ij ͗ i z j z ͒͘ ͑21͒
for all pairs ͑i , j͒ of nearest-neighbor honeycomb sites. In the classical limit, this may be approximated by

d ij = 1 2 ͓1 -M ij ͑i͒͑j͔͒, ͑22͒
where ͓Eq. ͑14͔͒ is the classical ground state. Close to the transition at ⌫ c = J ͱ 6, scales as ϳ ͱ ⌫ c -⌫ ͓Eq. ͑19͔͒ and the dimer density, thus, shows small deviations about 1 2 . To visualize the "dimerization" pattern in the vicinity of ⌫ c , the appropriate quantity is, therefore, a relative rescaled "dimer" density defined by D ij = 1 2 ͑d ij -1 2 ͒, and plotted in Fig. 4 for one of the 48 ground states. The obtained pattern is highly reminiscent of the ͱ 12ϫ ͱ 12 valence bond crystal ͑VBC͒ observed in the triangular lattice QDM. In particular, the 48 spin configurations give only 12 different dimer patterns, with a smaller unit cell containing 12 sites of the triangular lattice. It is also interesting to notice that ͑i͒ vanishes in the triangles located inside the large diamonds ͑marked with a dot in Fig. 4͒. From Eq. ͑10͒, the corresponding sites have a magnetization pointing exactly in the x direction, which means a total absence of vison and local correlations identical to that of the J = 0 paramagnetic ground state. The dimer density d ij is uniform inside these diamonds, in qualitative agreement with the intriguing observation made in Ref. 18.

Semiclassical analysis for ⌫ Õ J Ð ͱ 6

If J = 0, the ground state is the fully polarized state i x = 1. The first ͑degenerate͒ excited state is obtained by flipping one spin at some arbitrary triangle. To first order in J / ⌫, the spectrum is obtained by diagonalizing the perturbation H J in this subspace. This amounts to diagonalizing M, and leads to the eigenvalues ⑀ k i =2⌫ + Jm k i , a result already obtained previously. 17 To go further, beyond the limit J / ⌫Ӷ1, we perform a 1 / S semiclassical expansion. We generalize the spin-1 / 2 Hamiltonian to an arbitrary value S of the spin:

H = - J S 2 ͚ ͗i,j͘ M ij S i z S j z - ⌫ S ͚ i S i x ͑23͒
͓which reduces to Eq. ͑8͒ when S =1/ 2͔. Spin deviations away from the x directions are represented using Holstein-Primakoff bosons. To leading order in 1/S:

S i z = 1 2 ͱ 2S͑b i † + b i ͒, ͑24͒ S i x = S -b i † b i . ͑25͒
One truncates the Hamiltonian to quadratic order in b i , and diagonalizes it through a Bogoliubov transformation:

b i † = ͚ j U ij a j † + ͚ j V ij a j . ͑26͒
For the a j † operators to be bosonic creation operators, the matrices U and V must satisfy

U † U -V † V = 1. It is convenient to introduce a unitary matrix ⍀ which transforms M into a diagonal matrix M ˜: M = ⍀M ˜⍀ † , ͑27͒ M ˜kp = m k ␦ kp .

͑28͒

In this new basis, we can look for diagonal solutions for U and V:

U = ⍀U ˜⍀ † , V = ⍀V ˜⍀ † , ͑29͒ U ˜kp = u k ␦ kp , V ˜kp = v k ␦ kp .

͑30͒

Since M is real, we can choose ⍀ O͑N͒. After some algebra, one finds that the u k and v k which diagonalize H are given by

u k = cosh͑ k ͒, ͑31͒ v k = sinh͑ k ͒, ͑32͒ tanh͑2 k ͒ = - m k m k + 2⌫/J . ͑33͒
The energies of the Bogoliubov excitations are given by

⑀ k = J S ͱ ⌫ J ͩm k + ⌫ J ͪ.

͑34͒

As expected, the spectrum becomes gapless at ⌫ / J = ͱ 6 and we recover the localized spin flip energy ⑀ k Ӎ ⌫ / S when ⌫ → ϱ. The Bogoliubov spectrum is shown Fig. 5 for a few values of ⌫ at or above ⌫ c . One clearly sees that the dispersion is linear around the gapless points at the transition. This behavior is remarkably similar to that found numerically by Ralko et al. 19 for the QDM by Green's function quantum Monte Carlo. In the next section, we build on this resemblance to develop a variational approach to the vison spectrum of the QDM. Note that, as stated above, the Z 2 gauge theory and the QDM can only be rigorously mapped onto each other deep into the confined ͑VBC͒ phase, and the resemblance between their spectra in the deconfined ͑RVB͒ phases might seem fortuitous. A somewhat deeper connection will be described in the next section. In the Ising model, elementary excitations for large enough ⌫ / J are spin flips, induced by i z , which delocalize and get dressed under the effect of H J . In the equivalent dual gauge theory, this excitation is produced by the nonlocal string operator i z = ͟ lϽi l x . By analogy, it is natural to define a "point" vison creation operator for the QDM by 5,22

B B B B 0 0 0 0 4 4 4 4 A C A C A C A C 3 3 3 3 ε(κ) ε(κ) ε(κ) ε(κ) 1 1 1 1 A A A A 2
V i = ͑-1͒ N ˆ͑i͒ , N ˆ͑i͒ = ͚ lϽi n ˆl, ͑35͒
where the dimer operator n ˆl is defined by n ˆl = 1 if bond l is occupied, and 0 otherwise. Again, we consider here a finite lattice with open boundary conditions. Let us first see to which extent this operator is the analog of the vortex creation operator in Ising gauge theories.

A. Z 2 gauge structure of quantum dimer models

In gauge theories, the Wilson loop operator defined by

W ץ⍀ = ͟ lץ⍀ l z ͑36͒
plays a central role since it allows one to distinguish ͑in the absence of matter field, as here͒ the deconfined and confined phases depending on whether its ground state expectation value ͑or flux͒ tends to zero exponentially with the perimeter of the domain ⍀ ͑deconfined͒ or with the area of the domain ͑confined͒. For a Z 2 gauge theory, W ץ⍀ 2 = 1 and the flux going through ⍀ measured by W ץ⍀ can only take two values Ϯ1. In that respect, an interesting property of the operator i z = ͟ lϽi l

x is that it changes the flux between -1 and +1 if the site i is inside the domain. Then, i z j z commutes with the Wilson loop and does not change the flux unless i and j sit on opposite sides of the boundary. Since the deconfined phase can be accessed from the J Ӷ⌫ limit, the ground state can be thought of as an Ising paramagnet ͑ x = 1 and W ץ⍀ = 1 every-where͒ "dressed" perturbatively by successive applications of J i z j z , where i and j are nearest neighbors. Since J i z j z only changes the flux for pairs ij across the boundary, the expectation value of the Wilson loop operator behaves according to a perimeter law. A similar perturbative argument can be used to derive the area law in the confined phase.

Most of these standard Z 2 gauge theory results apply to the QDM, with one difficulty however. The Wilson loop operator cannot be defined in the same way since flipping the dimer occupation along a loop will often lead to an unphysical state that violates the condition of having exactly one dimer emanating from each site. Different ways to overcome this problem can be envisaged. One possibility is to accept that the Wilson loop operator gives zero when applied to states that are nonflippable along the loop. The property W ץ⍀ 2 = 1 is lost, but W ץ⍀ has eigenvalues 0, +1, or -1, and eigenstates with eigenvalues +1 or -1 are still interchanged when a vison operator ͓Eq. ͑35͔͒ is applied inside the domain. So confinement or deconfinement of visons is still expected to lead to area or perimeter laws.

Alternatively, starting from the observation that, when it does not give zero, the Wilson loop operator just shifts the dimers along the contour, one can try to define a ͑compli-cated͒ flux operator in the dimer space which shifts the dimer along a "fattened" contour which "adapts" locally to the dimer configuration it acts upon. Although tricky to define in practice, 27 this operator would have the advantage that its square is still equal to 1, preserving the manifest Z 2 structure of the theory.

Another way to underline the deep connection between the two models is to consider the Ising model as a soft-dimer model and introduce the projection operator onto the hardcore dimer Hilbert space, defined by

P ˆ= ͟ p ͑n ˆp -3͒͑n ˆp -5͒ ͑1 -3͒͑1 -5͒ , ͑37͒
where the operator n ˆp counts the number of frustrated bonds ͑dimers͒ around the plaquette p:

n ˆp = 1 2 ͚ i=1 6 ͑1 -M i,i+1 i z i+1 z ͒. ͑38͒
By construction, P ˆ= 1 on all the Ising configurations, which correspond to a valid hard-core dimer covering ͑of the triangular lattice͒ and P ˆ= 0 otherwise. Now, P ˆcommutes with ͟ i i x since all terms in P ˆcontain an even number of z operators. This means that P ˆconserves the total flux, so that a spin state with -1 flux becomes a dimer wave-function with an odd number of visons after projection.

B. Point vison

Equation ͑35͒ defines the simplest operator which changes the flux, that is, which creates a vison. So we may consider

͉i͘ = V͑i͉͒RK͘ = 1 ͱ N ͚ c ͑-1͒ N ˆ͑i͒ ͉c͘ ͑ 39a͒ = 1 ͱ N P ˆi z ͉͕ j x = + 1͖͘ ͑39b͒
as a first variational approximation to the true lowest eigenstate of H in the -1 flux sector ͑N is the total number of dimer coverings͒. The ground state ͉RK͘ is a zero-energy eigenstate of H, implying that the expectation value E K Ͻ 0 of the kinetic energy term exactly compensates the expectation value E P Ͼ 0 of the potential energy term. Indeed,

E P = 3Np 0 = -E K , ͑40͒
where 3N is the total number of rhombi ͑N the number of sites͒ and p 0 the probability to have two parallel dimers on a given rhombus in the classical dimer problem ͑with uniform measure over all dimer configurations͒. Using the Pfaffians method, 23,[START_REF] Krauth | Statistical Mechanics: Algorithms and Computations ͑Oxford[END_REF] one finds in the thermodynamic limit p 0 Ӎ 0.093 331 0104 . . . . ͑41͒

In ͉i͘ the expectation value of the potential energy term is the same as in ͉RK͘. In fact, all the observables which are diagonal in the dimer basis commute with z and, thus, have the same expectation value in the ground state ͉RK͘ and in the trial wave function ͉i͘. The increase of the energy is only kinetic. Because the kinetic terms corresponding to the three diamonds around i anticommute with i z , their expectation values have changed sign. Thus, we find ͗i͉H͉i͘ -͗RK͉H͉RK͘ = 6p 0 Ӎ 0.56. ͑42͒

C. Dispersing point vison

One can lower the energy by constructing plane waves. To compute the associated dispersion relation we have to evaluate the following matrix elements:

S ij 0 = ͗i͉j͘, ͑43͒
H ij 0 = ͗i͉H͉j͘. ͑44͒
We begin with the overlap matrix

S ij 0 = 1 N ͚ c ͗c͉͑-1͒ N ˆ͑i͒+N ˆ͑j͒ ͉c͘ ͑ 45͒
͓N ˆ͑i͒ and N ˆ͑j͒ are defined in Eq. ͑35͔͒ which simplifies to

S ij 0 = ⑀ ij 1 N ͚ c ͗c͉͑-1͒ N ˆ͑i,j͒ ͉c͘, ͑46͒ ⑀ ij = ͑-1͒ N 0 ͑i,j͒ , ͑47͒
where the local operator N ˆ͑i , j͒ counts the number of dimers across some path ⍀ ij connecting the triangles i and j ͑see Fig. 2͒, and N 0 ͑i , j͒ is equal to that number in the reference configuration c 0 , chosen with all the dimers horizontal ͑Fig. 2͒. This follows from two simple properties: ͗c ͉͑-1͒ N ˆ͑i͒+N ˆ͑j͒+N ˆ͑i,j͒ ͉ c͘ is independent of the configuration c, and ͗c 0 ͉͑-1͒ N ˆ͑i͒+N ˆ͑j͒ ͉ c 0 ͘ = 1. We finally write

S ij 0 = ⑀ ij ͗͑-1͒ N ˆ͑i,j͒ ͘, ͑48͒
where ͗¯͘ represents the average with equal weight over all dimer coverings. The average of any such diagonal observable can be computed using the Pfaffian of a ͑modified͒ Kasteleyn matrix. 23,[START_REF] Krauth | Statistical Mechanics: Algorithms and Computations ͑Oxford[END_REF] In the present case, the "string" observable ͑-1͒ N ˆ͑i,j͒ is coded in the Kasteleyn matrix by changing the signs of the matrix elements corresponding to bonds crossed by ⍀ ij , i.e., setting some bond fugacities to -1. To evaluate numerically such an expectation value, we construct the modified Kasteleyn matrices corresponding to a large enough triangular lattice, in which the path ⍀ ij is embedded. Finitesize effects decay exponentially with the system size, so that lattices with 28ϫ 28 sites ͑with periodic boundary condi-tions͒ can safely be used to evaluate S ij 0 with high accuracy up to distances d Ӎ 10 between triangles i and j.

The matrix elements S ij 0 , plotted in Fig. 7 as a function of the distance between i and j, decay exponentially ͑a result anticipated by Read and Chakraborty 22 ͒. These quantities have already been evaluated by Ioselevich et al. 5 using a classical Monte Carlo sampling. By construction, a product like ⑀ i 1 ,i 2 ⑀ i 2 ,i 3 ¯⑀i n ,i 1 ͑where i 1 , i 2 , ... ,i n form a closed loop of triangles͒ is equal to the parity of the number of sites enclosed in the loop. Thus, the signs of the matrix elements S ij 0 are similar to those of the hopping amplitude of a particle moving on the hexagonal lattice and subjected to a magnetic field corresponding to half a flux quantum per hexagon. 28 In such a case, the magnetic unit cell has to be doubled compared to the original lattice cell. Since the original unit cell contains one triangular site and two triangular plaquettes, the magnetic one contains four triangular plaquettes and, thus, four sites of the hexagonal lattice. In Fourier space, S 0 ͑k͒ is a 4ϫ 4 matrix, as the matrix M discussed previously. The same is also true for the Hamiltonian matrix elements described below.

The spectrum of the overlap matrix S 0 ͑k͒ is plotted in Fig. 8 ͑for k describing a representative path in the Brillouin zone͒. As an important result, the eigenvalues are strictly positive for all k.

Let us now turn to the matrix elements of the Hamiltonian. We wish to transform Eq. ͑44͒ into an expression which can be evaluated using the Pfaffians, that is, the expectation value of a diagonal observable in the dimer basis. First, we write H as a sum of projectors

H = 2 ͚ r 0 ⌸ ˆr0 , ͑49͒ FIG. 7.
͑Color online͒ Overlap ͉͗i ͉ j͉͘ between two point-vison states ͓defined in Eq. ͑39a͔͒ as a function the distance d ij between i and j. The ͑blue͒ squares correspond to Eq. ͑43͒ and the three ͑red͒ crosses correspond to Eq. ͑59͒ with ␣ = -0.8. The calculations are done using an exact evaluation of the Pfaffians on a 28-site lattice with periodic boundary conditions. Dashed line: guide for the eye corresponding to an exponential decay with the dimer-dimer correlation length -1 = 0.76 ͑Ref. 5͒.

⌸ ˆr0 = ͉ r 0 ͗͘ r 0 ͉, ͑50͒ |ψ r 0 = 1 √ 2 - ͑51͒
and expand Eq. ͑44͒ into

H ij 0 = 2 N ͚ c 1 ,c 2 ͚ r 0 ͑-1͒ N͑c 1 ,i͒+N͑c 2 ,j͒ ͗c 1 ͉⌸ ˆr0 ͉c 2 ͘, ͑52͒
with the notation N͑c , i͒ = ͗c ͉ N ˆ͑i͉͒c͘ = Ϯ 1. Because of the double sum ͚ c 1 ,c 2 , this does not yet have the form of a diagonal observable amenable to an evaluation with Pfaffians. To go further, we note that ͗c 1 ͉ ⌸ ˆr0 ͉ c 2 ͘ vanishes if c 1 and/or c 2 is not flippable around the rhombus r 0 . In addition, ͗c 1 ͉ ⌸ ˆr0 ͉ c 2 ͘ =0 if c 1 and c 2 differ anywhere outside r 0 . So we can restrict the sum to pairs of configurations ͑c , c ¯͒, which differ by a dimer flip in r 0 and which are identical elsewhere on the lattice:

H ij 0 = 2 N ͚ r 0 ͚ ͑c,c ¯͒ F r 0 ͑c͒=1 ͓͑-1͒ N͑c,i͒ ͗c͉ + ͑-1͒ N͑c ¯,i͒ ͗c ¯͉͔⌸ ˆr0 ͓͑-1͒ N͑c,j͒ ϫ͉c͘ + ͑-1͒ N͑c ¯,j͒ ͉c ¯͔͘, ͑53͒
where F r 0 ͑c͒ =1 if c is flippable on rhombus r 0 , and F r 0 ͑c͒ = 0 otherwise. c and c ¯only differ inside r 0 , so the signs ͑-1͒ N͑c,i͒ and ͑-1͒ N͑c ¯,i͒ are the same if i is not inside r 0 , and are opposite if i r 0 . Let us note ␦ i,r 0 =-1 if i r 0 , and ␦ i,r 0 = 1 otherwise.

This leads to

H ij 0 = 2 N ͚ r 0 ͚ ͑c,c ¯͒ F r 0 ͑c͒=1 ͑-1͒ N͑c,i͒+N͑c,j͒ ͓͗c͉ + ␦ i,r 0 ͗c ¯͉͔⌸ ˆr0 ͓͉c͘ + ␦ j,r 0 ͉c ¯͔͘.

͑54͒

Using the explicit form of ⌸ ˆr0 , we get

͓͗c͉ + ␦ i,r 0 ͗c ¯͉͔⌸ ˆr0 ͓͉c͘ + ␦ j,r 0 ͉c ¯͔͘ = 1 2 ͑1 -␦ i,r 0 ͒͑1 -␦ j,r 0 ͒ ͑55͒
and, finally,

H ij 0 = ͚ r 0 ͑1 -␦ i,r 0 ͒͑1 -␦ j,r 0 ͒ 1 2N ͚ c F r 0 ͑c͒=1 ͑-1͒ N͑c,i͒+N͑c,j͒ ͑56a͒ =⑀ ij͚ r 0 ͑1 -␦ i,r 0 ͒͑1 -␦ j,r 0 ͒ 1 2N ͚ c F r 0 ͑c͒=1 ͗c͉͑-1͒ N ˆ͑i,j͒ ͉c͘.
͑56b͒

The last expression is the average of a diagonal observable and can, thus, be evaluated using Pfaffians:

H ij 0 = 1 2 ⑀ ij͚ r 0 ͑1 -␦ i,r 0 ͒͑1 -␦ j,r 0 ͒͗͑-1͒ N ˆ͑i,j͒ F ˆr0 ͘, ͑57͒
where we used an operator notation for the "flippability" F ˆr0 ͉c͘ = F r 0 ͑c͉͒c͘. To get a nonzero contribution, there must be at least one rhombus r 0 containing i and j. So H ij 0 =0 if i and j are not nearest neighbors. Using modified Kasteleyn matrices in a way similar to that described for the overlap matrix S 0 , the nonzero matrix elements can be calculated. We find ͉H ij 0 ͉ =6p 0 for i = j; ͉H ij 0 ͉ =2p 0 when i and j are first neighbors.

Although the matrix elements of H are simple, the band structure is, however, not that of a simple tight-binding Hamiltonian, because of the nonorthogonality of the present variational vison states. In Fourier space, S 0 ͑k͒ and H 0 ͑k͒ are 4 ϫ 4 matrices. The spectrum is obtained by solving the generalized eigenvalue problem H 0 ͑k͒ k -E͑k͒S 0 ͑k͒ k = 0. The results are shown in Fig. 9 ͑␣ = 0 curve͒. The minimum of ⑀ k is found at ͑k x , k y ͒ = ͑ / 6, / 2͒, in agreement with the Monte FIG. 8. ͑Color online͒ Spectrum of the overlap matrix S ␣ ͑k͒ ͓Eqs. ͑43͒ and ͑59͔͒, for different values of ␣ ͑0, -0.2, -0.5, and -0.8͒, along the path A → B → C → A in the Brillouin zone ͑see Fig. 6͒. The bottom panel is a zoom on the lowest eigenvalues of S ␣ . In these calculations, the matrix elements S ij ␣ are neglected for triangles i and j at distance d ജ d max =10 ͑102th neighbor on the hexagonal lattice͒. Up to this distance, the finite-size lattice used in the calculation ͑28ϫ 28 sites͒ gives practically the infinite volume limit for S ij ␣ . The value of d max used here is large enough to ensure a good convergence of the spectrum since the curves obtained with a smaller truncation distance ͑d max Ӎ 8.47-shown here with the same colors͒ are almost superposed with that for d max = 10, except for the bottom of the spectrum at ␣ = -0.5 and ␣ = -0.8. The overlap spectrum turns out to be gapped for ␣ = 0 and ␣ = -0.2 ͑and probably at ␣ = -0.5 and ␣ = -0.8 too͒, indicating the linear independence of the vison states.

Carlo results. 6,18 The corresponding energy is ⑀ min = 0.16, which is significantly lower than the energy ͑0.56͒ found above for a completely localized point vison ͉i͘. The bottom of this variational point-vison band is, however, still high compared to Ivanov's 6 estimate ͑0.089͒ of the gap in the vison sector. The next section provides an improved family of variational states.

D. Dressed vison

The vison wave function of Eq. ͑39a͒ differs from the RK wave functions only through minus signs. It has the "accidental" property that dimer-dimer correlations are the same as in the ground state. 29 As a consequence, the excitation energy of such a state is carried only by the kinetic part of the Hamiltonian. Clearly, some local reweighting of the dimer configurations in the vicinity of the vortex core would allow an optimized balance between the potential and kinetic costs of the excitation, and would lead to an improved variational wave function. This amounts to "dressing" locally the initial point-vison state by some even-but fluctuatingnumber of additional point visons.

As a simple improvement of the vison wave function, we introduce a variational parameter ␣ to reweight the configurations depending on their flippability at the core of the vison. This gives the following vison state:

͉i,␣͘ = 1 ͱ N ͚ c ͑-1͒ N ˆ͑i͒ ͑1 + ␣F ˆi͉͒c͘, F ˆi = F ˆr1 ͑i͒ + F ˆr2 ͑i͒ + F ˆr3 ͑i͒ , ͑58͒
where F ˆi͉c͘ = ͉c͘ if the dimer configuration c is flippable around one of the three rhombi r 1 ͑i͒, r 2 ͑i͒, and r 3 ͑i͒ containing the triangle i, and F ˆi͉c͘ = 0 otherwise. As for the point vison of Eq. ͑39a͒, the vison states of Eq. ͑58͒ are not orthogonal and we have to evaluate their overlaps:

S ij ␣ = ͗i,␣͉j,␣͘. ͑59͒
Repeating the transformation leading to Eq. ͑48͒, we get

S ij ␣ = ⑀ ij ͗͑-1͒ N ˆ͑i,j͒ ͑1 + ␣F ˆi͒͑1 + ␣F ˆj͒͘, ͑60͒
which is an expectation value for a diagonal operator that we evaluate using Pfaffians. In addition to the sign changes due to ͑-1͒ N ˆ͑i,j͒ , some entries of the Kasteleyn matrix have to be modified to incorporate the flippability operators F ˆr. More precisely, counting only the coverings which satisfy F ˆr =1 is done by "isolating" the rhombus r, that is, by switching to zero in the Kasteleyn matrix the 14 bonds which connect the sites of rhombus r to their neighbors outside r.

Beyond some distance between the vison cores i and j ͑fourth neighbor on the hexagonal lattice͒, no rhombus can touch simultaneously both triangles. In that case, it can be shown that S ij ␣ = S ij 0 is independent of ␣. The eigenvalues of the overlap matrix S ␣ ͑k͒ are displayed in Fig. 8 for a few selected values of ␣. Although a full convergence as a function of the truncation distance d max ͑see caption of Fig. 8͒ has not been obtained, we believe that there is a finite gap for all the values of ␣ shown here, and that the dressed vison states are linearly independent.

The evaluation of the Hamiltonian matrix elements

H ij ␣ = ͗i,␣͉H͉j,␣͘ ͑ 61͒
for dressed vison can still be done using Pfaffians, but the algebraic manipulations are slightly more lengthy than the previous ones, and we refer the reader to Appendix B. We find nonzero matrix elements H ij ␣ up to ͑and including͒ the ninth neighbor ͑compared to first neighbor for point vison͒.

The results for the variational dispersion relation are shown Fig. 9. The qualitative shape of the lowest band is almost unchanged compared to the point-vison states ͑␣ =0͒, except for an almost uniform shift which lowers the gap. The minimum of E k is found at ͑k x , k y ͒ = ͑ / 6, / 2͒ = B for a variational parameter ␣ Ӎ -0.8, and the corresponding energy is E min = 0.119, about 33% higher than the exact value.

The fact that simple wave functions like those of Eq. ͑39a͒ and ͑58͒ reproduce qualitatively the shape of the exact dispersion relation is presumably due to the fact that the dimer-dimer correlation length is rather small at the RK point ͑of order of one lattice spacing 1 ͒. As a consequence, the exact vison states only differ from Eq. ͑39a͒ at short distances from the vortex core, and the long-distance part ͑string͒, responsible for the flux per hexagon, is essentially 6͒. The curve for ␣ = 0 corresponds to the variational energies of point-vison states ͓Eq. ͑39a͔͒. The absolute minimum is E min = 0.119, found at ͑k x , k y ͒ = ͑ / 6, / 2͒ = B for a value Ӎ-0.8 of the variational parameter ␣. The exact value of the gap ͑0.089͒ ͑Ref. 6͒ is marked as a dotted horizontal line. In this calculation, the matrix elements S ij ␣ are neglected for triangles i and j at distance d ജ d max = 10. This value is large enough to ensure a perfect convergence of the spectrum below E Շ 0.6, as can be checked from the fact that the dispersion curves obtained with a smaller truncation distance ͑d max Ӎ 8.47-shown here with the same colors͒ are practically identical. Some higher energy states, however, are not fully converged.

exact. This is, of course, no longer true away from the RK point and in the direction of the crystal, where the correlation length rapidly grows.

In an attempt to extend this variational approximation away from the RK point, we computed the gap of the dressed vison states in perturbation theory to first order in ͑1-V͒. However, at this order, the gap turns out to close very slowly and does not lead to a meaningful estimate for the critical V at the liquid-crystal transition. This failure is closely related to the remark above: the size of the core of the vison presumably grows rapidly away from the RK point, a feature which cannot be accounted for with the present variational states.

IV. CONCLUSIONS

In this paper, we have developed two simple approaches to describe the vison excitations of the QDM on the triangular lattice. The first one is based on a soft-dimer version of the model, which exactly takes the form of a Z 2 gauge theory. We have shown that a semiclassical spin-wave approximation to the spectrum of the dual Ising theory captures the important fact that the disappearance of the RVB liquid is due to a vison condensation. 30 It also reproduces the qualitative shape of the vison spectrum in the disordered phase and, more importantly, at the transition ͑linear spectrum at the correct points of the Brillouin zone͒, as well as the spatial pattern of the ordered crystalline state.

The second approach is a variational approximation to the vison wave functions at the RK point of the ͑hard-core͒ QDM. It reproduces semiquantitatively the vison dispersion relation, and provides a simple picture for the vison wave functions which goes beyond the naive point-vison approximation.

Beyond the problem of vison dispersion and condensation, we expect these approaches to be useful for other problems related to vison excitations in QDM, in particular, their mutual interaction and their interaction with vacancies and mobile holes. This is left for future investigations.

The goal of this appendix is to derive an effective Hamiltonian for the model of Sec. II A defined by the Hamiltonian

H = H J + H ⌫ = -J ͚ l l x -⌫ ͚ i ͟ l͑i͒ l͑i͒ z ͑A1͒
in the limit ⌫ / J Ӷ 1 to fourth order in ⌫ / J. For simplicity, let us define the unperturbed Hamiltonian H 0 ϵ H J and the per-turbation V ϵ H ⌫ . Since the Hilbert space of the model is restricted by the constraint that the number of dimers starting from a given site is odd, the ground state manifold of the unperturbed Hamiltonian for positive J consists of all states having exactly one dimer emanating from each site. Let us denote by P 0 the projector onto the Hilbert space generated by these states, and by S the resolvent defined by

S = - 1 -P 0 H 0 -E 0 . ͑A2͒
It is easy to check that V changes the parity of the number of dimers. Indeed, the term of V acting on the triangle i transforms the states with 0 and 3 ͑1 and 2͒ dimers around the triangle i into each other. This implies that P 0 VP 0 = 0 and, more generally, that the effective Hamiltonian only contains even powers of the perturbation. Thanks to the property P 0 VP 0 = 0, the fourth order contribution reduces to three terms, and the effective Hamiltonian up to fourth order reads

H eff = P 0 VSVP 0 + P 0 VSVSVSVP 0 - 1 2 P 0 ͑VS 2 VP 0 VSV + VSVP 0 VS 2 V͒P 0 .
Up to a constant, this effective Hamiltonian can be written as a QDM acting on four-and six-site plaquettes. The four-site Hamiltonian has the form of the regular RK model with t = ⌫ 2 / J -⌫ 4 / J 3 and V = ⌫ 4 / 2J 3 . The six-site Hamiltonian only consists of kinetic terms that flip the dimers around the three possible types of six-site plaquettes shown in Fig. 10, with amplitude -3⌫ 4 / 4J 3 for types ͑1͒ and ͑2͒, and with amplitude -⌫ 4 / J 3 for type ͑3͒.

APPENDIX B: HOPPING AMPLITUDE FOR THE DRESSED VISONS

For the dressed vison states ͓Eq. ͑58͔͒, the equivalent of Eq. ͑54͒ is

H ij ␣ = 2 N ͚ r 0 ͚ ͑c,c ¯͒ F r 0 ͑c͒=1
͑-1͒ N͑c,i͒+N͑c,j͒ ͓͗c͉͑1 + ␣F ˆi͒ + ␦ i,r 0 ͗c ¯͉͑1 + ␣F ˆi͔͒⌸ ˆr0 ͓͑1 + ␣F ˆj͉͒c͘ + ␦ j,r 0 ͑1 + ␣F ˆj͉͒c ¯͔͘.

͑B1͒

Using the explicit form of the projector ⌸ ˆr, we get ͑-1͒ N͑c,i͒+N͑c,j͒ ͕1 + ␣F i ͑c͒ -␦ i,r 0 ͓1 + ␣F i ͑c ¯͔͖͕͒1 + ␣F j ͑c͒ -␦ j,r 0 ͓1 + ␣F j ͑c ¯͔͖͒. ͑B2͒

Unlike Eq. ͑56b͒, both the coverings c and c ¯enter the expression. This does not have the form of a diagonal observable, and the terms F r ͑c ¯͒ ͓with r ͕r 1 ͑i͒ , r 2 ͑i͒ , r 3 ͑i͒ , r 1 ͑j͒ , r 2 ͑j͒ , r 3 ͑j͖͔͒ need to be eliminated to allow for an evaluation with Pfaffians. By inspecting the possible relative positions of two rhombi r 0 and r, one arrives at the following two relations:

͑1͒ If two rhombi r 0 and r have zero, one, three, or four sites in common, F r ͑c ¯͒ = F r ͑c͒ for any pair of configurations ͑c , c ¯͒ which differ by a flip around the rhombus r 0 . ͑2͒ It r 0 and r have two sites in common, let us call b the bond of r which does not touch r 0 . Then we have F r ͑c ¯͒ = D b ͑c͓͒1-F r ͑c͔͒, where D b ͑c͒ =1 if b is occupied by a dimer of c, and 0 otherwise.

We may combine the two cases above into some compact notation F r ͑c ¯͒ = A r,r 0 ͑c͒, valid for any pair of configurations ͑c , c ¯͒ which differ by a flip around the rhombus r 0 . Accordingly, we define an operator A ˆi,r 0 = A ˆr1 ͑i͒,r 0 + A ˆr2 ͑i͒,r 0 + A ˆr3 ͑i͒,r 0 for each triangle i and rhombus r 0 . The matrix element of Eq. ͑B1͒ is now expressed as the expectation value of a diagonal operator:

H ij ␣ = ⑀ ij 2 ͚ r 0
͗͑-1͒ N ˆ͑i,j͒ F ˆr0 ͑1 + ␣F ˆi -␦ i,r 0 ͑1 + ␣A i,r 0 ͒͒͑1 + ␣F ˆj -␦ j,r 0 ͑1 + ␣A j,r 0 ͒͒͘.

͑B3͒

This expression has to be expanded into a polynomial in F ând D ˆoperators before each term can be evaluated thanks to the Pfaffian of an appropriate Kasteleyn matrix. As before the ͑-1͒ N ˆ͑i,j͒ introduces some sign changes, and each F ˆ͑or D ˆ͒ operator requires isolating the corresponding rhombus ͑or bond͒ by switching to zero the corresponding matrix elements. Several tens of terms typically appear for each pair of triangles ͑ij͒, and an automated treatment by computer had to be coded to obtain the hopping amplitudes. The results of Fig. 9 represent several hundreds of CPU hours using the software MAPLE to generate all the correlators and evaluate the corresponding Pfaffians on a finite lattice with 28ϫ 28 sites. *gregoire.misguich@cea.fr Ising models and the Z 2 gauge theories appearing in the description of Mott insulators. 27 In the kagome lattice ͑and on any lattice made of corner sharing triangles, in general͒, this can be done explicitly in a simple way ͑Ref. 2͒. 28 The fact that the Hamiltonian describing the vison hopping is not translation invariant is due to the fact that the dimer liquid has a nontrivial projective symmetry group ͑Ref. 25͒. 29 In the kagome-lattice model of Ref. 2, this is an exact eigenstate. 30 In the frustrated Ising model formulation of the dimer model, x measures the presence ͑-1͒ or absence ͑+1͒ of a vison, and z creates or annihilates a vison ͓see Eq. ͑39b͒, for instance͔. The classical ground state turns out to have ͗ z ͘ = 0 in the dimer liquid phase and ͗ z ͘ 0 in the crystal phase ͑Sec. II B 2͒. So, the vison creation or annihilation operator acquires a finite expectation value at the transition, the usual signature of a particle condensation. This is almost identical to the standard hard-coreboson ↔ spin-1 2 correspondence. If the particle density is represented by 1 2 ͑1-x ͒ ͑as here for visons͒, magnetic long-range order in the z ͑or y͒ direction ͑for the spin variables͒ is equivalent to Bose condensation ͑off-diagonal long-range order in the Boson operators͒. An important difference is, however, that the number of visons is not conserved. Only their parity is con-

FIG. 2 .FIG. 3 .

 23 FIG.2. ͑Color online͒ Dashed lines: Horizontal paths used in Eq. ͑6͒ connecting the boundary of the system to triangles i and j. Dotted segment: Path ⍀ ij . In this example, ⍀ ij crosses one dimer of the reference configuration c 0 ͑ellipses͒ and, therefore, ⑀ ij =-1 ͓Eq. ͑47͔͒.

FIG. 4 .

 4 FIG. 4. ͑Color online͒ Plot of the rescaled relative dimer density D ij in one of the 48 classical ground states of the FFIM in transverse field, for a transverse field ⌫ just below ⌫ c . Thin red bonds represent D ij =0 ͑corresponding to the highest dimer density d ij = 0.5͒. Blue bonds are for D ij Ͻ 0, that is, a lower dimer density, with a width proportional to ͉D ij ͉. D ij takes only four different values: 0, -0.259, -0.518, and -0.776.

2 2 2 FIG. 5 .FIG. 6 .

 256 FIG. 5. ͑Color online͒ Spin-wave dispersion relation S⑀ k / J ͓Eq. ͑34͔͒ for the path A → B → C → A in the Brillouin zone ͑see Fig. 6͒ and different values of ⌫ / J: 3 ͑top͒, 2.75, 2.5, and ͱ 6 Ӎ 2.448 ͑bot- tom͒. At the critical point ⌫ / J = ͱ 6, the spectrum is linear around k = ͑ / 6, / 2͒ and k = ͑5 / 6, / 2͒.

FIG. 9 .

 9 FIG. 9. ͑Color online͒ Variational vison ͓Eq. ͑58͔͒ dispersion relation for different values of ␣, along the path A → B → C → A in the Brillouin zone ͑see Fig. 6͒. The curve for ␣ = 0 corresponds to the variational energies of point-vison states ͓Eq. ͑39a͔͒. The absolute minimum is E min = 0.119, found at ͑k x , k y ͒ = ͑ / 6, / 2͒ = B for a value Ӎ-0.8 of the variational parameter ␣. The exact value of the gap ͑0.089͒ ͑Ref. 6͒ is marked as a dotted horizontal line. In this calculation, the matrix elements S ij ␣ are neglected for triangles i and

FIG. 10 .

 10 FIG.10. The three types of six-site plaquettes around which dimer shifts are generated at fourth order in ⌫ / J.

ACKNOWLEDGMENTS

We are very grateful to F. Becca, M. Ferrero, D. Ivanov, V. Pasquier, and A. Ralko for numerous insightful discussions. F.M. acknowledges the financial support of the Swiss National Fund, of MaNEP, and of the Région Midi-Pyrénées through its Chaire d'Excellence Pierre de Fermat program, and thanks the SPhT ͑Saclay͒ and the Laboratoire de Physique Théorique ͑Toulouse͒ for their hospitality.