Clusterization transition between cluster Mott insulators on a breathing Kagom\'{e} lattice

Motivated by recent experimental progress on various cluster Mott insulators, we study an extended Hubbard model on a breathing Kagom\'{e} lattice with a single electron orbital and $1/6$ electron filling. Two distinct types of cluster localization are found in the cluster Mott regime due to the presence of the electron repulsion between neighboring sites, rather than from the on-site Hubbard interaction in the conventional Mott insulators. We introduce a unified parton construction framework to accommodate both type of cluster Mott insulating phase as well as a trivial Ferm liquid metal and discuss the phase transitions in the phase diagram. It is shown that, in one of the cluster localization phases, the strong inter-site repulsion results into locally metallic behavior within one of two triangular clusters on the breathing Kagom\'{e} lattice. We further comment on experimental relevance to existing Mo-based cluster magnets.


I. INTRODUCTION
Cluster Mott insulators (CMIs) seem to become a new frontier for exploring the emergent correlated physics [1][2][3]. The cluster magnet 1T-TaS 2 develops a commensurate charge density wave order at about 120K, and the enlarged unit cell due to the charge order then has a Davidstar shape with 13 lattice sites [4,5]. The enlarged unit cell traps one unpaired electron that is Mott localized on the cluster unit of the David star, and the system in the commensurate charge density wave state forms a CMI in two spatial dimensions. The surging field of twistronics, that was initiated from the twisted bilayer graphenes [6][7][8][9][10][11], potentially can be another example of realizing CMIs where the electrons are localized on the large moiré unit cell. These moiré unit cells are often one to two orders of magnitude larger than the lattice constant of the original untwisted crystals. The twisting procedure provides a new knob to tune the physical properties of the underlying systems. The common ingredient shared by these systems is the large cluster unit for the electronic degrees of freedom, and the longer range interactions ought to be considered [1][2][3][12][13][14][15][16][17][18][19][20][21]. This ingredient leads to distinct and interesting features and experimental consequences in different realizations of cluster localization.
In expectation of the large internal electronic degrees of freedom inside each cluster, we explore the rich phase diagram of the CMIs. We present the observation by * These authors contributed equally.
showing the existence of distinct cluster localizations and study the phase transition between the CMIs on the breathing Kagomé lattice. This is partly motivated by the experiments on various Mo-based two-dimensional cluster magnets [22][23][24][25][26]. We study a 1/6-filled extended Hubbard model with the nearest-neighbor repulsions on a breathing Kagomé lattice. The Mott insulating physics in this partially filled system arises from the large nearest-neighbor repulsions [1,12,27] and localization of the electrons in the triangular cluster units. Due to the asymmetry between the up-and down-triangles and the resulting difference in the interactions and hoppings, two different cluster localizations are expected. We first show that, for the case of cluster localization on only one type of triangular units (e.g. the up ones) in the strong breath- arXiv:2011.02813v1 [cond-mat.str-el] 5 Nov 2020 ing limit, the ground state is smoothly connected to the one for the triangular lattice Hubbard model at half filling. We then explore the phase transition between two distinct cluster localizations and further address the consequences on the spin physics. In terms of the lattice gauge theory formulation, this transition is identified as a Higgs transition. The correspondence between the lattice gauge theory formulation and the physical variables are clarified.
The remainder of the paper is organized as follows. We begin in Sec. II by introducing the extended Hubbard model on the breathing Kagomé lattice with a single electron orbital and 1/6 electron filling. In Sec III we explore the strong breathing limit and discuss the type-I CMI and the Mott transition. We go beyond the strong breathing limit in Sec. IV and study the type-II CMI as well as its emergent U(1) c gauge structure. We further introduce a unified parton construction in Sec. V to reveal the rich physics of both type-I and type-II CMIs at the mean-field level. We establish the generic phase diagram and discuss the phase transition between two distinct CMIs. We conclude by discussing the experimental relevance and consequence about Mo-based cluster magnets in Sec. VI.

II. EXTENDED HUBBARD MODEL
We start from the extended Hubbard model on the breathing Kagomé lattice, where c † iσ (c iσ ) creates (annihilates) an electron with spin σ at the lattice site i, and t ij = t 1 (t 2 ), V ij = V 1 (V 2 ) for ij on the up-(down-) triangles. Here n i ≡ σ c † iσ c iσ defines the electron occupation number at the lattice site i. We are interested in the regime with one electron in each Kagomé lattice unit cell, and thus the electron filling for this Hubbard model is 1/6 [22,24]. This model is expected to capture the essential physics of the Mobased cluster magnets [2,12]. For the fractional 1/6 filling here, the Mott localization is driven by the inter-site repulsions (V 1 , V 2 ) rather than the on-site Hubbard-U interaction and the electrons are localized in the (elementary) triangles of the Kagomé lattice instead of the lattice sites. Due to the asymmetry between the up-and downtriangles, the Mott localization in the up-triangles and down-triangles does not occur simultaneously. Setting the Hubbard-U as the largest energy scale, we study the properties of the model and explore its phase diagram.
The kinetic part of the model can be readily diagonalized and the electrons form the following three bands, where k 1 ≡ k · b 1 , k 2 ≡ k · b 2 , and b 1 , b 2 are two elementary lattice vectors of the underlying Bravais lattice (see Fig. 1). These three electron bands are wellseparated from each other and only touch at certain discrete momentum points. In particular, E 1 (k) and E 2 (k) have Dirac-point band touchings at the Brillouin zone corners when t 1 = t 2 . With the 1/6 electron filling, the electrons fill half the lowest band E 1 (k) and the ground state of the kinetic part is a Fermi liquid (FL) metal.

III. STRONG BREATHING LIMIT AND TYPE-I CMI
Electron correlations are considered on top of the kinetic part. A strong Hubbard-U merely suppresses double occupation on a single lattice site and cannot cause localization due to the fractional filling here. What replaces is the cluster localization from the inter-site interactions. We first explore the strong breathing limit with V 2 = 0 and study the cluster localization driven by the remaining interaction V 1 in the framework of the slave-rotor construction [28][29][30][31]. A strong repulsion V 1 penalizes the double occupancy on the up-triangles and would drive a Mott transition from FL metal to a CMI. Because the number of the up-triangles is equal to the total electron number, there is exactly one electron in each up-triangle in this cluster Mott regime. To describe different phases and study the Mott transitions, we first employ the standard slave-rotor representation for the electron operator where the bosonic rotor (e iθi ) carries the electron charge and the fermionic spinon (f † iσ ) carries the spin quantum number. To constrain the enlarged Hilbert space, we introduce an angular momentum variable L z i , where L z i is conjugate to the rotor variable with Since the interaction U is assumed to be the largest, in the large U limit the double electron occupation is always suppressed. Hence, the angular variable L z i primarily takes L z i = 1/2 (−1/2) for a singly-occupied (empty) site. Decoupling the electron hopping into the spinon and the rotor sectors, we obtain the Hamiltonians for the spin and charge sectors where r labels the center of up-triangle or equivalently the unit cell of the lattice, and h is a Lagrangian multiplier to enforce the Hilbert space constraint. The Hamiltonian is invariant under Here we have introduced an angular momentum operator L r as where the lattice site is labeled by the combination of the unit cell r and the sublattice index µ. L r measures the total electron occupation on the up-triangle at r. Moreover, from L rµ , we find that L r can take −1, 0, 1, 2. Finally, because L z i = ±1/2, the last term in the second line of Eq. (7) reduces to a constant and can be dropped. It is convenient to define the conjugate variable for L r . We introduce a super-rotor operator e ±iΘr whose physical meaning is to create and annihilate an electron charge in the up-triangle at r. Clearly, we have As [θ rµ − θ rν , L r ] = 0, the hopping terms inside the uptriangles commute with the V 1 interaction and we can set θ rµ ≡ Θ r . The hopping terms inside the down-triangles, that describe the electron tunneling from the neighboring up-triangles, do not commute with the V 1 interaction. Increasing V 1 penalizes the kinetic energy gain through hoppings on the down-triangle bonds and causes the electron cluster localization in the up-triangles. Nevertheless, the electrons remain mobile inside each up-triangle and can gain kinetic energy through hoppings within the up-triangle. Thus, the system is locally "metallic" within each up-triangle and remains so even when the interaction V 1 becomes dominant. Using the local metallic condition (θ rµ ≡ Θ r ) to optimize the intra-up-triangle hopping, we obtain a reduced rotor Hamiltonian that is defined on the triangular lattice formed by the centers of the up-triangles, where rr labels two neighboring up-triangles and h = h + 5V 1 /2, J 2 is defined on down-triangles. The relevant degrees of freedom for the Mott transition is the super-rotor mode e iΘr . When it is condensed and e iΘr = 0, we obtain a FL metal. When it is gapped with e iΘr = 0, a CMI with electrons localized on all up-triangles is obtained, and we refer this CMI as type-I CMI in Fig. 2(a). In this type-I CMI, there exists charge coherence within the up-triangle as it is "locally metallic". The gauge field fluctuations within the up-triangles become massive from the Higgs mechanism. The gauge fluctuations on the links between two neighboring up-triangles remain gapless and we represent it by a rr for two up-triangles at r and r . The reduced rotor HamiltonianH c and the spinon Hamiltonian In the type-I CMI, the spinon mean-field Hamiltonian H s describes the spinon hopping at the mean-field level. The spinon bands are identical to the electronic ones, E µ (k), except for the modified hopping. Thus, the spinons fill a half of the lowest spinon band, leading to a spinon Fermi surface. The resulting spin sector is a U(1) quantum spin liquid (QSL) with a spinon Fermi surface. It is generally believed that, the U(1) QSL is in the deconfined phase due to the spinon Fermi surface that suppresses the instanton events. When the super-rotor mode is condensed, the U(1) gauge field picks up a mass via the Higgs' mechanism, and the charge rotor and fermionic spinons are then combined back to the original electron.
Here we solve the charge sector HamiltonianH c and the spinon Hamiltonian H s self-consistently for the phase diagram and Mott transition. Following the standard procedure, we implement the coherent state path integral for the super-rotor variables φ † r ≡ e iΘr and φ r ≡ e −iΘr . By integrating out the field L r , we obtain the partition function, with the effective action We have dropped the term with parameterh that is required to vanish since r L r = 0. The Lagrange multiplier λ r is also introduced in the partition function to enforce the unimodular constraint |φ r | = 1 for each uptriangle. We take a uniform saddle point approximation by setting λ r = λ and further integrate out the φ fields. Finally we get the following saddle point equation in the Mott insulating phase, where S BZ is the area of the first Brillouin zone of the triangular lattice and ω k is the dispersion of the superrotor mode with When λ = 6J 2 , the dispersion ω k becomes gapless. That means the super-rotor mode is condensed. Combining   2. (a) The slave-rotor mean-field phase diagram at V2 = 0. We exclude the 120-degree state in the strong coupling limit (V1 t2). Inset describes the free and uncorrelated motion of the electrons inside the up-triangles, and the direction is arbitrarily chosen. (b,c) The phase diagram of the extended Hubbard model for different parameters.
this condensation condition with the super-rotor saddle point equation Eq. (15) and the spinon-sector mean-field theory, we construct the phase diagram in the strong breathing limit as shown in Fig. 2(a).
Here we do not consider the possibility of magnetic ordering in the strong Mott regime. For a small (large) V 1 /t 2 , we obtain a Fermi liquid metal (or a U(1) QSL with a spinon Fermi surface). The Mott transition is continuous and of the quantum XY type in the meanfield theory, and is expected to be so even after including the U(1) gauge fluctuations [32]. The phase boundary of the Mott transition is understood as follows. For smaller (larger) t 1 /t 2 , the electrons gain more (less) kinetic energy from the t 2 hopping or the inter-up-triangle hopping, and thus, a larger (smaller) critical V 1 /t 2 is needed to localize the electrons in the up-triangles. In particular, in the limit of t 1 /t 2 → ∞, our model with V 2 = 0 and 1/6 electron filling is equivalent to a triangular lattice Hubbard model at half-filling where the triangular lattice is formed by the up-triangles. Therefore, the U(1) QSL with a Fermi surface in the type-I CMI is smoothly connected to the one proposed for the triangular lattice Hubbard model at half-filling [29,30].

IV. EMERGENT U(1)c GAUGE STRUCTURE IN TYPE-II CMI
As V 2 gradually increases from zero, the free motion of electrons inside the up-triangles becomes less favorable energetically because this motion creates double occupancy configurations on the down-triangles. Thus, at a critical V 2 , the electron number on each down-triangle is also fixed to be one, and we experience the cluster localization on both types of triangles. This new cluster Mott state is referred as type-II CMI. The slave-rotor representation in this phase fails to capture the proper physics and we should introduce a new parton construction. Be-fore that, we need to first understand the low-energy physics of the charge sector, especially in the type-II CMI. We will show the charge localization pattern in the type-II CMI leads to an emergent compact U(1) lattice gauge theory description for the charge-sector quantum fluctuations. In the slave-rotor formalism, the chargesector Hamiltonian with V 2 interaction is given by where we have dropped the U interaction term because L i = ±1/2 only gives a constant for the L 2 i term in the large U limit. Up to a mapping from the rotor operators to the spin ladder operators, e ±iθi = L ± i , this charge sector Hamiltonian is equivalent to a Kagomé lattice spin-1/2 XXZ model in the presence of an external magnetic field. We can recast the charge sector Hamiltonian as Here the effective spin-1/2 ladder operators L ± i satisfy and the effective "magnetic" field reads B eff ≡ h + 3(V 1 + V 2 ). The 1/6 electron filling can be regarded as the total "magnetization" condition N s −1 i L z i = −1/6, where N s is the total number of Kagomé lattice sites.
When the inter-site repulsion V 1,2 dominate over the hoppings t 1,2 , the system enters into the type-II CMI phase where the cluster localization appears on both types of triangles. In terms of the effective spin L z i , the electron charge localization condition in the type-II CMI is Therefore, the allowed effective spin configuration is "2down 1-up" in every triangle. These allowed classical spin configurations are extensively degenerate. But the degeneracy will be further lifted after involving the transverse effective spin exchanges. Physically, the effective interactions are originated from the collective hopping processes of electrons (shown in Fig. 3) and can be obtained from a third order degenerate perturbation theory. The resulting ring exchange Hamiltonian has the form of where " " refers to the elementary hexagon of the Kagomé lattice, J ring = 6J 3 1 /V 2 2 + 6J 3 2 /V 2 1 is the ring exchange parameter and "1, 2, 3, 4, 5, 6" are the six vertices on the corner of the elementary hexagon on the Kagomé lattice (see Fig. 3).
We now demonstrate that the effective Hamiltonian H c,ring can be mapped into a compact U(1) lattice gauge theory on the dual honeycomb lattice. As shown in Fig. 1(a), this dual honeycomb lattice is formed by the centers of up-and down-triangles, labeled as r and r respectively. We follow the previous work Ref. 33 and introduce the lattice U(1) gauge fields (E, A) by defining where r ∈ u, E rr = −E r r , and A rr = −A r r . The fields E and A are identified as the electric field and the vector gauge field of the compact U(1) lattice gauge theory and satisfy [E r,r+eµ , A r,r+eµ ] = −i. With this identification, the local "2-down 1-up" charge localization condition in Eq. (20) is interpreted as the "Gauss' law" for the emergent U(1) lattice gauge theory. The effective ring exchange Hamiltonian H c,ring reduces to a gauge "magnetic" field term on the dual honeycomb lattice, where ∆ × A is a lattice curl defined on the " " that refers to the elementary hexagon on the dual honeycomb lattice. As this internal gauge structure emerges at the low energies in the charge sector, we refer this gauge field as the U(1) c gauge field. The fate of U(1) c gauge field can become confining as it is in two spatial dimensions. However, as the gapless spinon matter is involved, it is likely that the instanton events can still get suppressed. Thus, even though the plaquette charge order is expected in previous works [12,[34][35][36][37][38], the gauge deconfinement can coexist with the charge order. This is very much like the AFM * phase where the spin quantum number fractionalization and the antiferromagnetic order coexist [39,40]. The more detailed structure inside the charge sector of the type-II CMI is not the focus of this work. We are more concerned about the localization pattern and thus assume the charge fractionalization in the type-II CMI. A strict analysis requires non-perturbative computations involving quantum fluctuations and is beyond the meanfield theory. In this work, we ignore the instanton effect and focus on constructing the mean-field phase diagram for the extended Hubbard model.

V. MEAN-FIELD THEORY FOR THE TRANSITION BETWEEN TYPE-I TO TYPE-II CMIS
In this section, we go beyond the strong breathing limit and build a generic framework which can support both type-I and type-II CMIs. As mentioned in the last section, the slave-rotor representation is incapable of describing the type-II CMI state where cluster localization occurs on both types of triangles. Therefore, we first introduce a new parton representation based on the emergent gauge structure in type-II CMI and then establish the phase diagram at the mean-field level. We also discuss the properties in each clusterization phase and the phase transitions between them.

A. Slave-particle Construction and Mean-field Theory
To study the transition between two distinct cluster localization states, we return to the charge sector Hamiltonian in Eq. (7) by adding the V 2 interaction, Since the electron is not localized on a lattice site in the CMIs, the rotor variable e iθi is insufficient to describe all the phases and phase transitions, except for the special limits for type-I CMI that we have analyzed. To fix the problem, we extend the slave-rotor representation to a new parton construction for the electron operator [1,31,41], where e µ connects the up-triangle center r and the neighboring down-triangle centers r + e µ , Φ † r (Φ r+eµ ) creates (annihilates) the bosonic charge excitation in the triangle at r (r + e µ ), and l ± r,r+eµ ≡ |l ± r,r+eµ |e ±iAr,r+e µ is an open string operator of the U(1) c gauge field in the charge sector connecting the charge excitations in the neighboring triangles at r and r + e µ . Under the U(1) c gauge transformation, Φ † r → Φ † r e iχr , Φ r+eµ → Φ r+eµ e −iχr+e µ , and A r,r+eµ → A r,r+eµ e −iχr+iχr+e µ . To constrain the Hilbert space of the parton construction, one defines the following operator [31], that measures the local U(1) c (electric) gauge charge, and for the remaining part of the paper, r refers to the centers of both up (denoted as 'u') and down (denoted as 'd') triangles. Here, η r = +1 (−1) for r ∈ u (r ∈ d) and l r,r+ηreµ = L z r,r+ηreµ . We further supplement this definition with a Hilbert space constraint [31,41], such that the physical Hilbert space is recovered. For the type-II CMI, Q r = 0 for every triangle. Due to the single electron occupancy on all triangles for type-II CMI, the electron motions are correlated in type-II CMI instead of the free electron motion in the inset of Fig. 2(a) for a type-I CMI. This correlated electron motion leads to the emergent U(1) c gauge structure here [12] and the plaquette charge order whose consequences on the spin sectors are explained in previous works [12,[34][35][36][37][38]. This is not the focus of this work where we are more concerned about the distinct types of cluster localization. Using the new parton construction, the Hubbard model becomes which is supplemented with the Hilbert space constraint.
In mean-field treatment, we decouple the kinetic terms and H u c for the charge sector in the up-triangles, H d c for the charge sector in the down-triangles, H s for the spinon, and H l for the U(1) gauge link, Type-II CMI Φr = 0 for r ∈ u, d.
FL metal Φr = 0 for r ∈ u, Φr = 0 for r ∈ d. Lagrange multiplier λ r is used to implement the unimodular constraint for the Φ field at each r site. The effective action S i for the up-and down-triangle subsystems are where rr refers to the nearest-neighbor sites on each type of triangle subsystem. The rest of the treatment on each subsystems is identical to what we did to the superrotor mode in Sec. III and then we can find the critical V i /J i at which the bosons are condensed. The resemblance between the above actions and the the action of Eq. (14) indicates the close connection between this slaveparticle approach used here and the slave-rotor formulation used in Sec. III. In the strong breathing limit with V 2 = 0, this two approaches should give qualitatively the same results. But quantitatively, the current approach, through the string parameters, takes into account of the reduction of the spinon or electron bandwidth due to the on-site Hubbard interaction. As a result, we expect that it could give a more reliable phase diagram especially for the FL metal phase. The generic mean-field phase diagrams for different choices of couplings are depicted as Fig. 2(b,c). The four phases correspond to different behaviors of the charge bosons (see Tab. I). When the charge bosons from both up-and down-triangles are condensed, the FL metal is realized. When they are both gapped and uncondensed, we have the type-II CMI. When the charge bosons from one triangle are condensed and the other is uncondensed, we have the type-I CMI. Here the subindex 'u' or 'd' to the type-I CMI indicates which triangles the electrons are localized in. In mean-field theory, because the charge is a higher energy degree of freedom, the spinon sector was treated as a spectator, rather than the driving force.
We turn to explain the phase boundaries in Fig. 2(b,c). As we increase V 2 /t 1 , the effective electron hopping on the up-triangle bonds gets suppressed which effectively enhances the kinetic energy gain through the downtriangle bonds. Thus, a larger V 1 /t 2 is required to drive a Mott transition. A similar argument applies to the boundary between type-I u and type-II CMIs. A larger V 2 /t 1 is needed to compete with the kinetic energy gain on the up-triangle bonds for a larger V 1 /t 2 in type-I u CMI and to drive a transition to type-II CMI. For t 1 > t 2 , electrons are more likely to be localized in the up-triangles to gain the intra-cluster kinetic energy. Thus, a smaller V 1 /t 2 is needed to drive a Mott transition and a larger V 2 is needed to drive the system from type-I u to type-II CMIs.
The phase transition between the tpye-I and type-II CMIs can also be understood in the charge boson picture. To be concrete, we focus on the transition from the type-II CMI to the tpye-I u CMI and the extension to the tpye-I d CMI is direct. In the type-II CMI, charge bosons from both up-and down-triangles are gapped and uncondensed. With the decreasing of the inter-site repulsion V 2 /t 1 on the down-triangle bonds, the charge bosons become condensed on the down-triangle subsystem and the U(1) s gauge field also acquires a mass concurrently. The two fractionally-charged charge bosons Φ from two types of triangles then are combined back into the original unit-charged charge rotor e iθ . The large inter-site repulsion V 1 /t t still preserves the single electron occupancy on the up-triangle subsystem. Thus, the charge rotor e iθ is well-defined on the center of an up-triangle and within the up-triangle, the localized electron can move more or less freely. In this sense, the condensation of the charge bosons from the down-triangles leads to the local "metallic" clusters in the up-triangles. After the charge boson condensation, there is no charge fractionalization in the type-I CMI, but the spin-charge separation still survives. Because of the local "metallic" clusters, only the U(1) s gauge field living on the down-triangle bonds that connect the up-triangles remains active and continues to fluctuate at the low energies. The low-energy physics is described by the spinon Fermi surface coupled with a fluctuating U(1) s gauge field, leading to a U(1) QSL in the triangular lattice formed by the up-triangles.
All transitions discussed above are continuous at meanfield level, except the transition between type-I u and type-I d CMIs that is strongly first order. Beyond mean-field theory, the transition between FL metal and type-I CMIs will remain continuous and quantum XY type [32,42] while the transition into type-II CMI may depend on the detailed charge structure inside type-II CMI. Moreover, our mean field theory does not capture the charge quantum fluctuation inside the type-II CMI as described by the compact U(1) c gauge theory in Sec. IV, but does obtain qualitatively correct phase boundaries.

VI. DISCUSSION
We discuss the experimental relevance and consequences about the Mo-based cluster magnets. These compounds, M 2 Mo 3 O 8 (M = Mg, Mn, Fe, Co, Ni, Zn, Cd), LiRMo 3 O 8 (R = rare earth) and other related variants [43][44][45][46], incorporate the Mo 3 O 13 cluster unit, and the physical properties of most materials have not been carefully studied so far. According to our theory, more anisotropic systems with a stronger breathing tend to favor the type-I CMI. Li 2 InMo 3 O 8 is more anisotropic than LiZn 2 Mo 3 O 8 from the lattice parameters. For LiZn 2 Mo 3 O 8 , the spin susceptibility shows a "1/3 anomaly" and double Curie regimes [22,26]. This was attributed to the plaquette charge order in the type-II CMI that reconstructs the spin sector. In contrast, Li 2 InMo 3 O 8 is characterized by one Curie regime with the Curie temperature Θ CW = −207K down to 25K [43]. The Curie constant is consistent with one unpaired spin-1/2 moment per Mo 3 O 13 cluster in the type-I CMI. Below 25K, the spin susceptibility of Li 2 InMo 3 O 8 saturates to a constant, which is consistent with the expectation from a spinon Fermi surface U(1) QSL. Besides the structural and spin susceptibility data, however, very little is known about Li 2 InMo 3 O 8 . It is also likely that this system is located in the 120-degree order state of the model. Thus, more experiments are needed to confirm the absence of magnetic ordering in Li 2 InMo 3 O 8 and also to explore the magnetic properties of ScZnMo 3 O 8 and other cluster magnets [47]. From the numerical aspect, first principal calculation is carried out recently and supports the experimental findings and theoretical understandings in Mo 3 O 8 magnets. Namely, LiZn 2 Mo 3 O 8 is shown to exhibit the plaquette order with one dangling spin, Li 2 InMo 3 O 8 is a CMI with 120-degree order and Li 2 ScMo 3 O 8 displays a spin liquid behavior [48].