Quantum spin liquid and cluster Mott insulator phases in the Mo3O8 magnets

We unveil the microscopic origin of largely debated magnetism in the Mo3O8 quantum systems. Upon considering an extended Hubbard model at 1/6 filling on the anisotropic kagomé lattice formed by the Mo atoms, we argue that its ground state is determined by the competition between kinetic energy and intersite Coulomb interactions, which is controlled by the trimerisation of the kagomé lattice into the Mo3O13 clusters, and the sign of hopping parameters, specifying the electron localisation at such clusters. Based on first-principles calculations, we show that the strong interaction limit reveals a plaquette charge order with unpaired spins at the resonating hexagons that can be realised in LiZn2Mo3O8, and whose origin is solely related to the opposite signs of intracluster and intercluster hoppings, in contrast to all previous scenarios. On the other hand, both Li2InMo3O8 and Li2ScMo3O8 are demonstrated to fall into the weak interaction limit where the electrons are well localised at the Mo3O13 clusters. While the former is found to exhibit long-range antiferromagnetic order, the latter is more likely to reveal short-range order with quantum spin liquid-like excitations. Our results not only reproduce most of the experimentally observed features of the Mo3O8 systems, but will also help to describe various properties in other quantum cluster magnets.

During the past few years, the Mo 3 O 8 cluster magnets have attracted a great deal of both experimental and theoretical attention as a candidate to host QSL. In these compounds, the Mo atoms arranged in anisotropic kagomé layers are trimerised, and the [Mo 3 O 13 ] 15− clusters form a triangular lattice, as shown in Fig. 1a 16,17 . Six out of seven valence electrons in the Mo 3 O 13 cluster are responsible for a strong intracluster metal-metal bonding, and the seventh electron remains unpaired occupying a totally symmetric molecular a 1 state. LiZn 2 Mo 3 O 8 was first reported to exhibit a QSL behaviour [18][19][20] . The magnetic susceptibility of LiZn 2 Mo 3 O 8 has been experimentally shown to follow a Curie-Weiss law with low-and hightemperature regimes transitioning at T C~9 6 K, whose Curie constants are related as C L ≈ C H /3 and where the disappearance of 2/3 of the paramagnetic spins was attributed to valence bond condensation on the triangular lattice of the Mo 3 O 13 clusters. In a first attempt to explain these unusual features, the authors of ref. 21 suggested the formation of an emergent honeycomb lattice due to opposite rotations of the Mo 3 O 13 clusters effectively decoupling the central cluster with an orphan paramagnetic spin. Another scenario was outlined in ref. 22 , where a plaquette charge order (PCO) existing in a Mott insulator on the anisotropic kagomé lattice at 1/6 filling was conjectured to host a U(1) QSL state with the spinon Fermi surface that is reconstructed at low temperatures filling 2/3 of the spinon states.
However, the adequacy of the proposed mechanisms was questioned with recently synthesised Li 2 InMo 3 O 8 and Li 2 ScMo 3 O 8 , both featuring magnetic moments well localised at the Mo 3 O 13 clusters. While the former was identified with a Neél 120 ∘ magnetic order at T N = 12 K, for the latter no magnetic ordering has been observed down to 0.5 K 23,24 . Instead, muon spin rotation and inelastic neutron scattering measurements suggested that Li 2 ScMo 3 O 8 shows a short-range magnetic order below 4 K with QSL-like excitations. Despite having similar crystal structures, these Mo 3 O 8 systems manifest essentially unalike magnetic properties, whose enigmatic origin remains an unsolved problem.
Abstract electronic and spin models abounding with fascinating theories have been indispensable in advancing our understanding of QSLs. However, due to the complexity of the problem, firstprinciples calculations, which were proven to be an extremely valuable theoretical tool in a broad range of subfields in materials science, stay in the shadow of many hypothetical proposals that are often disconnected from real systems, seriously hampering the progress. In this work, upon revising a single-orbital extended Hubbard model on the anisotropic kagomé lattice at 1/6 filling, we uncover two different regimes of the plaquette charge ordered and cluster Mott insulator states governed by the interplay of kinetic energy and intersite Coulomb interactions, that were overlooked in all previous studies. By means of first-principles calculations, we will determine the material-specific model parameters for the family of the Mo 3 O 8 quantum magnets and demonstrate that these states can indeed be realised in LiZn 2 Mo 3 O 8 , Li 2 ScMo 3 O 8 and Li 2 InMo 3 O 8 . In particular, we will show that the opposite sign of hopping parameters is peculiar to the trimerised kagomé layers of the Mo 3 O 8 systems, while their strength specifies the degree of electron localisation at the Mo 3 O 13 clusters. Fig. 1b is the extended Hubbard model on the kagomé lattice at 1/6 filling:

The model of interest shown in
where c yσ m (c σ m ) creates (annihilates) an electron with spin σ at site m, n σ m ¼ c yσ m c σ m is the density operator (n m ¼ n " m þ n # m ), t and t 0 stand for intracluster and intercluster hopping parameters defined on the T and T 0 triangles, respectively, and the interaction terms include the on-site U, intracluster V and intercluster V 0 Coulomb repulsions (see Supplementary Note 1). Taking a shorter bond length in T , the T triangles correspond to the Mo 3 O 13 clusters accommodating one unpaired electron and forming a triangular lattice.

First-principles calculations
The calculated band structures of LiZn 2 Mo 3 O 8 , Li 2 ScMo 3 O 8 and Li 2 InMo 3 O 8 are shown in Fig. 2a, indicating the a 2 and e 2 states below the Fermi level, which are responsible for the Mo-Mo bonding in the Mo 3 O 13 cluster, and the nonbonding a 1 and e 1 states well separated from the rest of the spectrum. According to our results, the splitting Δ of the a 1 and e 1 levels varies significantly within the systems, being relatively small in LiZn 2 Mo 3 O 8 and gradually increasing in Li 2 ScMo 3 O 8 and Li 2 In-Mo 3 O 8 . Given that one electron occupies these states, the value of Δ suggests the degree of electron localisation at the molecular a 1 level, which is, in turn, responsible for the formation of localised magnetic moments. As will be shown below, this fact is one of the main ingredients to understand distinct magnetic behaviour in these materials. The extended Hubbard model, Eq. (1), for the systems in question was constructed in the basis of Wannier functions corresponding to the a 1 and e 1 bands, as shown in Fig. 2b. The full set of model parameters is given in Table 1 (see Supplementary Note 4). Due to the trimerisation of the kagomé lattice, V 0 <V ( U and jt 0 j < jtj hold for each system. But importantly t and t 0 always have opposite signs, and t < 0 while t 0 > 0. This can be related to the fact that in the Mo 3 O 13 cluster the direct d-d (always negative) hopping dominates due to shorter Mo-Mo bonds. Because this term decays rapidly with the metal-metal distance (~1/r 5 , ref. 25 ), the hopping process via common oxygens having the opposite sign starts to prevail between the clusters, and t 0 turns out to be positive. In a more general way, the sign alternation of hopping parameters can be regarded to occur due to the formation of the Mo 3 O 13 clusters, thus enforcing the localisation of a single unpaired electron at the molecular a 1 level.
The competition between intracluster and intercluster interactions along with the electron filling and geometrical frustration brings about different regimes of electron localisation. Taking into account the calculated values of t=V 0 and t 0 =V, we will address two limits of Eq. (1).

Plaquette charge order
Let us consider jtj ( V 0 and jtj 0 ( V. Due to 1/6 filling, the Hubbard U cannot localise electrons at the lattice sites, and as a result they move without encountering any double occupancy. Since U is not operative, it is the intersite V and V 0 that are responsible for electron localisation leading to a highly degenerate charge ordered state, where each corner-sharing triangle hosts exactly one electron. This degeneracy is further lifted by hopping parameters that induce collective tunnelling processes, when the electrons hop either clockwise or counter-clockwise along the T and T 0 bonds stabilising a charge pattern with three electrons at the hexagons, as shown in Fig. 3a, b. To lowest order in t=V 0 and t 0 =V, it corresponds to the quantum dimer model for two plaquette states where the sum runs over all hexagons 26,27 . When mapped onto the dual hexagonal lattice, the ground state of H ⎔ for spinless electrons is described by the PCO (Fig. 3c) with an emergent triangular lattice of resonating hexagons 28,29 . Calculations within exact diagonalization for finite geometries show that the PCO is the ground state in the case of spinful electrons (see Supplementary Note 2), that will be regarded as the strong interaction limit of Eq. (1).
One can further include antiferromagnetic spin fluctuations between next-nearest neighbours in each hexagon, and t nn is the corresponding hopping (Fig. 3a). Assuming that the PCO effectively decouples hexagons, H D ¼ H ⎔ þ H S for a single hexagon can be solved exactly yielding four fourfold degenerate states (see Supplementary Note 2). When g 1 and g 2 have opposite signs (g 1 > 0 and g 2 < 0), regardless of the value of J nn the ground state of H D is represented by: withg ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , ψ 3 j i and ψ 4 j i can be obtained by applying time-reversal symmetry to the above states. As seen from Eq. (2) and schematically shown in Fig. 4a, the resulting ground state of the resonating hexagon is given as a superposition of the plaquette states A j i and B j i, each having valence bonds resonating in the hexagon and leaving one spin unpaired, whose location is smeared out in the hexagon owing to the resonating nature of valence bonds. Such an unusual entanglement with dangling spins originates solely from the asymmetry of tunnelling processes that, in turn, facilitates partial singlet pairing between the resonating electrons, while the unpaired spins behave paramagnetically in a thermodynamic limit.
From our first-principles calculations, a small splitting Δ of the a 1 and e 1 levels is found in LiZn 2 Mo 3 O 8 preventing the electrons from being localised at the molecular states, and the calculated t=V 0 and t 0 =V suggest the PCO with unpaired spins at the resonating hexagons as a possible ground state for LiZn 2 Mo 3 O 8 .
Given g 1 = 13.5 meV, g 2 = −40.1 meV and J nn = 1.4 meV, the calculated T C~9 2.0 K between two paramagnetic regimes is in excellent agreement with neutron powder diffraction data 18 (see Supplementary Fig. 5).
Originally, the PCO state on the kagomé lattice in the strong interaction limit was studied for spinless fermions at 1/3 filling 27,30 . For spinful electrons, ferromagnetism has been proposed in ref. 31 as one of the ground states at 1/6 and 5/6 fillings. This scenario turns out to be justified only when both t and t 0 have the same sign, so that the Hamiltonian of the quantum dimer model H ⎔ can be chosen to contain exclusively positive entries owing to the bipartite lattice, and the Perron-Frobenius theorem can be used for a finite geometry. When the hopping parameters have opposite signs, the Perron-Frobenius theorem is no longer applicable, as is also confirmed by our calculations for a single cluster and finite geometries with periodic boundary conditions, where antisymmetric configurations with minimal total spin are found to be more favourable (see Supplementary Fig. 3). Several studies have also considered the same model for spinful electrons at 1/3 filling with hopping parameters of the same sign 32,33 , where the PCO was predicted in the limit of strong intersite Coulomb interactions. However, no possible scenario for the two paramagnetic regimes has been proposed.
The PCO state was reconsidered in ref. 22 , where the disappearance of 2/3 of the spins was attributed to a conjectural U(1) QSL that reconstructs the spinon Fermi surface at low temperatures partially filling the spinon states, inert to external magnetic field. A striking difference with the present study is that the hopping parameters in ref. 22 were assumed to have the same sign. In this case, one can find that the ground state for g 1 < 0 and g 2 < 0 on a single hexagon is given by a symmetric configuration of the plaquette states, ð A j i þ B j iÞ= ffiffi ffi 2 p (see Supplementary Note 2), and the PCO itself does not feature any pairing of spins with unpaired electrons. In fact, the emergence of the U(1) QSL suggested in ref. 22 was introduced phenomenologically by "artificially" modifying the ring tunnelling processes to mimic the PCO-driven breaking of translational symmetry. In contrast, in the present case with g 1 < 0 and g 2 > 0 the PCO state is represented by antisymmetric configurations of A j i and B j i, where each plaquette is, in turn, given as an antisymmetric configuration of spins forming valence bonds with one dangling spin in a hexagon. Finally, the unusual pairing given by Eq. (2) can also be realised when g 1 and g 2 have the same sign (g 1 < 0 and g 2 < 0) and the antiferromagnetic coupling between next-nearest neighbours is large, J nn > 2 3 ðÀg 1 À g 2 ÀgÞ, as was suggested in refs. 22,34 . Although the PCO can be destroyed by strong spin fluctuations 32 , the calculated thermodynamic properties for a single hexagon shown in Fig. 4b clearly demonstrate that two paramagnetic regimes possess a much higher T C when g 1 and g 2 have opposite signs. Moreover, the negligibly small J nn obtained for LiZn 2 Mo 3 O 8 shows that the partial spin pairing in the resonating hexagons is driven solely by g 1 > 0 and g 2 < 0. Thus, the PCO state presented in our study arises when t and t 0 have opposite signs, and the partial pairing of spins in resonating hexagons comes about in a straightforward way from the asymmetry of tunnelling processes.
Signatures of the partial spin pairing driven by the PCO state can be directly traced from the spin-spin correlation functions (see Supplementary Fig. 4 for the density-density and plaquette-plaquette correlation functions). As visualised in Fig. 5a, the correlation function for the PCO state with g 1 < 0 and g 2 < 0 is essentially positive, indicating no spin pairing, and the calculated structure factor has maxima at the Γ and K points of the extended Brillouin zone. On the contrary, the correlation function for the PCO state with g 1 < 0 and g 2 > 0 shown in Fig. 5b changes the sign between next-nearest neighbours, and the structure factor reveals maxima at the K points of the original Brillouin zone. This data can further be used in the analysis of inelastic neutron scattering experiments to verify our findings.

Cluster Hubbard model
As |t| increases, the electrons start moving freely within the T triangle, and the number of electrons at the adjacent T 0 triangles fluctuates. When jtj $ V 0 , the perturbation theory considered above breaks down, and the electrons minimise their energy by forming bound "molecular" states. As a result, the original model in Eq. (1) can be reformulated as a three-orbital extended Hubbard model on the triangular lattice formed by the T triangles: with Δ = 3t. As follows, H CF has the form of crystal field that splits the electronic states at the T triangle into the single a 1 and double degenerate e 1 states with energy levels 2Δ 3 and À Δ 3 , respectively: a 1 j i ¼ 1 ffiffi j iÞ with w = e 2πi/3 . Importantly, the a 1 state is occupied when Δ < 0 (t < 0). Despite the weak interaction limit, the electrons are localised at the T triangles due to the large splitting of molecular levels. We refer to this state as a cluster Mott insulator as opposed to the PCO phase, where the electron delocalisation, in turn, is entirely driven by intersite Coulomb interactions. When both t < 0 and t 0 > 0 are large, the electrons localised at the T triangles can develop long-range magnetic order. In this limit, the on-site e U ¼ Uþ2V 3 comes back into play and forbids any double occupancy at the T triangles, and the corresponding spin model H 4 ¼ P hiji J 4 S i Á S j on the triangular lattice can be derived to second order in t 0 =U and t 0 =Δ (see Supplementary Note 3): which can be both ferro-and antiferromagnetic, that can explain why some of the recently found Mo 3 O 8 systems are ferromagnetic insulators 35 Fig. 6. Consequently, although the electrons tend to localise at the Mo 3 O 13 clusters, Li 2 ScMo 3 O 8 is more likely to fall into the intermediate regime, where any long-range magnetic order is suppressed by quantum fluctuations down to low temperatures. Since the number of electrons at the T 0 triangles is allowed to fluctuate when t=V 0 and t 0 =V are not strong, we conclude that the magnetic order in Li 2 ScMo 3 O 8 is short range with possible QSL-like excitations.

DISCUSSION
Having considered an extended Hubbard model on the anisotropic kagomé lattice at 1/6 filling as the low-energy model for the Mo 3 O 8 cluster magnets, we showed that it features two different limits: a PCO of resonating hexagons with valence bond condensation and orphan spins, as realised in quantum paramagnet LiZn 2 Mo 3 O 8 , and a cluster Mott insulator with the electrons localised at the kagomé triangles, as revealed in Li 2 InMo 3 O 8 and Li 2 ScMo 3 O 8 showing a Néel-type antiferromagnetic order and QSL behaviour, respectively. Based on firstprinciples calculations, we demonstrated that their manifestation can be attributed to the trimerisation of the kagomé lattice controlling the competition between kinetic energy and intersite Coulomb interactions and specifying the character of electron localisation at the Mo 3 O 13 clusters, that can help to unravel a largely speculated origin of magnetism in these systems.
We believe that the opposite sign of t and t 0 is a peculiar aspect of the trimerised kagomé lattice inherent to the Mo 3 O 8 quantum magnets that reflects the bonding character of the Mo 3 O 13  clusters and plays an important role in their different magnetic behaviour. According to the general Jahn-Teller theorem, a lattice trimerisation should lift the ground-state degeneracy so that a single electron resides at the a 1 orbital of the T triangle forming a one-dimensional representation of the point group, that only occurs when t < 0 and t 0 > 0.
The key point of our study is that the magnetism of the Mo 3 O 8 systems can be described within a single model, where both strong localisation and itineracy appear as the opposing limits of competing interactions. Our calculations provide a good agreement with experimental data for both Li 2 InMo 3 O 8 and Li 2 Sc-Mo 3 O 8 23,24 and can explain the difference in their magnetic properties. The results obtained for LiZn 2 Mo 3 O 8 also agree with inelastic neutron scattering data for powder samples 19,20 , reporting a gapless spectrum of magnetic excitations and short-range spin correlations, which are inherent to the PCO state, as well as the absence of detectable sharp excitations that suggests a disordered valence bond solid or a resonating valence bond state. Available momentum dependencies of the inelastic magnetic scattering intensity reveal scatterings concentrated at small momenta q < 1.0 Å −1 with the maximum at q~0.41 Å −1 at 1.7 K 19 , which is within the range of the first Brillouin zone (the length of the in-plane reciprocal lattice vector is~1.25 Å −1 for the experimental structure of LiZn 2 Mo 3 O 8 36 ) and is also consistent with our calculations of the static spin structure factor indicating the maximum weight at the K point of the first Brillouin zone, as shown in Fig. 5c (~0.72 Å −1 for the experimental structure of LiZn 2 Mo 3 O 8 ). That being said, the possible origin of paramagnetism in LiZn 2 Mo 3 O 8 can still be speculated.
Our first-principles calculations would possibly rule out some previous scenarios proposed within the same model for decoupling 2/3 of the spins in LiZn 2 Mo 3 O 8 at low temperatures, such as a U(1) QSL state with the spinon Fermi surface 22 and strong coupling between next-nearest neighbours 34,37 . Cooperative rotations of the Mo 3 O 13 clusters resulting in an emergent hexagonal lattice 21 assume electrons to be well localised at the clusters and are at variance with our results in the strong interaction limit. The presence of mobile Li ions at~50 − 100 K can also locally contribute to the electron localisation [18][19][20] . However, the peaks from the Zn/Li sites observed in Li NMR measurements have different maxima below the paramagnetic phase transition, and the Li NMR data from ref. 20 show a minimal spin concentration from impurities, indicating that the origin of magnetic response in LiZn 2 Mo 3 O 8 arises from the Mo layers. Nevertheless, the dynamical and disorder effects cannot be unequivocally excluded from the consideration and should be carefully addressed in further studies. Importantly, the low-energy models derived in this study suggest that the delocalised nature of electrons in LiZn 2 Mo 3 O 8 plays a very intricate role in defining magnetic properties advancing short-range spin dynamics over long-range ordering 38 . We believe that further experimental studies on single crystals are important to verify the origin of unusual magnetic behaviour in LiZn 2 Mo 3 O 8 , as was previously done to identify fractional spinon excitations of the quantum antiferromagnetic chains 39 and continuous spin excitations of the triangular-lattice QSL 40 .
Finally, it is known that spin-1 2 systems with an odd number of electrons can reveal both long-range order and short-range correlations with topological excitations 41 , and various scenarios with unusual magnetic properties at low temperature are expected to be realised in other cluster magnets. For example, Na 3 A 2 (MoO 4 ) 2 Mo 3 O 8 (A = In, Sc) was experimentally shown to feature a paramagnetic behaviour similar to LiZn 2 Mo 3 O 8 42 , or Nb 3 Cl 8 with the same number of electrons and similar crystal structure was found to be completely nonmagnetic below 90 K, in contrast to the Mo 3 O 8 systems 43,44 . The analysis presented in this work can be generally applied to these and other trimerised quantum systems, such as Li 2 In 1−x Sc x Mo 3 O 8 37 and 8 45 , and further experimental and theoretical studies of these materials hold promise for discovering novel quantum effects.

First-principles calculations
Electronic structure calculations were performed within local density approximation 46 and projected augmented waves 47 , as implemented in the Vienna ab initio simulation package (VASP) 48 . Band structure calculations including spin-orbit coupling were carried out by using norm-conserving pseudopotentials as implemented in the Quantum ESPRESSO package (QE) 49 . The plane wave cutoff was set to 600 and 1360 eV for VASP and QE, respectively, the Brillouin zone was sampled by a 10 × 10 × 5 Monkhorst-Pack k-point mesh 50 , and the convergence criteria for the total energy calculations was 10 −9 eV. High-symmetry k-points used for the band structure calculations are K ¼ ð 1 3 ; 1 3 ; 0Þ, M ¼ ð 1 2 ; 0; 0Þ, A ¼ ð0; 0; 1 2 Þ, H ¼ ð 1 3 ; 1 3 ; 1 2 Þ, L ¼ ð 1 2 ; 0; 1 2 Þ. To calculate model parameters of Eq. (1), we employed the maximally localised Wannier functions as implemented in the wannier90 package 51,52 . The parameters of Coulomb interactions were calculated within constrained random-phase approximation in the basis of Wannier functions 53,54 .
The crystal structure parameters of Li 2 InMo 3 O 8 , Li 2 ScMo 3 O 8 and LiZn 2 Mo 3 O 8 adopted from refs. 23,36 are given in Supplementary Tables 1  and 2. For LiZn 2 Mo 3 O 8 , the disorder of the Li and Zn atoms was fixed, and their arrangement was checked to give a small effect on the final set of model parameters. It is worth noting that there is a small splitting of the a 1 and e 1 states corresponding to the inequivalent Mo layers due to the Li/Zn disorder. Both octahedral sites having a higher Li occupancy 20 and tetrahedral sites can contribute to the local splitting of molecular levels leading to a small change of the model parameters.
The crystal structures were visualised with VESTA 55 .

Exact diagonalisation
Exact diagonalisation of the quantum dimer model H D was performed on finite clusters containing N = 27 and 36 sites with periodic boundary conditions. The spin-spin correlation function is computed as: C s ðr 0 ; r k Þ ¼ hS 0 Á S k i À hS 0 i Á hS k i; where 〈...〉 is the thermal average. The corresponding structure factor is calculated on the extended Brillouin zone (double of the original one): e ÀiqÁðrk Àr0Þ C s ðr 0 ; r k Þ: Specific heat and spin susceptibility for H D on a single hexagon are calculated as: respectively, where β is inverse temperature, and m z is the total magnetic moment of the plaquette state.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.