Ezgi Berberoglu

Christian T Stoeck

Sebastian Kozerke

Martin Genet

VALIDATION OF CARDIAC STRAIN ESTIMATION FROM 3D-TAGGED MAGNETIC RESONANCE IMAGES USING FINITE ELEMENT IMAGE CORRELATION

Validation of cardiac strain estimation from 3D-tagged magnetic resonance images using finite element image correlation

Introduction

Recent clinical and experimental studies have revealed that compared to global properties, regional changes allow for diagnosis and intervention in many cardiac diseases [START_REF] Genet | Tissue Tracking Technology for Assessing Cardiac Mechanics[END_REF]. Speckletracking echocardiography and cardiovascular magnetic resonance (CMR) feature tracking are the most common techniques to assess regional cardiac motion. In addition, CMR tagging techniques are available in the research setting. Although global strain components are shown to be reproducible for both CMR approaches, regional heterogeneity is not detected robustly [START_REF] Genet | Tissue Tracking Technology for Assessing Cardiac Mechanics[END_REF]. Moreover, strain computation has been shown to underestimate the radial component [START_REF] Tobon-Gomez | Benchmarking Framework for Myocardial Tracking and Deformation Algorithms: An Open Access Database[END_REF]. The objective of the present work is to extend and validate the recently proposed finite element digital image correlation method by Genet et al. for strain quantification from 3D tagged images [START_REF] Genet | Finite Element Digital Image Correlation for Cardiac Strain Analysis from 3D Whole-Heart Tagging[END_REF].

Methods

The finite element-based digital image correlation tool used in this study [START_REF] Genet | Finite Element Digital Image Correlation for Cardiac Strain Analysis from 3D Whole-Heart Tagging[END_REF] utilizes a mechanical regularizer, penalizing any deviation from mechanical equilibrium for strain computation. In this work, what we improve is to project the images onto the mesh, and then interpolate and differentiate, instead of simply interpolating image intensities and gradients at the mesh integration points [START_REF] Genet | Finite Element Digital Image Correlation for Cardiac Strain Analysis from 3D Whole-Heart Tagging[END_REF]. This allows for a higher robustness with respect to image discretization by taking into account only the voxels within the mesh. Then, synthetic 3Dtagged images are generated from a reference biomechanical model of the left ventricle [START_REF] Imperiale | Constitutive Parameter Estimation Using Tagged-MRI Data[END_REF][START_REF] Sainte-Marie | Modeling and Estimation of the Cardiac Electromechanical Activity[END_REF]. The proposed image correlation method is run on the synthetic tagged images, and the obtained strains are compared to the ground truth of the biomechanical model. Additionally, the effect of image resolution and noise are investigated to understand the limits of strain analysis on tagged images.

Results

We validate our image correlation tool on tagged images after demonstrating that strain components for synthetic images match with the reference results (Figure 1). Then, we run analysis on resampled images and observed a deviation from the reference strain as we increase the downsampling factor. Compared to the reference, using downsampling factors of 2, 4, and 16 results in a deviation of peak radial strain by 2%, 6%, and 10%, respectively while the other strain components do not change significantly. We also observe a similar trend when we run our simulations using images with different signal-to-noise ratio.


PhD Laurent Hoffer

PhD Christophe Muller

PhD Xavier Morelli
email: xavier.morelli@inserm.fr

Gisbert Schneider
email: gisbert.schneider@pharma.ethz.ch

Matthias Rarey
email: rarey@zbh.uni-hamburg.de

Philippe Roche
email: philippe.roche@inserm.fr

Chemistry-Driven Hit-To-Lead Optimization Guided by Structure-Based Approaches

Keywords: hit-to-lead optimization, structure-based drug design, de novo design, library design, virtual screening Biological and Medicinal Chemistry, Drug design, Structure-activity relationships, Virtual screening hit-to-lead optimization, structure-based drug design, de novo design, library design, virtual screening

For several decades, hit identification for drug discovery has been facilitated by developments in both fragment-based and high-throughput screening technologies. However, a major bottleneck in drug discovery projects continues to be the optimization of primary hits from screening campaigns to derive lead compounds. Computational chemistry or molecular modeling can play an important role during this hit-to-lead (H2L) stage by both suggesting putative optimizations and decreasing the number of compounds to be synthesized and evaluated. However, it is also crucial to consider the feasibility of organically synthesizing these virtually designed compounds. Furthermore, the generated molecules should have reasonable physicochemical properties and be medicinally relevant. This review focuses on chemistry-driven and structure-based computational methods that can be used to tackle the difficult problem of H2L optimization, with emphasis being placed on the strategy developed in our laboratory.

Introduction

The first step in a drug discovery project involves the identification of hit compounds that typically exhibit weak to moderate affinity for the biological target. [1] These primary hits are usually discovered either via the high-throughput screening (HTS) of large collections of diverse molecules that are of medium complexity or the use of a fragment-based technology in which small chemical libraries of low-molecular-weight fragments are considered. [2] There is a current trend toward fragment-based drug discovery (FBDD), since the optimization phase of hit compounds from HTS can be fastidious while also maintaining a reasonable molecular weight. [3] Moreover, FBDD approaches can cover a larger chemical space in the optimization phase, which leads to structural novelty and a higher probability of success. An increasing number of success stories have been reported in the past decade, and more than 30 compounds derived from fragments are currently in clinics or on the market. [4] Following orthogonal validation, the confirmed fragments/hits are improved in a subsequent stage called hit-to-lead (H2L). Typical H2L involves chemical modifications around the validated hit to optimize its affinity for the target to become a lead compound. [5] These optimization phases can be completed by exploiting a trial-and-error strategy, and several cycles are usually needed to reach a suitable affinity. Moreover, other features may be preferred, such as improving the physicochemical properties of the compounds and enabling or maintaining a degree of selectivity with respect to undesired related targets. In practice, successful H2L optimization can improve the binding constants by several orders of magnitude. However, because of the numerous parameters that must be considered, the H2L process can also be unsuccessful, time-consuming, and highly expensive. [6] 4 Various computer-based approaches have been developed to overcome bottlenecks during H2L. [7] For instance, ligand-based methods, such as the Quantitative Structure-Activity Relationship (QSAR), can be used to optimize a series of compounds by exploiting available experimental data and calculated descriptors. [8] However, this review will focus on computational structure-based approaches [9] , such as molecular docking, de novo design (DnD) and pharmacophore screening. These methods rely on structural data to generate models, and success stories have recently been reported. [10] In practice, X-ray crystallography is the primary approach to determine the binding mode of a compound in its binding site at the atomic level. [11] Such structural data greatly facilitates the H2L optimization by clearly identifying nearby protein sub-pockets around the engaged fragment that could be used in the optimization phase.

The most intuitive approach for optimizing a fragment using an in silico structure-based approach is the transposition of the growing paradigm from FBDD to the molecular modeling context: The affinity of a compound for its target is increased by adding chemical moieties that are able to create new favorable contacts while maintaining its original binding mode (Figure 1A). Two additional H2L strategies used in FBDD, namely, linking and merging, can also be virtually mimicked (Figures 1B &1C). [12] In contrast to the growing concept where a single fragment is required, merging and linking consist in the covalent assembly of two nonoverlapping fragments, either directly or via a spacer of variable length. [13] Merging and linking strategies are less frequently used than growing and are also more challenging because they require two hits and the conservation of their original orientation after their fusion into a single compound. However, when these strict criteria are validated, merging and linking can lead to outstanding improvements in affinity. [14] In practice, these two approaches are primarily applicable to fragments because merging or linking two compounds with moderate to high molecular weight would result in molecules that are too large to be used as probes or drugs.

Computational-based tools that are able to address growing, linking or merging strategies are of primary interest during the H2L stage. [12,13,15] For instance, a virtual-focused library that is around the validated fragment can be generated within the first step. Subsequently, this library can be virtually screened using a constrained docking strategy to mimic the growing paradigm and maintain the original interactions. The best putative optimizations are finally selected in the top hit list using a scoring function. Alternatively, classical "unconstrained docking" of the focused library can be performed. Pharmacophore filters [16] or interaction fingerprints [17] are subsequently used to extract putative optimizations that maintain the original binding mode.

Finally, de novo design algorithms [18,19] , with the binding mode of the validated fragment serving as the starting seed, can also be used to sample the cavity and generate putative optimized compounds.

A critical point is determining how to handle the creation of new covalent bonds made during the virtual H2L optimization step. A first attempt to generate reasonable structures has been made by the RECAP [20] or BRICS [21] methods. For instance, molecules in RECAP are fragmented around specific bond types, and new terminal atoms are flagged to capture their previous chemical environment. Virtual fragments with complementary flags can subsequently be merged to design new virtual compounds. However, there is no assessment regarding their synthesis, since neither organic chemistry rules nor the availability of building blocks (BBs) were used during the process. This approach leads to either time-consuming follow-up efforts to devise synthesis routes or even the inability to produce the virtually generated compounds.

Consequently, it is of critical importance to address the synthesis tractability of virtually generated compounds in the context of prospective drug discovery projects. [18,22] In this effort, generic reaction schemes "reactants → products" can be encoded to mimic common organic synthesis routes. [23,24] These chemistry rules can subsequently be retained during the design of new molecules to facilitate the synthesis stage.

Several remarkable computational structure-based H2L methods were described in the past decade, but they lacked this fundamental requirement of ensuring the synthesizability of the generated compounds for use in drug discovery projects. [25,26] De novo design or alternative structure-based H2L methods, which approached this critical issue, are presented in this review, with an emphasis on the in silico steps of the DOTS strategy developed in our laboratory. [27]

Virtual H2L methods relying on the "de novo design" concept

By definition, de novo design tools can automatically build compounds "from scratch" within a binding site of known 3D structure using predefined sets of substructures and rules governing their linkage. [18,19] In theory, DnD methods enable the exploration of a considerably larger chemical space, which is in contrast to classic virtual screening (VS), where only commercially available compounds are tested. Although conceptually similar, DnD tools can be distinguished by their algorithms, ranking/scoring functions, convergence criteria, and branch-pruning strategies. This methodological approach, developed in the early 1990s, was initially appealing but is currently not widely practiced by chemoinformaticians. [28,29] Indeed, DnD methods often exhibit several critical drawbacks that limit their application to prospective cases. The primary drawbacks include the following: 1) low reliability in predicting affinity, 2) potentially poor physicochemical properties of the designed compounds and, most importantly, 3) issues regarding the compounds' synthetic tractability. [18] As a response to the first criticism, in silico H2L methods that rely on DnD approaches, as introduced below, were upgraded to use a pre-positioned fragment in the binding site as the original seed. Starting from an experimentally validated fragment with a known binding mode is likely to increase the reliability of their predictions.

The third criticism of DnD approaches -the major drawback -concerns the difficulty in synthesizing certain of the suggested compounds because the rules governing bond creation during the process did not take into account any organic chemistry knowledge. This primary issue was approached by incorporating virtual reaction schemes in the workflow to ensure the synthetic accessibility with more confidence. Both ligand-based and structure-based DnD tools were developed to address this problem. DOGS [24] and SYNOPSIS [30] are popular examples of such chemistry-driven ligand-based DnD methods. However, other structure-based DnD approaches that rely on both substructure seed prepositioned in the binding site and encoded chemistry knowledge exist and are discussed below.

The SPROUT program was one of the first published DnD software in the 1990s, but several major upgrades have been recently added. [29] This program is now able to take a validated hit or fragment as seed and to incorporate chemistry-based rules during the virtual design process. SPROUT was successfully used to design inhibitors of dihydroorotate dehydrogenase enzyme [31] , and more recently BACE-1 inhibitors, based on a synthetic Suzuki reaction scheme. [32] LigBuilder [33] is another validated tool that was recently updated to tackle the usual drawbacks from DnD software. The current version (2.0), is able to perform both ab initio design and lead optimization of a compound of interest while considering both synthesis accessibility and filter-based drug-likeness of the designed compounds. [34] LigBuilder was successfully applied to generate nanomolar inhibitors of Cyclophilin A. [35] Similarly, the AutoGrow tool [36] was also updated to include organic chemistry rules within the operators of its evolutionary-based algorithm. [37] In this instance, the popular Autodock Vina docking engine is used for both conformational sampling and scoring stages. [38] Physicochemical filters were also added to reject non drug-like compounds.

Beccari et al developed a suite of programs called LiGen for DnD that both handles the chemical rules for designing accessible compounds and efficient docking with pharmacophore constraints. [39] Each module can be used either separately or combined in a complex workflow procedure. This toolbox is reported to tackle most issues in the DnD field while being efficient from a computational point of view. However, there is no case study to highlight its ability in either retrospective or prospective cases.

Cheron et al reported the development of OpenGrowth for the computer-based H2L optimization of compounds under binding site constraints. [40] In this study, the probability that a given fragment is connected to another one in a reference drug dataset is used to guide the design process. At the end of the process, designed molecules exhibit both reasonable synthetic accessibility and good physicochemical properties, although organic synthesis routes are not explicitly taken into account. A new version, including the implementation of an explicit synthetic accessibility score during the design process, is planned.

It should be noted that many DnD tools were only validated in retrospective case studies. However, all these programs exhibit promise and should be useful for prospective projects aiming to design active and accessible compounds.

Alternative approaches for in silico H2L optimization

Chevillard et al optimized several fragment-like compounds into low micromolar ligands for the β2-adrenergic receptor target. [41] The PINGUI method relies on both structural data and a series of 58 encoded organic chemistry rules published by Hartenfeller et al for designing putative accessible ligands by merging an original fragment with compatible BBs. [23] In practice, the SEED program is used to dock pre-processed BBs (called surrogates) where the reactive center is modified according to the considered synthesis reaction. [42] For instance, a methylamine group will replace the aldehyde function from the BB if the reductive amination synthesis scheme is selected. In a second step, putatively interesting surrogates, where the reactive center is sufficiently close to the one from the original hit without any overlap with other atoms, are identified. Next, final products corresponding to the coupling of the hit and selected surrogates are generated and docked using the DOCK software. [43] A last modeling step involving the refinement and rescoring of selected poses using the SZYBKI method is performed. [44] The LeadOp+R method was developed to perform structure-based H2L optimization with synthetic accessibility. [45] This method relies on approximately 200 encoded chemical reactions while allowing multi-steps design. The pipeline starts with a query structure pre-positioned in the binding site and user-defined preferred ligand-receptor interactions. New 1D/2D compounds are created by combining the current structure with BBs by using encoded chemistry rules without knowledge of the binding site at this stage. 3D conformers are subsequently generated for each compound, superimposed to the shared reference substructure in the binding site, and evaluated for their ability to make additional favorable contacts without any clashes. The cycle stops when required interactions between both entities are fulfilled. Several physicochemical filters are used to discard compounds with undesirable properties. They successfully applied their approach on two retrospective projects: several known potent inhibitors were designed for each target, while suggested synthetic routes shared steps with published ones.

An original strategy, relying on both pharmacophore and docking concepts, allowed Schulz et al to discover several covalent inhibitors of the enteroviral 3C protease target. [46] First, a new pharmacophore feature type in the LigandScout software was developed to look for predefined reactive functions in the vicinity of a cysteine residue. [47] Next, a 3D-pharmacophore was designed to catch fundamental interactions with the enteroviral 3C protease while looking for covalent binders using this new feature. Then, a library of fragment-like compounds was screened against the 3D-pharmacophore. Several fragment hits were identified and experimentally validated using mass spectrometry. The best one was further investigated but exhibited some instability. To bypass this issue, a scaffold hopping strategy, using simple SMARTS-based substructure search, was successfully employed to find alternative binding cores that still contain the required features. A protocol was later developed for the rational design of optimized analogs of the best new hit, while maintaining the fundamental interactions with the protein. Thus, a virtual library was generated by coupling the hit and commercially available BBs using encoded chemistry reaction rules. LigandScout and GOLD programs were used to identify putative optimizations that would occupy adjacent sub-pockets. Finally, the authors reported both reversible and irreversible inhibitors of the enteroviral 3C protease target.

Evers et al from the Sanofi-Aventis company reported the CROSS method for either rescaffolding or the optimization of compounds using explicit handling of organic reactions. [48] The CROSS approach relies on the BROOD software to quickly identify pre-processed fragments that could replace an undesirable core by using "exit vectors" and 3D-shape analysis. [49] The main advantage is that these pre-processed fragments are directly connected to available BBs and specific chemical reactions for easier continued investigation. In addition, the use of chemical protection groups is also allowed to facilitate the results in the synthesis steps. Finally, generated compounds are post-processed using molecular docking and ADME-Tox predictions. While the main usage is rescaffolding, the CROSS method can also be employed in virtual H2L optimization using either a growing or linking strategy. Indeed, linking and rescaffolding are similar concepts, where a linker is selected to connect two moieties while maintaining their original orientation.

AutoCouple was recently reported as a useful tool to expand the chemical space in hit optimization. [50] In the first step, a diversity-oriented library is designed by virtually coupling one hit that includes a reactive function, with a list of commercially available BBs using encoded chemistry rules. In a second step, the rDock program is used to screen the library while adding constraints on the original moiety to maintain the reference binding mode during the conformational sampling stage. [51] The authors applied this strategy to the design of potent CBP bromodomain inhibitors. An acetyl benzene moiety served as the reference substructure able to mimic the acetylated lysine from histone tails, which are recognized by this epigenetic reader. Different organic reactions were considered to build the virtual library that was docked by rDock in the CBP binding site with constraints on the reference moiety. The best poses were minimized using the CHARMM molecular mechanic program [52] during a post-processing stage, before the final selection of target compounds to be synthesized was made. This study led to the discovery of several nanomolar inhibitors for the CBP target. The predicted binding mode of several compounds was successfully confirmed by solving the structure of the complex using X-ray crystallography.

Finally, the Diversity-Oriented Target-focused Synthesis (DOTS) is an integrated strategy developed in our laboratory for generic H2L optimization relying on the growing paradigm. [27] This strategy involves molecular modeling (chemical library design and structure-based VS) and robotic-based experimental stages (diversity-oriented de novo synthesis and in vitro evaluation and validation). The virtual steps are an upgrade from a former computational optimization method [26] relying on the RECAP algorithm [20] , which did not consider the synthetic accessibility of the designed compounds. The general DOTS workflow can be summarized as follows (Figure 2): 1) hit identification and characterization of its binding mode using structural biophysics method, such as X-ray crystallography; 2) design of a virtual focused library around the hit using a database of commercially available BBs, encoded organic chemical rules, and post-processing of the library to extract a diverse subset of representative compounds that also possess reasonable physicochemical properties without any undesirable functions to medicinal chemists; 3) constrained VS of the library with the S4MPLE tool [26,53] to identify the best putative optimizations that create additional favorable contacts while maintaining the original binding mode; 4) parallel synthesis of the compounds using a chemistry robot; and 5) in vitro evaluation with a robotic screening workstation. The in silico part of DOTS relies on real chemical knowledge and allows for the production of theoretically accessible and diverse compounds, while exploring the chemical space around the hit, and matching regular physicochemical and medicinal chemistry-like features. The DOTS strategy was successfully applied to the optimization of a previously reported xanthine core that binds the first bromodomain (BD1) of the BRD4 protein. [54] Several sub-micromolar inhibitors were developed and validated in one cycle of optimization, with the best one displaying a K d value of 190 nM. 13 X-ray crystallography was used to solve the structure of the best inhibitor in complex with the BRD4(BD1) protein, and confirmed the predicted binding mode. Several features, including linking optimization and design of covalent inhibitors, are currently under development in order to address all available H2L strategies.

Summary and Outlook section

This review contributes to an overview of structure-based computational H2L methods that tackle the synthetic accessibility of virtually generated compounds. The different programs discussed in this review, and their main characteristics, are summarized in Table 1. Most of them currently encode real chemistry knowledge to design molecules that can be synthesized with high probability in one or two steps. Such approaches help scientists design accessible and optimized compounds that perfectly fit the binding site of the target.

Despite the shared ability to handle virtual chemistry, these in silico structure-based H2L methods rely on various strategies. De novo design algorithms start from a given compound as the original seed to perform the optimization process, while alternative approaches use docking or pharmacophore concepts to identify promising compounds from a virtual focused library that was designed around the hit to optimize.

In certain cases, a single reaction scheme is considered, but it is still possible to explore a large chemical space around the starting hit due to the high number of commercially available BBs. Moreover, target compounds selected using these approaches can also be produced in parallel using robotic platforms when a single reaction is considered, as exemplified in the DOTS approach. inhibitors. [55] B. Linking example leading to the discovery of a Pyruvate kinase (PKM2) activator. [56] C. Merging example towards the development of small molecule inhibitors of Mycobacterium tuberculosis transcriptional repressor protein (EthR). [57] Figure 2. Schematic workflow of the DOTS strategy. 1) Following hit identification, the binding mode is characterized using structural biophysics methods such as X-ray crystallography. 2) A virtual focused library is conceived by combining a database of functionalized BBs with an activated form of the original hit using SMARTS-encoded medicinal chemistry-relevant organic synthesis rules. The raw library is then filtered to extract a diverse set of compounds with reasonable physicochemical properties. 3) The focused library is virtually screened under constraints with S4MPLE to identify compounds that create additional favorable contacts while maintaining the original binding mode. 4) The selected compounds are synthesized using an automated parallel chemistry robot. 5) Compounds are finally evaluated in vitro with a robotic screening workstation and the best molecules are further validated using orthogonal methods.

Figure 1 .

 1 Figure 1. Selected examples of fragment growing, linking and merging strategies. 2D

Acknowledgments

We would like to thank Prof. A. Varnek and Dr. D. Rognan for the organization of the 2018 International Summer School in Chemoinformatics.

Computational H2L can be useful to optimize hits in a time-efficient and cost-effective manner, as highlighted by successful cases described in this review. These in silico methods clearly should play a larger role in drug discovery in both academic and pharmaceutical environments.

Author Comments: Table 1. Hit to lead algorithms discussed in this review. The different approaches can be divided into two major categories according to the method used to optimize the original hit, de novo design (chapter 2) and virtual screening (chapter 3). All these methods rely on explicit synthesis rules to design new compounds except for OpenGrowth and they all use medicinal chemistry--like filters to discard molecules with undesirable properties.


François Vernay

Arnaud Ralko

Federico Becca

Frédéric Mila

Francois Vernay

Identification of an RVB liquid phase in a quantum dimer model with competing kinetic terms

Keywords: number͑s͒: 75.10.Jm, 05.50.ϩq, 05.30.Ϫd

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

I. INTRODUCTION

After their first derivation by Rokhsar and Kivelson in 1988 in the context of cuprates, 1 the hard-core quantum dimer models ͑QDM͒ have attracted significant attention. The phase diagrams of the QDM on the square and triangular lattices have been investigated in great details, 2,3 and, following the pioneering work of Moessner and Sondhi on the triangular lattice, 4 the very existence of a stable resonating valence bond ͑RVB͒ phase has been unambiguously demonstrated. 5 The presence of a liquid phase with deconfined vision excitations 6 has also been established for a toy model living on the Kagome lattice. 7 However, the relationship between QDMs and Mott insulators, the physical systems for which they were proposed in the first place, is not straightforward. It is well established by now that the ground state of the S = 1 2 Heisenberg model on the square and triangular lattices exhibits long-range magnetic order of Néel and 120 degree type, respectively, and this type of order cannot be reached within the variational basis of Rokhsar and Kivelson, which consists of short-range singlet dimers. For a QDM to be a good effective model, one should thus identify models for which the subspace of shortrange dimer coverings on a certain lattice is a good variational basis.

The first example of such a case was provided by the S = 1 2 Heisenberg model on the trimerized Kagome lattice. 8 Indeed, an effective spin-chirality model living on a triangular lattice can be derived, and, at the level of a mean field decoupling between spin and chirality, the ground state manifold consists of all dimer coverings on the triangular lattice. Going beyond mean field is thus expected to lead to a relevant effective QDM. Using Rokhsar and Kivelson's prescription, which consists in truncating the Hamiltonian and inverting the overlap matrix within the basis of dimer coverings, Zhitomirsky has derived such an effective Hamiltonian 9 and shown that the main competition is between kinetic terms involving loops of length 4 and 6, respectively, and not a competition between a kinetic and a potential term, as for the QDM derived by Rokhsar and Kivelson. The next logical step would be to study the properties of this effective QDM. This is far from easy however. We know from the experience with the standard QDM model on the triangular lattice that the clusters reachable with exact diagonalizations are much too small to allow any significant conclusion regarding the presence of an RVB phase, and since there is no convention leading only to negative off-diagonal matrix elements, it is impossible to perform quantum Monte Carlo simulations.

Recently, two of us came across another model, for which the low energy sector consists of almost degenerate singlet coverings on the triangular lattice. This model is a Kugel-Khomskii 10 model that was derived in the context of LiNiO 2 , and the mean-field equations that describe the decoupling of the spin and orbital degrees of freedom possess an infinite number of locally stable solutions. These solutions are almost degenerate and correspond to spin singlet ͑and orbital triplet͒ dimers on the triangular lattice. 11 Following Rokhsar and Kivelson's prescription, an effective QDM can also be derived ͑see the Appendix͒. As for the S = 1 2 Heisenberg model on the trimerized Kagome lattice, it consists of a competition between kinetic terms, with two important differences however. The main term of length 6 lives on loops that have a shape of large triangles, a term absent in the other case. But more importantly, the off-diagonal matrix elements are all negative.

Since the competition between kinetic processes was never investigated before, we have decided to concentrate on the minimal model obtained by keeping only the dominant term of length 6 for clarity. We have checked that the properties of the complete effective model are similar. This minimal model is described by the Hamiltonian

H = -t ͚ ͉͑//͗͘ --͉ + H.c.͒ -tЈ ͚ ͩͯ/ \ -ʹͳ / - \ ͯ+H. c. 0 0 ͪ + V ͚ ͉͑//͗͘//͉ + ͉ --͗͘ --͉͒,

͑1͒

where the sums run over the four-site and six-site loops with all possible orientations. Although the repulsion is a higher order process, we have included a repulsion term in the Hamiltonian, and we will treat its amplitude V as a free parameter to be able to make contact with the Rokhsar-Kivelson model on the triangular lattice. The hopping amplitudes t and tЈ are negative, and although the ratio tЈ / t is in principle fixed by their expression in the perturbative expansion, we will also treat it as a free parameter.

Our central goal in this paper is to determine the nature of the ground state as a function of tЈ / t. With respect to what we already know about QDMs, the main question is whether a competition between kinetic terms can also lead to a liquid phase. As we shall see, the answer to that question is positive, a liquid phase being present in a finite region of the phase diagram in the tЈ-V plane.

The paper is organized as follows. In Sec II, we briefly review the basic preliminaries used in the rest of the paper. The results obtained with exact diagonalizations are presented in Sec. III, those obtained with quantum Monte Carlo in Sec. IV, and the conclusions in Sec. V. The perturbation calculation which has motivated the investigation of this QDM is finally presented in the Appendix.

II. THE METHOD

In this section, we present a brief introduction to the numerical methods, to the clusters used in the analysis, and to the physical concepts underlying the determination of the phase diagram. More details can be found in Ref. 5.

Let us first discuss the shape of the finite-size clusters. In general, a finite cluster is defined by two vectors T 1 and T 2 and, in order to have the symmetries for rotations by 2 /3, they must satisfy 12

T 1 = lu 1 + mu 2 , T 2 = -mu 1 + ͑l + m͒u 2 ,
where l and m are integers and u 1 = ͑1,0͒ and u 2 = ͑1/2, ͱ 3/2͒ are the unitary vectors defining the triangular lattice. The number of sites in the cluster is N = l 2 + m 2 + lm. In order to have also the axial symmetry, and therefore all the symmetries of the infinite lattice, we must take either lm =0 or l = m. The first possibility corresponds to type-A clusters ͑with the notation of Ref. 5͒, with N = l 2 sites; the second one gives rise to type-B clusters, with N =3ϫ l ϫ l ͑for examples, of both cases, see Fig. 1͒. Since for tЈ / t = 0 and V / t = 0 the ground state belongs to a crystalline phase with a 12-site unit cell, 4,5 we will restrict in this work to clusters with a number of sites multiple of 12, in order not to frustrate this order. To limit the finite size effects related to the geometry of the clusters, we will concentrate on type-B clusters, which are always compatible with this order. Note that the six-dimer loop kinetic term does not introduce further restrictions since all it requires is to be able to accommodate six-site unit cells.

A very important concept in the QDM is the existence of topological sectors. Indeed, in the triangular lattice, the Hilbert space is split into four disconnected topological sectors on a torus defined by the parity of the number of dimers cutting two lines that go around the two axes of the torus and denoted by ͑p , q͒ with p , q =0 ͑respectively, 1͒ if this number is even ͑respectively, odd͒. One can convince oneself, by a direct inspection of the effect of the four-site and six-site terms, that these numbers are conserved quantities under the action of the Hamiltonian of Eq. ͑1͒. More generally, this is a consequence of the fact that the topological sectors are not coupled by any local perturbation. These topological sectors are extremely useful to distinguish between valence bond solids and spin liquids. Indeed, valence bond solids are only consistent with some topological sectors, whereas RVB spinliquid phases are characterized by topological degeneracy. Therefore, the main goal will be to investigate whether, in the thermodynamic limit, the topological sectors are degenerate or not. In that respect, it is useful to remember that the ͑0,1͒ and ͑1,0͒ sectors are always degenerate with either ͑0,0͒ or ͑1,1͒ ͑depending on the cluster geometry͒ since they contain the same configurations rotated by an angle /3. 5 One can thus, without any loss of generality, restrict oneself to the analysis of the ͑0,0͒ and ͑1,1͒ sectors. Therefore, we define the absolute value of the topological gap as

⌬E = ͉E 00 -E 11 ͉, ͑2͒
where E 00 and E 11 are the total ground-state energies for the topological sectors with p = q = 0 and p = q = 1, respectively. This gap is expected to scale to zero with the cluster size in the RVB phase. Finally, in order to detect a possible dimer order, we also consider the static dimer-dimer correlations

D i,j ͑r -rЈ͒ = ͗D i ͑r͒D j ͑rЈ͒͘, ͑3͒
where D i ͑r͒ is the dimer operator defined as follows: It is a diagonal operator in the configuration space that equals 1 if there is a dimer from the site r to the site r + a i , with

a 1 = ͑1,0͒, a 2 = ͑1/2, ͱ 3/2͒, or a 3 = ͑-1/2, ͱ 3/2͒ and van- ishes otherwise.
The method used is the same one as that used by Ralko and collaborators 5 to determine the phase diagram of the Rokhsar-Kivelson QDM on the triangular lattice and our investigations are Lanczos diagonalizations and Green's function Monte Carlo ͑GFMC͒ simulations. In particular, the GFMC is a zero-temperature stochastic technique based on the power method: Starting from a given wave function and by applying powers of the Hamiltonian, the ground state is statistically sampled to extract its energy and equal-time correlation functions. In principle, as in other Monte Carlo algorithms, in order to reduce the statistical fluctuations, it is very useful to consider the importance sampling, through the definition of a suitable guiding function. Unfortunately, when dealing with dimer models, it is very hard to implement an accurate and, at the same time, efficient guiding function for the crystalline phases. This problem is particularly relevant when the six-site term becomes dominant. In this case, our simulations suffer from wild statistical fluctuations, deteriorating the convergence of the GFMC. As a consequence, we are not able to reach the largest available size, 432-site cluster, for all parameters tЈ / t and V / t. By contrast, given the simple form of the spin-liquid ground state ͑that reduces to a superposition of all the configurations with the same weight at the Rokshar-Kivelson point͒, in the disordered region we can use the guiding function with all equal weights and obtain very small fluctuations ͑no fluctuation at the Rokshar-Kivelson point͒ and, therefore, excellent results with zero computational effort. Combining these facts, the loss of the GFMC convergence can be interpreted as a signal for the appearance of a crystalline phase. Of course, this is not a quantitative criterion, but, as it will be shown in the following, it gives reasonable insight into the emergence of a dimer order.

Finally it should be mentioned that, since the different topological sectors are completely decoupled ͑each dimer configuration belonging to one and only one of them͒ and cannot be connected by the terms contained in the Hamiltonian, within the GFMC it is possible to work in a given topological sector, making it possible to extract the groundstate properties of each of them.

III. EXACT DIAGONALIZATIONS

To get a first idea of the properties of the model, we start with the results we have obtained with exact diagonalizations of finite clusters for the model of Eq. ͑1͒ with V / t = 0. Let us first begin with the ground-state energy for the 12-and 36site clusters for both the topological sectors ͑0,0͒ and ͑1,1͒ ͑see Fig. 2͒. Note that the 12-site cluster is of type B, whereas the 36-site one is of type A. We have that, for both sizes, a level crossing occurs for tЈ / t ϳ 2. Below that value, the ground state is in ͑1,1͒ topological sector, in agreement with the earlier results of Ref. 5 for V / t = 0. The main difference between the 12-site and the 36-site clusters is that, for the 36-site cluster, the topological ground-state energies stay very close in a large parameter range for tЈ / t ϳ 2. This could suggest that, upon increasing the size, this level crossing might evolve into a phase where these energies are rigorously degenerate, giving rise to a liquid phase without any crystalline order.

In order to give an idea of the various phases, we report in Fig. 3 the dimer-dimer correlations of Eq. ͑3͒ for the 36-site cluster below, at, and above the level crossing, which takes place at tЈ / t = 2.2 for this size. For small tЈ / t, the correlations show a pattern similar to that of the intermediate ͱ 12ϫ ͱ 12 phase of the standard QDM model, already shown in Ref. 5. It should be stressed that, since in these calculations the Hamiltonian has the translational symmetry, the 12-site unit cell is not directly visible from Fig. 3, and in order to have a clearer evidence one should break the symmetry by hand. Nonetheless, as it has been shown in Ref. 5, these results are in perfect agreement with the existence of a ͱ 12ϫ ͱ 12 phase with a crystalline ground state in the thermodynamic limit. When tЈ / t is large, another pattern arises, which has never been observed in the standard QDM, and which presents a kind of six-site triangle ordering. In this case, the dominant kinetic term involving six sites ͓see the second term of the Hamiltonian ͑1͔͒ induces a dimer pattern with the same symmetry, possibly inducing a ground state with a six-site unit cell in the thermodynamic limit. Also in this case the translational symmetry of the ground state partially masks the existence of a regular dimer pattern. Unfortunately, we will not be able to confirm this prediction since, as stated before, the GFMC algorithm has serious problems of convergence inside this phase and for large clusters.

In any case, for intermediate values of tЈ / t, the correlations decay very rapidly with the distance and are close to those obtained in the liquid phase of the standard QDM with V / t Շ 1 and tЈ / t =0. 5 This fact gives another evidence of the possible existence of an RVB phase between two ordered phases in the model with competing kinetic terms and without the dimer repulsion. Of course by considering the exact diagonalization results only it is impossible to give definite statements on the stabilization of this liquid phase and, therefore, in the following section, we will consider a more systematic study of the topological gap, in order to unveil the existence of a wide disordered region that develops from the Rokshar-Kivelson point V / t = 1 of the standard QDM and survives up to V / t = 0 and finite tЈ / t. FIG. 2. ͑Color online͒ Difference between E 00 and E 11 , the total ground-state energies of the topological sectors with p = q = 0 and p = q = 1, respectively. The results are found by exact diagonalizations for clusters with 12 and 36 sites. FIG. 3. ͑Color online͒ Dimer-dimer correlations on the 36-site cluster. The dimer of reference is the thickest one in the up-right corner. The thicker the line, the farther the value of the correlation from the uniform distribution equal value 1 / 36. Solid lines are used when the correlation is higher than 1 / 36, and dashed lines when it is lower.

IV. GREEN'S FUNCTION MONTE CARLO

In this section, we use the GFMC method to extend the results of the preceding section to larger clusters, with up to 432 sites, and to map out the phase diagram in the tЈ-V plane.

A. The case of V / t =0

Let us first describe the results we have obtained for V / t = 0 and consider the behavior of the topological gap given in Eq. ͑2͒. In Fig. 4, for clarity, we divided the results in two sets for 0 ഛ tЈ / t ഛ 0.6 and tЈ / t ജ 1. The first remarkable feature is that, for most parameters, the gap decreases between 12 and 48 sites, regardless of its behavior for larger sizes. Therefore, the possibility to study much larger clusters is crucial for this analysis. Indeed, there is a clear change of behavior for tЈ / t ϳ 1.6, a value beyond which our extrapolations give solid evidences in favor of a vanishing gap in the thermodynamic limit. On the other hand, for smaller tЈ / t ratios, we have a clear evidence that ⌬E increases for large sizes. Therefore, we come to the important conclusion that, the crystalline ͱ 12ϫ ͱ 12 phase is destroyed by increasing the amplitude of the six-site term and the system is eventually driven into a liquid RVB phase.

To have a further confirmation of the existence of this disordered phase, we have calculated the dimer-dimer correlations. The results obtained on the 108-site cluster for tЈ / t = 2 are reported in Fig. 5, where the correlation functions for parallel dimers along the same row are plotted as a function of the distance. Given the small number of clusters available, a precise size scaling of the order parameter is not possible and also a meaningful estimation of the correlation length is very hard. Nevertheless, the behavior is definitely consistent with an exponential decay and the uncorrelated value of 1 / 36 is approached very rapidly, as expected in a liquid phase without any crystalline order.

Unfortunately, the larger tЈ / t region is numerically far more difficult to access. Indeed, as stated above, although the GFMC is in principle numerically exact, we have not been able to find an efficient guiding function to perform the importance sampling and wild statistical fluctuations prevent us to reach a safe convergence for large clusters. In practice, we have access to clusters up to 108 sites, that are still too small to predict the thermodynamic behavior. For instance, for V / t = 0, the convergence stops before one can observe any increase of the topological gap, and the criterion used for the phase transition on the other side of the RVB phase cannot be used any more. However, the lack of convergence is a clear sign that one enters a new ͑crystalline͒ phase. So, if the change of behavior of the topological gap with the size cannot be observed, we take as a definition of the boundary for the phase transition the parameters for which the convergence is not good any more. This is expected to be semiquantitative, and indeed, as we shall see in the next section when studying the full diagram, the points obtained with this criterion agree reasonably well with those obtained with the reopening of the topological gap.

In summary, although the region where the six-site term dominates over the usual four-site dimer flip is not accessible by using GFMC, based on our numerical results for small tЈ / t, we can safely argue that the crystalline ͱ 12ϫ ͱ 12 phase is destabilized by increasing the six-site kinetic term, leading to a true disordered ground state with topological degeneracy.

B. Phase diagram in the tЈ-V plane

In this section we prove that the RVB phase found in the preceding paragraph for V / t = 0 is connected to the one obtained for the standard QDM, i.e., close to the Rokshar-Kivelson point and tЈ / t = 0. In order to do that, we have investigated a generalization of the previous model that also includes a repulsion V between dimers facing each other, see Eq. ͑1͒. For tЈ / t = 0, this model reduces to the standard QDM, which has an RVB phase for 0.75Շ V / t Շ 1. 4,5 To map out the complete phase diagram of this model, we have done the same analysis of the preceding paragraph for different values of V / t between 0 and 1. As an example, we show in Fig. 6 the finite-size scaling of the topological gap for tЈ / t = 1 and several values of V / t. The global behavior is the same as before and we have clear evidence that the topological gap present at V / t = 0 persists up to V / t ϳ 0.25, and that it opens again for V / t ϳ 0.8. In that case, as in many other cases, it turned out to be possible to actually observe the opening of the gap upon leaving again the RVB phase before the convergence problems were too strong.

The resulting phase diagram is depicted in Fig. 7, where different symbols have been used depending on whether the boundary was determined from the increase of the correlation length with the size or from the loss of convergence: Solid circles when the closing of the topological gap was observable and solid diamonds when the convergence of the GFMC was lost for large systems. Interestingly enough, these different criteria build a relatively smooth line, a good indication that it can be interpreted as a phase boundary. Note that we did not perform simulations for V / t Ͼ 1 where, for a vanishing tЈ / t, a crystalline phase with staggered dimer order is stabilized. So we cannot exclude that the boundary extends beyond V / t Ͼ 1 for small tЈ / t. Remarkably, the RVB phase we found for V / t = 0 is connected to the RVB phase reported before for the standard QDM, i.e., tЈ / t = 0, the total RVB phase building up a large stripe that encompasses a significant portion of the phase diagram. We have also calculated static correlation function for several values of the parameter, but they merely confirm the identification of the phases and are not reported for brevity.

V. CONCLUSIONS

Coming back to our long-term motivation, namely to find an RVB phase in a realistic model of Mott insulators, this paper contains significant results of two sorts. First of all, we have shown ͑see the Appendix͒ that, starting from the quasidegenerate mean-field ground state of a Kugel-Khomskii spin-orbital model, one can construct a QDM with two remarkable properties: It describes a competition between two kinetic terms of comparable magnitude, and all off-diagonal matrix elements in the dimer basis are negative. This has allowed us to implement the GFMC and to investigate the results of the competition between these terms. It turns out that the competition between these kinetic terms leads to the disappearance of the ͱ 12ϫ ͱ 12 crystalline order when tЈ / t ϳ 1.6. This transition is similar to the transition into the RVB phase that happens in the standard QDM upon approaching the Rokhsar-Kivelson point. Indeed, the two phases can be connected into a single RVB phase in the context of a generalized QDM. As far as the numerical investigation of the model is concerned, the main open issue is to pin down the nature of the phase that occurs when the six-dimer kinetic term dominates. Unfortunately, the GFMC suffers from severe statistical fluctuations whenever the guiding function is not accurate, i.e., for clusters larger that 108 sites and large tЈ / t. Therefore, we cannot make any definite statements on the phase where the six-site term dominates. Another interesting question is of course the nature of the quantum phase transitions between these phases ͑continuous or first order͒. We are currently working on that rather subtle issue in the context of the standard QDM.

The general features of the RVB phase are consistent with the phenomenology of LiNiO 2 , which exhibits neither orbital nor magnetic long-range order. A number of points deserve further investigation however. The precise form of the QDM does not seem to be an issue: The actual model that can be derived along the Rokhsar-Kivelson lines has more terms ͑see the Appendix͒, but preliminary results show that the RVB phase is present in that model as well. The fact that the RVB phase does not contain the point tЈ / t = 1.34 derived in the Appendix is not really an issue either. First of all, this ratio was determined for vanishing Hund's rule coupling and one vanishing hopping integral, and its precise value in that case should at best be taken as an indication of its order of magnitude in the actual system. Besides, the other six-site terms pull the RVB region down to smaller values of the relative ratio of the six-dimer term to the four-dimer term. What would deserve more attention is the validity of the expansion that leads to the effective QDM. The small parameter of the expansion is not that small ͑␣ =1/ ͱ 2, see the Appendix͒, and it would be very useful to better understand to which extent such an expansion can be controlled. Nevertheless, the present results strongly suggest that the presence of an RVB liquid phase between competing ordered phases is a generic feature of QDM, and that to identify such a phase in a realistic Mott insulator via an effective QDM is very promising.

͉⌽ ij ͘ = ͉ ij ͘ ͉ ij ͘, ͑A4͒
where ͉ ij ͘ and ͉ ij ͘ are, respectively, the spin and orbital components. The spin component is the usual singlet given by

͉ ij ͘ = ␣͉͑↑ i ↓ j ͘ -͉↓ i ↑ j ͒͘, ͑A5͒
regardless of the orientation of the bond. We have denoted by ␣ the normalization coefficient, whose explicit value is of course given by 1 / ͱ 2, to be able to keep track of the order in ␣ of various overlaps. The orbital part depends on the bond and is given by

͉ 12 ͘ = ͉a 1 ͉͘a 2 ͘, ͑A6͒ ͉ 13 ͘ = ͯ-1 2 a 1 - ͱ 3 2 b 1 ʹͯ -1 2 a 3 - ͱ 3 2 b 3 ʹ , ͑A7͒ ͉ 23 ͘ = ͯ-1 2 a 2 + ͱ 3 2 b 2 ʹͯ -1 2 a 3 + ͱ 3 2 b 3 ʹ , ͑A8͒
with the convention of Fig. 8, and with ͉a͘ = ͉d 3z 2 -r 2͘ and ͉b͘ = ͉d x 2 -y 2͘. We also use the convention that the wave function has a plus sign if the two sites are in the order defined by the arrows of Fig. 8, and a minus sign otherwise. With these definitions, the matrix element of the Hamiltonian on each bond is given by FIG. 8. Sign convention for a triangle.

͗⌽ ij ͉H ij ͉⌽ ij ͘ = ⑀2␣ 2 = ⑀͗⌽ ij ͉⌽ ij ͘, ͑A9͒
where ⑀ is defined by

⑀ = -4U ˜th 2 U ˜2 -J 2 Ͻ 0. ͑A10͒
According to the mean-field analysis of Ref. 11, the wave functions obtained as tensor products of these wave functions on all dimer coverings on the triangular lattice should constitute a good variational basis if t h ЈӶ t h and J Ӷ U ˜. In the following, we consider for simplicity the case t h Ј= 0 and J =0.

͑A12͒

The inverse of the square root of the overlap matrix cannot be calculated exactly, but this can be done approximately in the context of an expansion in powers of ␣. Indeed, the overlap matrix can be expanded as S = I + 2A␣ 4 + 2B␣ 6 + O͑␣ 8 ͒, ͑A13͒

which leads to ͑S͒ -1/2 = I -A␣ 4 -B␣ 6 + O͑␣ 8 ͒. ͑A14͒

In these expressions, the matrices A and B only have nonvanishing matrix elements, equal to 1, between configurations that are the same except on a four-dimer loop ͑see Fig. 9͒ or on one of the three types of six-dimer loops ͑see Fig. 10͒ where the matrices C and D have nonvanishing matrix elements under the same conditions as matrices A and B. Since all these expansions start with ␣ 4 , it is clear that the first contribution to the diagonal part of H eff will be of order ␣ 8 . So, to order ␣ 6 , the effective Hamiltonian will only have off-diagonal matrix elements. Moreover, to this order, these matrix elements are simply given by ͗m͉H ˜͉n͘ = ͗⌿ m ͉H ˜͉⌿ n ͘. They are tabulated in Table I for configurations that are the same except on a four-dimer loop ͑Fig. 9͒ or on one of the three types of six-dimer loops ͑Fig. 10͒. All these matrix elements are negative. Among the six-dimer loop terms, the matrix element of type ͑3͒ is the largest, and its ratio to the four-dimer loop term is 1.34.

FIG. 1 .

 1 FIG. 1. Example of type-A ͑left͒ and type-B ͑right͒ clusters with 16 and 12 sites, respectively.

FIG. 4 .

 4 FIG. 4. ͑Color online͒ Topological gaps for V / t = 0 as a function of 1 / ͱ N, where N is the number of sites and for different values of tЈ / t. Upper panel: Small values of tЈ / t, where the gap opens for large clusters. Lower panel: Larger values of tЈ / t. For tЈ / t տ 1.6, the finite-size gap closes upon increasing the cluster size, signaling a liquid phase.

FIG. 6 .

 6 FIG.6. ͑Color online͒ Topological gap for tЈ / t = 1 and for various V / t as a function of 1 / ͱ N, where N is the number of sites.Small and large values of V / t have been shown separately in the upper and lower panels for clarity.

FIG. 9 .

 9 FIG. 9. A four-dimer loop.

Following

  Rokhsar and Kivelson,1 the idea is now to perform a unitary transformation to derive the effective QDM. If one defines the overlap matrix by S mn = ͗⌿ m ͉ ⌿ n ͘, where ͉⌿ m ͘ and ͉⌿ n ͘ are dimer coverings, the states defined by͉m͘ = ͚ n ͑S -1/2 ͒ m,n ͉⌿ n ͘ ͑ A11͒constitute an orthonormal basis and the matrix elements of the Hamiltonian in this basis are given by H mn eff ϵ ͗m͉H͉n͘ = ͚ kl ͑S -1/2 ͒ mk ͗⌿ k ͉H͉⌿ l ͑͘S -1/2 ͒ ln .

  , respectively. Similarly, the Hamiltonian matrix H ˜defined by H ˜mn = ͗⌿ m ͉H͉⌿ n ͘ -⑀N d ␦ mn , where N d is the number of dimers, has an expansion in powers of ␣ that reads H ˜= C␣ 4 + D␣ 6 + O͑␣ 8 ͒, ͑A15͒

TABLE I .

 I Overlap and matrix element for two configurations differing by a loop-4, a loop-6 of type ͑1͒-͑3͒.

		loop-4	loop-6 ͑1͒	loop-6 ͑2͒	loop-6 ͑3͒
	͗m͉H͉n͘	-	66⑀␣ 4 16 2	-	103⑀␣ 6 16 2	-	24⑀␣ 6 16 2	-	177⑀␣ 6 16 2
		Ӎ-0.0644⑀	Ӎ-0.0503⑀	Ӎ-0.0012⑀	Ӎ-0.0864⑀

ACKNOWLEDGMENTS

The authors acknowledge useful discussions with P. Fazekas, M. Ferrero, and K. Penc. This work was supported by the Swiss National Fund and by MaNEP. One of the authors ͑F.B.͒ is supported by CNR-INFM and MIUR ͑COFIN 2005͒.

APPENDIX: DERIVATION OF THE MODEL

LiNiO 2 is a layered compound in which the Ni 3+ ions are in a low spin S = 1 2 state with a twofold orbital degeneracy. A fairly general description of this system is given by a Kugel-Khomskii Hamiltonian defined in terms of two hopping integrals t h and t h Ј, the on-site Coulomb repulsion U ˜and the Hund's coupling J which, on a given bond, takes the form 11

with the usual definitions for the projectors on the singlet and triplet states of a pair of spins,

The vector n ij z depends on the type of bond. With the convention of Fig. 8, they are given by n 12 z = ͑0,0,1͒,

͑A3͒

The operators T i are pseudospin operators acting on the orbitals.

On a given bond ͑see Fig. 8͒, a dimer is defined by the following wave function: