Tripartite entangled plaquette state in a cluster magnet

TL;DR

This paper uses quantum Monte Carlo simulations to identify a tripartite entangled plaquette state in a 2D anisotropic Kagome lattice spin model, revealing its phase transitions and potential relevance to quantum materials and optical lattice systems.

Contribution

It demonstrates the existence of a tripartite entangled plaquette state in a Kagome lattice XXZ model and explores its phase diagram and transitions, connecting to experimental systems.

Findings

Identified a tripartite entangled plaquette state preserving all Hamiltonian symmetries.

Mapped out the phase diagram and phase transitions involving this state.

Connected the findings to charge dynamics in a Hubbard model and optical lattice realizations.

Abstract

Using large-scale quantum Monte Carlo simulations we show that a spin-$1/2$ XXZ model on a two-dimensional anisotropic Kagome lattice exhibits a tripartite entangled plaquette state that preserves all of the Hamiltonian symmetries. It is connected via phase boundaries to a ferromagnet and a valence-bond solid that break U(1) and lattice translation symmetries, respectively. We study the phase diagram of the model in detail, in particular the transitions to the tripartite entangled plaquette state, which are consistent with conventional order-disorder transitions. Our results can be interpreted as a description of the charge sector dynamics of a Hubbard model applied to the description of the spin liquid candidate ${\mathrm{LiZn}}_{2}{\mathrm{Mo}}_{3}{\mathrm{O}}_{8}$, as well as a model of strongly correlated bosonic atoms loaded onto highly tunable {\it trimerized} optical Kagome lattices.

Figures

Equations

$${\mathcal{H}}_{\text{c}}$$

$$S_{\text{q}}^{\alpha }=\frac{1}{N_{\text{s}}}\sum _{{\text{r}}_i}{e}^{i\text{q}\cdot \left({\text{r}}_i+\text{upalpha}\right)}{S}_{{\text{r}}_i+\text{upalpha}}^z,$$

$$n_{\text{q}}^{\alpha \beta }=\frac{1}{N_{\text{s}}}\sum _{{\text{r}}_i{\text{r}}_j}{e}^{i\text{q}\cdot \left[\left({\text{r}}_i+\text{upalpha}\right)-\left({\text{r}}_j+\text{upbeta}\right)\right]}⟨S_{{\text{r}}_i+\text{upalpha}}^{+}{S}_{{\text{r}}_j+\text{upbeta}}^{-}⟩.$$

$$B{B}_{\text{q}}^{\alpha \beta }=\frac{1}{N_{\text{s}}}\sum _{{\text{r}}_a{\text{r}}_b}{e}^{i\text{q}\left({\text{r}}_a-{\text{r}}_b\right)}⟨B_{{\text{r}}_a}^{\alpha }{B}_{{\text{r}}_b}^{\beta }⟩,$$

$${\chi }_t=3\beta \left[⟨{\left(\frac{1}{3}\sum _{i\in t}{n}_i\right)}^2⟩-{⟨\left(\frac{1}{3}\sum _{i\in t}{n}_i\right)⟩}^2\right],$$

$$⟨S_i^z{S}_j^{+}{S}_k^{-}{S}_{i^{\prime }}^z{S}_{j^{\prime }}^{+}{S}_{k^{\prime }}^{-}⟩$$

$$=\frac{4}{{\left(\beta J_{\pm }\right)}^2}⟨\left(n-1\right)\sum _t{S}_i^z[t]S_{i^{\prime }}^z[t_{\text{next}}]⟩,$$

Full text

I INTRODUCTION

Frustration in the context of magnetism refers to the phenomena where interactions between magnetic moments compete at the microscopic level, usually due to a combination of antiferromagnetic exchange and lattice geometry. Unlike unfrustrated systems, where symmetry breaking in the groundstate prevails, the inability of frustrated magnets to satisfy each and every microscopic interaction can lead to the emergence of exotic groundstate phases, such as valence bond solids,Tamura et al. (2006); Matan et al. (2010); Sachdev (1999); Vasiliev et al. (2015) spin liquids,Pratt et al. (2011); Han et al. (2012); Balents (2010); Yan et al. (2011) classicalBramwell and Gingras (2001) and quantum spin ices,Gingras and McClarty (2014); Ross et al. (2011) among many others.

Geometric frustration typically arises in magnetic moments that are localized on single ions, through their spatial arrangement and that of their exchange interactions. However, moments on ionic crystal lattices are often susceptible to structural distortions, orbital couplings, or mixing between magnetic and non-magnetic layers. These and other perturbations may disrupt the formation of delicate exotic phases, leading in many cases to conventional ordering. This difficulty may be alleviated in the recently-discovered geometrically frustrated antiferromagnets where the magnetic moments are localized on small transition-metal clusters, rather than being localized on a single ion.Sheckelton et al. (2012) These materials have been shown to avoid the key limitations of ion-localized moments mentioned above, making them natural candidates to search for exotic states of matter.Sheckelton et al. (2012); Mourigal et al. (2014); Sheckelton et al. (2014); Haraguchi et al. (2015); Chen et al. (2014, 2016) In a recent development, Sheckelton et alSheckelton et al. (2012); Mourigal et al. (2014); Sheckelton et al. (2014) found that the molecular magnet on the triangular lattice ${\mathrm{LiZn}}_{2}{\mathrm{Mo}}_{3}{\mathrm{O}}_{8}$ exhibits spin liquid behavior with low-energy spin correlations consistent with the highly-coveted resonating valence-bond solid state.Anderson (1987) Similarly, evidence of strong quantum fluctuations and spin liquid behavior was found in a related cluster magnet ${\mathrm{Li}}_{2}{\mathrm{ScMo}}_{3}{\mathrm{O}}_{8}$, while the isomorphic compound ${\mathrm{Li}}_{2}{\mathrm{InMo}}_{3}{\mathrm{O}}_{8}$ was found to develop long-range 120*∘* magnetic order.Haraguchi et al. (2015)

Attempts to elucidate the microscopic origin of the experimental observations in ${\mathrm{LiZn}}_{2}{\mathrm{Mo}}_{3}{\mathrm{O}}_{8}$ have included a model of lattice distortions leading to an emergent honeycomb lattice where the spins form a quantum spin liquid,Flint and Lee (2013) as well as a purely electronic description based on a $1/6$-filled extended Hubbard model with nearest-neighbor repulsion on a trimerized Kagome lattice.Chen et al. (2016) The later work suggested that the ground state of ${\mathrm{LiZn}}_{2}{\mathrm{Mo}}_{3}{\mathrm{O}}_{8}$ may be a U(1) spin liquid with plaquette charge order and a spinon Fermi surface, whose finite-temperature properties may explain the two surprising Curie-Weiss regimes observed in the experimental data.Chen et al. (2016)

The U(1) spin liquid state arises from a generic procedure where a mean-field decoupling of the charge and spin degrees of freedom in terms of a slave-rotor representation of the electron operatorsFlorens and Georges (2004) is performed. In such an approach, the electron systems are mapped onto a spinon Hamiltonian coupled to a bosonic lattice model of the charge sector via mean-field parameters. The intuition behind such an approach is in the observation that in certain strongly-coupled electron systems, the dynamics of the spin and charge degrees of freedom is markedly different. Therefore, the electron may be better understood as being composed of separate charge and spin variables. To determine the fate the ground state of the overall fermionic system, both spinon and bosonic Hamiltonians have to be solved simultaneously. Since the resulting bosonic Hamiltonian associated with the charge sector is generically strongly interacting, such a problem can be solved, for instance, via a standard mean-field decoupling. However, other approaches that include some spatial correlations, e.g., quantum Monte Carlo (QMC), exact diagonalization, or density-matrix renormalization group methods, are desirable and clearly improve the quality of the description of the many-body electron problem under study.Zhao and Paramekanti (2007)

In this work, we use large-scale QMC simulations to show that a two-dimensional spin-$1/2$ model on an anisotropic Kagome lattice–proposed as the description of the charge sector of an extended Hubbard model applied to the cluster magnet ${\mathrm{LiZn}}_{2}{\mathrm{Mo}}_{3}{\mathrm{O}}_{8}$–exhibits three different phases: a ferromagnet (FM), a fractionally-filled tripartite entangled plaquette state (W) and a valence-bond solid (VBS). These phases are arranged in the phase diagram presented in Fig.1(a). We find that the FM-to-W insulator transition is continuous and belongs to the three-dimensional $O(2)$ universality class if the magnetization is kept constant across the transition;Fisher et al. (1989) instead, if the magnetization varies, the transition is generically mean-field.Fisher et al. (1989) We also find evidence indicating that the transition between the W state and the VBS, up to the system sizes accessible with QMC, appears to be first order, in agreement with the theoretical expectation that it belongs to the three-dimensional three-state clock universality.Herrmann (1979) Finally, we re-examine the transition between the FM and the VBS in terms of a recently proposed scenario of quantum criticality with two length scales Shao et al. (2016). We find that the presence of irregular system-size dependent oscillations in the observables prevents us from drawing a firm conclusion about the applicability of this scenario to the FM-to-VBS transition.

II MODEL AND QUANTUM MONTE CARLO SIMULATIONS

We begin by motivating a Hamiltonian for the charge sector of the cluster magnet ${\mathrm{LiZn}}_{2}{\mathrm{Mo}}_{3}{\mathrm{O}}_{8}$. This material is formed by small triangular ${\mathrm{Mo}}_{3}{\mathrm{O}}_{13}$ units, where each triangular plaquette accommodates one unpaired electron. The small ${\mathrm{Mo}}_{3}{\mathrm{O}}_{13}$ triangular units are located on the sites of a triangular lattice shown in Fig. 1(b). Thus the ${\mathrm{Mo}}$ atoms can be thought of as forming a trimerized Kagome lattice. Following Ref. Chen et al., 2016, we consider the electrons hopping on the Kagome lattice as described by an extended Hubbard model at 1/6 electron filling. In this model, both the on-site and the nearest-neighbor Coulomb interactions are included in addition to the electron hopping. The authors of Ref. Chen et al., 2016 employ a standard slave-rotor representation of the constituent fermions.Florens and Georges (2004) The electron operator is reformulated as the product of a U(1) charge rotor variable and a fermionic spinon, $c^{{\dagger}}_{{\bf r}\sigma}=e^{i\theta_{\bf r}}f^{{\dagger}}_{{\bf r}\sigma}$, where the bosonic rotor $e^{i\theta_{\bf r}}$ (the fermionic spinon $f^{{\dagger}}_{{\bf r}\sigma}$) creates an electron charge (a spinon with spin $\sigma$) at lattice site ${\bf r}$. Since the Hilbert space has been enlarged, one has to introduce a constraint to get back to the physical Hilbert space through an angular momentum variable: $S^{z}_{\bf r}=\sum_{\sigma}f^{{\dagger}}_{{\bf r}\sigma}f_{{\bf r}\sigma}-1/2$. The local electron Hilbert space is thus represented as $\lvert 0\rangle_{c}={\lvert 0\rangle}_{f}{\lvert{S^{z}=-1/2}\rangle}_{\theta}$, ${\lvert\uparrow\rangle}_{c}={\lvert\uparrow\rangle}_{f}{\lvert S^{z}=1/2\rangle}_{\theta}$, ${\lvert\downarrow\rangle}_{c}={\lvert\downarrow\rangle}_{f}{\lvert S^{z}=1/2\rangle}_{\theta}$, and ${\lvert\uparrow\downarrow\rangle}_{c}={\lvert\uparrow\downarrow\rangle}_{f}{\lvert S^{z}=3/2\rangle}_{\theta}$. When the on-site Coulomb interaction is dominant, the double electron occupancy on a single site is forbidden and the operator $S^{z}$ describes an effective spin-1/2 angular momentum operator that is conjugate to the charge rotor variable, i.e., $[\theta_{\bf r},S^{z}_{{\bf r}^{\prime}}]=i\delta_{{\bf r},{\bf r}^{\prime}}$. Thus the rotor operators can be identified as the spin ladder operators $S^{\pm}_{\bf r}=e^{\pm i\theta_{\bf r}}$. Using a mean-field decoupling, the original Hubbard model is transformed into two Hamiltonians for the spinon and charge sectors coupled via mean-field parameters. Under the slave-rotor reformulation, the charge sector of the extended Hubbard model is now described by an effective spin-1/2 model on an anisotropic Kagome lattice with

[TABLE]

where the Ising exchange interaction (transverse exchange interaction) accounts for the nearest-neighbor Coulomb interaction (electron tunneling). The index $t=1,2$ indicates the different up (1, shaded with solid lines) and down (2, white with dashed lines) triangles forming the Kagome lattice, as shown in Fig.1(b). The anisotropy of the lattice is encoded in the coupling constants $J_{\pm t}$, $J_{zt}$, whereas the effective magnetic field $h_{\text{eff}}$ controls the average effective magnetization. We will consider systems at both strictly conserved magnetization as well as at fixed magnetic fields. We set $J_{\pm 1}=J_{\pm 2}=J_{\pm}$ as the reference energy scale of the problem. This model can also be thought of as one describing hardcore bosonsMatsuda and Matsubara (1957) loaded on a trimerized Kagome lattice. In the hardcore-boson language, the FM phase corresponds to a superfluid (SF), while the W phases correspond to fractionally-filled Mott insulators where a hardcore bosons are localized on the triangles with the largest $J_{z}$. Such a system of hardcore bosons on trimerized lattice geometry can be realized using superlattice techniques in ultracold gases.Santos et al. (2004)

The Hamiltonian in Eq.(1) cannot be solved analytically, but large-scale QMC simulations are allowed in the sign-problem free regime where $J_{\pm t}>0$. This is the natural choice if we interpret Eq.(1) as a system of hardcore bosons. We develop a finite-temperature Stochastic Series Expansion Sandvik (1999); Syljuåsen and Sandvik (2002); Melko (2007) (SSE) QMC algorithm with directed loop updates. We map out the phase diagram of the model through measurements such as magnetization $m_{z}=\langle\hat{m}\rangle=\langle\frac{1}{V}\sum_{\text{\bf r}}S_{\text{\bf r}}^{z}\rangle$, uniform spin susceptibility $\chi_{z}=\frac{V}{T}\left(\langle\hat{m}^{2}\rangle-\langle\hat{m}\rangle^{2}\right)$, superfluid stiffness $\rho_{s}$,Pollock and Ceperley (1987) and a set of diagonalSandvik (1999) and off-diagonalDorneich and Troyer (2001) spin structure factors. Figure 1(a) shows the QMC phase diagram for the model of Eq. (1) extracted from the finite-temperature and the finite-size scaling of the superfluid stiffness and the diagonal structure factor, performed up to lattice sizes of $V=L\times L\times 3=60\times 60\times 3$ and inverse temperature of up to $\beta=J_{\pm}/T=60$. The phase diagram presented in Fig.1(a) is obtained at strictly enforced fixed magnetization $m_{z}=-1/6$ that corresponds to the 1/6 electron filling in LiZn2Mo3O8. To do this efficiently, we first tune the magnetic field $h_{\text{eff}}$ such that the average magnetization is as close as possible to the desired magnetization sector. We subsequently run simulations that are still grand-canonical but whose measurements are taken only at configurations that are in the desired magnetization sector.

The grand-canonical phase diagram of the isotropic $J_{z1}=J_{z2}$ case has been explored in Ref. Isakov et al., 2006 where the authors found two phases: a FM and a VBS phase. It was found that the FM phase is characterized by long-range in-plane magnetic order with wave vector $\text{\bf q}={\bf 0}$, finite superfluid stiffness $\rho_{s}>0$, and uniform susceptibility $\chi_{z}>0$. The VBS was studied in great detail in Ref. Isakov et al., 2006: it is a gapped, translationally broken three-fold degenerateNikolić and Senthil (2005) phase, where spins resonantly flip in each hexagon marked by red circles in Fig.(1)(b). The remaining spins anti-align along the $z$ direction and their wave vector is $\text{\bf q}=\text{\bf K}=(2\pi/3,0)$.Isakov et al. (2006) The average magnetization is $m_{z}=-1/6$ and $\chi_{z}=0$.

III The tripartite entangled plaquette state

The gapped, tripartite entangled plaquette states W preserve all the symmetries of the Hamiltonian and are not present in the isotropic $J_{z1}=J_{z2}$ caseIsakov et al. (2006). We find that, in the strong triangle limit, the W states reduce to simple product states where the three spins on the strong triangle (up $\text{W}_{\bigtriangleup}$ if $J_{z1}\gg J_{z2}$ or down $\text{W}_{\bigtriangledown}$ if $J_{z2}\gg J_{z1}$) form tripartite entangled three-qubit W states.Santos et al. (2004); Damski et al. (2005); Dür et al. (2000) The three-qubit states are given by $|\text{W}\rangle=\left(|{\uparrow\downarrow\downarrow}\rangle+|{\downarrow\uparrow\downarrow}\rangle+|{\downarrow\downarrow\uparrow}\rangle\right)/\sqrt{3}$.

To understand the properties of the W phases and the phase transitions to the nearby FM and VBS phases, we measure several correlation functions. We measure the diagonal spin structure factor $S^{\alpha\beta}_{\text{\bf q}}/N_{\text{s}}=\langle S^{\alpha}_{\text{\bf q}}S^{\beta}_{-\text{\bf q}}\rangle-\langle S^{\alpha}_{\text{\bf q}}\rangle\langle S^{\beta}_{-\text{\bf q}}\rangle$, where

[TABLE]

and $N_{s}$ is the total number of unit cells. Here, $\text{\bf r}_{i}$ points to the sites of the underlying triangular Bravais lattice and the vector $\bm{\upalpha}$ refers to the position of each site within the unit cell with respect to the vector $\text{\bf r}_{i}$. Its purpose is to detect the long-range diagonal order: if the system magnetically orders then $S_{\text{\bf q}}=\sum_{\alpha}S^{\alpha\alpha}_{\text{\bf q}}$ will scale with system size for at least one value of q. We also consider the off-diagonal spin structure factor Dorneich and Troyer (2001)

[TABLE]

Using this quantity, we also study the equivalent of the condensate fraction in Bose systems,Penrose and Onsager (1956); Giamarchi et al. (2008) defined as the ratio of the “zero-momentum occupation” to the volume of the system $f_{0}=n_{\text{\bf 0}}/V=\sum\limits_{\alpha}n^{\alpha,\alpha}_{\text{\bf 0}}/V$. Finally, we consider the bond-bond structure factor using a four-point correlation function

[TABLE]

where $B^{\alpha}_{\text{\bf r}_{a}}=S_{i_{a\alpha}}^{+}S_{j_{a\alpha}}^{-}+S_{i_{a\alpha}}^{-}S_{j_{a\alpha}}^{+}$. Nearest neighbor sites $i_{a\alpha}$ and $j_{a\alpha}$ belong to one of 6 bonds $\alpha$ in a unit cell located at position $\text{\bf r}_{a}$. If bond order develops then $BB_{\text{\bf q}}=\sum_{\alpha}BB^{\alpha\alpha}_{\text{\bf q}}$ should scale with system size for at least one value of q, with which we define $B_{\text{\bf q}}=BB_{\text{\bf q}}/V$.

The spin structure factors of the $\text{W}_{\bigtriangledown}$ phase are presented in Fig.(2). We find none of the structure factors (off-diagonal $n_{\text{\bf q}}$ Fig.(2)a, diagonal $S_{\text{\bf q}}$ Fig.(2)b, and bond $BB_{\text{\bf q}}$ Fig.(2)c) displays peaks that signal the translational symmetry breaking. This is confirmed in the corresponding size scaling of selected peaks in Fig.(2)d through Fig.(2)e. For instance, $S_{\text{\bf K}}$ and $BB_{\text{\bf K}}$, which are non-vanishing in the VBS phase, quickly go to zero in the $\text{W}_{\bigtriangledown}$ phase, within the precision of our calculations.

We have removed the zero-momentum peak of the bond-bond structure factor since $B_{\text{\bf 0}}$ remains finite as $V\to\infty$ in all phases. However, a closer look at the different sublattice $BB^{\alpha\beta}_{\text{\bf 0}}$ could reveal a potential rotational symmetry breaking in the unit cell of the Kagome lattice. In Fig.(3) we present results for the finite-size scaling of $BB^{\alpha\beta}_{\text{\bf 0}}$ for the different bonds in the unit cell and their extrapolations to the thermodynamic limit. In the inset the extrapolations to the thermodynamic limit reveal that there is no symmetry breaking in the unit cell within the error bars of our simulation. In the hardcore boson language, the real-space pattern extracted from the extrapolations shows that the bosons are predominantly delocalized along the three sites of the strong triangles (1 boson per down triangle, in this case) since this pattern can be understood as the real space distribution of the square of the average kinetic energy along the bonds in the unit cell. A direct consequence of this observation is that charge fluctuations in the weak triangle are larger than in the strong one. We quantify these fluctuations through a “local” uniform susceptibility

[TABLE]

where $n_{i}=S^{z}_{i}+1/2$. In particular, in the $\text{W}_{\bigtriangledown}$ phase for $J_{z2}/J_{\pm}=4.5$, $J_{z1}/J_{\pm}=2$, $\beta=J_{\pm}/60$, and $V=60\times 60\times 3$ we find that $\chi_{\bigtriangledown}=0.855\pm 0.001$, whereas a significantly more fluctuating weak triangle $\chi_{\bigtriangleup}=7.386\pm 0.002$ is found. Lastly, we examine the chiral-chiral correlation function $\langle E_{t}E_{t^{\prime}}\rangle$, where $E_{t}={\mathbf{S}_{i_{1}}}\left({\mathbf{S}_{i_{2}}}\times{\mathbf{S}_{i_{3}}}\right)$ ($i_{1},i_{2},i_{3}\in t$). Although the time reversal symmetry is explicitly broken by the finite magnetization in the model, chiral order may still potentially develop. In the limit where $J_{z2}\gg J_{z1}$, spin chirality may develop by correlating the fluctuations of the spins on the strong triangles as a way to offset the diagonal energy terms on the weak triangle. We find, however, that the chiral correlations are only enhanced as the $\text{W}_{\bigtriangledown}$ is approached, but the correlations still decay exponentially fast, as in the FM phase. In the appendix we demonstrate how to measure the chiral-chiral correlations in SSE.

IV PHASE TRANSITIONS

IV.1 FM-to-W phase transition

The phase diagram presented in Fig.1(a) exhibits three types of phase boundaries: FM-to-W, FM-to-VBS, and W-to-VBS phase transitions. We first examine the transition between the W and the FM phase. Since the tripartite entangled plaquette state W is a fully symmetric, then such a quantum critical line is expected to be a conventional $(d+1)$-dimensional order-disorder transition (with $d=2$). We consider two situations. First, the W phase is approached at strictly conserved magnetization. Second, the W phase is approached at a fixed magnetic field $h_{\text{eff}}$. The transition at fixed magnetization is anticipated to be of the $(d+1)$-dimensional $O(2)$ vector model type.Fisher et al. (1989) If the W phase is approached at a fixed $h_{\text{eff}}$, the transition corresponds to the appearance of a dilute fluid of excess quasiparticles or holes on top of the W states where the magnetization deviates from $m=-1/6$, and is generically Gaussian.Fisher et al. (1989) For these types of continuous critical points in two dimensions, the superfluid stiffness scales as $\rho_{s}L^{z}=F_{\rho_{s}}\left(L^{1/\nu}\left(J-J_{c}\right),\beta/L^{z}\right)$,Fisher et al. (1973, 1989); Sandvik (2010) where $F_{\rho_{S}}$ is a scaling function, $z$ the dynamical critical exponent, $\nu$ the correlation length exponent, and $J-J_{c}$ the distance to the critical point in terms of the control parameter $J$. In Fig.4 we scrutinize the finite-size scaling of the superfluid stiffness at the transition between the FM and the $W_{\vartriangle}$ as a function of $J_{z1}/J_{\pm}$ at fixed $J_{z2}/J_{\pm}=2$, $m=-1/6$, and $\beta/L=1$. To produce the plots in Fig.4 we have used the critical exponents of $(d+1)$-dimensional $O(2)$ vector model, i.e., $z=1$, $\nu=0.6717\pm 0.0001$, and $\eta=0.0381\pm 0.0002$.Campostrini et al. (2006) In Fig.4(a) we show the rescaled superfluid stiffness $\rho_{s}L$ as a function of $J_{z1}/J_{\pm}$, which becomes system-size independent at the critical point $J_{z1c}/J_{\pm}=3.3325\pm 0.0001$, as implied by the scaling relation. We numerically extract $F_{\rho_{s}}$ by plotting the $\rho_{s}L$ as a function of $\left(J_{z1}/J_{\pm}-J_{z1c}/J_{\pm}\right)L^{1/\nu}$, where a clear collapse is seen in Fig.4(b). Finally, in the inset in Fig.4(b), we analyze the size scaling of the condensate fraction $f_{0}$, which vanishes as $f_{0}\sim L^{-\left(\eta+1\right)}$ at the critical point.Carrasquilla and Rigol (2012) We find that our data for $f_{0}$ are well described by a straight line when plot as a function of $L^{-\left(\eta+1\right)}$ with the $\eta$ obtained in Ref.Campostrini et al., 2006. Our results are consistent with a transition described by the $(d+1)$-dimensional $O(2)$ vector model.

Similarly, we investigate the generic transition at fixed $h_{\text{eff}}/J_{\pm}=-2.935$, $J_{z2}/J_{\pm}=2$, and $\beta/L^{z}=0.1$. The critical exponents that characterize this transition are $z=2$, $\eta=0.0$, $\nu=0.5,$Fisher et al. (1989) which we use in the following to locate the generic critical point. In Fig.5(a) we show the rescaled superfluid stiffness $\rho_{s}L^{2}$ as a function of $J_{z1}/J_{\pm}$, which becomes system-size independent at the critical point $J_{z1c}/J_{\pm}=3.9428\pm 0.0004$. In Fig.5(b) we extract the scaling function $F_{\rho_{s}}$ across the Gaussian critical point, where we observe again a clear collapse to a unique curve.

In summary, the phase transitions from the FM toward the W states are consistent with the picture of $(d+1)$-dimensional order-disorder transitions and with the picture we have described in Section III, namely, that the W states are insulators adiabatically connected to product states of tripartite entangled plaquette states.

IV.2 FM-to-VBS transition

We now turn our attention to the FM-to-VBS quantum phase transition. This transition has been numerically studied in Ref.Isakov et al., 2006 and Ref.Damle and Senthil, 2006 using grand-canonical SSE Monte Carlo. The authors in Ref.Damle and Senthil, 2006 reported evidence of a continuous phase transition based on data for the Binder cumulant and the superfluid stiffness consistent with a non-Ginzburg-Landau deconfined critical point, but did not rule out the possibility of a weak first-order phase transition. On the other hand, Ref.Isakov et al., 2006 reported finite-size scaling of the superfluid stiffness, structure factor, and kinetic energy histograms on larger system sizes. The authors found evidence of a first-order quantum phase transition. Their strongest evidence in favor of a first-order transition was based on extremely low-temperature histograms of the kinetic energy which exhibited a double peak structure signaling coexistence at the critical point. Here we revisit this critical point using canonical measurements of the superfluid stiffness and structure factor supplemented with several finite-size scaling analyses. First, we examine the superfluid stiffness. Assuming a continuous phase transition where $z=1$, we compute the rescaled superfluid stiffness $\rho_{s}L$ as a function of $J_{z2}/J_{\pm}=J_{z1}/J_{\pm}=J_{z}/J_{\pm}$ with $\beta/L=1$ and explore two finite-size scaling scenarios. First, we consider a conventional scaling scenario described by a divergent length scale $\xi\propto\delta^{-\nu}$, where $\delta=J_{z}-J_{zc}$ controls the distance to the quantum critical point and $\nu$ is the correlation length exponent. Assuming a dynamical exponent $z=1$, for a system of linear size $L$ close to $\delta=0$ the superfluid stiffness is singular and scales as $\rho_{s}(\delta,L)=L^{-z}F_{\rho_{s}}\left(\delta L^{1/\nu},L^{-\omega}\right)$, where $\omega$ is a correction-to-scaling exponent. This relation means that at the critical point $\rho_{s}L=a+bL^{-\omega}$ (with $a$ and $b$ constants), and that $\rho_{s}L-bL^{-\omega}$ becomes system-size independent. In Fig.6(a) we plot the relation $\rho_{s}L-bL^{-\omega}$ vs $J_{z}/J_{\pm}$ for our numerical estimates of $\rho_{s}$, which become approximately system-size independent around $J_{z}/J_{\pm}=3.845\pm 0.004$. Note that the crossing point tends to move slowly towards lower values of $J_{z}/J_{\pm}$ and that the value of $\rho_{s}L\equiv\rho_{s}^{c}L$ between two subsequent system sizes is slowly diverging with the system size, which means that corrections to scaling are significant. A more detailed picture arises by considering the size scaling around the critical point: in Fig.6(b) we plot $\rho_{s}L$ vs linear system size $L$ for different values of $J_{z}/J_{\pm}$ near the critical point. On increasing $J_{z}/J_{\pm}$, we notice the development of strong system-size dependent oscillations in the superfluid stiffness which become stronger as the VBS phase is approached. The period of the oscillations can be traced back to the translational symmetry breaking of the VBS phase since the minima of the $\rho_{s}$ oscillations appear at system sizes that exactly accommodate the wave vector of the VBS pattern, i.e., $\text{\bf q}=\text{\bf K}=(2\pi/3,0)$ (see thin vertical lines in Fig.6(b)). Because of the oscillations, fitting the data (to either all the data or to just the local minima or maxima ) to the scaling form $\rho_{s}L=a+bL^{-\omega}$ produce estimates for $a,b$ and $\omega$ with error bars in the first significant digits. Apart from the conventional possibilities discussed above, the slow divergence of the stiffness near the critical point opens up the possibility for yet another scenario, i.e., one with a continuous transition with two diverging length scales, as recently proposed in Ref. Shao et al., 2016. A prediction from the two-length scale scenario is that if there is a second large-$L$ scale controlled by $\delta L^{1/\nu^{\prime}}$, the superfluid stiffness behaves as $\rho_{s}L=L^{1-\nu/\nu^{\prime}}\left(a+bL^{-\omega}\right)$ at the critical point,Shao et al. (2016) which we now test. In Fig.6(c) we display $\rho_{s}L^{z\nu/\nu^{\prime}}-b/L^{\omega}$, which, once again, becomes approximately system-size independent near the estimate of the critical point $J_{z}/J_{\pm}=3.845\pm 0.004$. The parameters used in Fig.6(c) are obtained from fitting our data for $\rho_{s}\left(L,J_{z}/J_{\pm}\right)$ to the form predicted by the two-length scale scenario $\rho_{s}L=L^{1-\nu/\nu^{\prime}}\left(a+bL^{-\omega}\right)$ presented in Fig.6(d). We find that $\nu/\nu^{\prime}\approx 0.4\pm 0.2$, while $\omega\approx 2\pm 2$. Even though the quality of the fits to the two-length scaling forms is notably better than the conventional continuous scaling, the significance of some of the fitting parameters we obtain is again compromised by the large error bars resulting from the oscillations in the data.

We also discuss the finite-size scaling of the structure factor $S_{\text{\bf q}}$ at momentum $\text{\bf q}=\text{\bf K}=(2\pi/3,0)$, which scales to a finite value in the thermodynamic limit inside the VBS phase and to zero in the FM phase. Assuming a continuous phase transition, the structure factor scales as $S_{\text{\bf K}}L^{z+\eta-2}=F_{S}\left(\left(J-J_{c}\right)L^{1/\nu},\beta/L^{z}\right)$.Isakov et al. (2006) In Fig.7 we analyze numerical data for $S_{\text{\bf K}}$ using size scaling. Assuming $z=1$, Fig.7(a) displays the rescaled $S_{\text{\bf K}}L^{z+\eta-2}$ vs $J_{z}/J_{\pm}$, while in Fig.7(b) we attempt at numerically obtaining the scaling function $F_{S}$ by plotting $S_{\text{\bf K}}L^{z+\eta-2}$ vs $(J_{z}/J_{\pm}-J_{zc}/J_{\pm})L^{1/\nu}$. To produce the collapse in Fig.7(b) we find that $\eta=0.10\pm 0.01$ and $\nu=0.40\pm 0.02$. Even though the collapse looks compellingly consistent with criticality, the absence of a crossing in Fig.7(a) and the appearance of a jump in $S_{\text{\bf K}}L^{z+\eta-2}$ suggests that the transition may indeed be first order. We finalize our analysis by mentioning that we have also performed simulations (not shown) for the same transition across the $J_{z2}/J_{\pm}=J_{z1}/J_{\pm}-0.18$ line reaching the same conclusions.

To conclude this section, we have reexamined the FM-to-VBS transition using a canonical algorithm with which we obtain data for the superfluid stiffness and structure factor near the critical point. Apart from the conventional scenarios, and since this transition may potentially exhibit non-classical behavior with the presence of fractionalized excitations, we have also explored the possibility of a scaling with two length scales, which has been recently introduced in Ref. Shao et al., 2016 in the context of a $J-Q$ model exhibiting such phenomenon. Although our data seem consistent with that scenario, the large error bars associated with fitting numerical data to such predictions highlight the complications inherent to the application of this type of analysis to models like ours (Eq. 1). Furthermore, the scaling of the structure factor, which shows drifting crossings between system sizes, is perhaps consistent with the idea that this transition is first order as suggested in Ref.Isakov et al., 2006. Our numerical estimate for the critical point based on the canonical measurements $J_{zc}/J_{\pm}=3.845\pm 0.004$ is slightly below the previous grand-canonical resultsIsakov et al. (2006) $J_{zc}/J_{\pm}=3.898\pm 0.001$.

V VBS-to-W transition

The phase diagram presented in Fig.1(a) features a VBS-to-W-state quantum phase transition which we now inspect. The VBS present in our model is characterized by a three-fold degenerate ground state where translational symmetry is broken and signaled by the formation of Bragg peaks in the diagonal and bond structure factors at $\text{\bf q}=\text{\bf K}=(2\pi/3,0)$. We thus expect that the related $(d+1)$ classical model showing the same universality class is the three-dimensional three-state, $Z_{3}$, clock model, which is equivalent to the three-dimensional three-state Potts model.Wu (1982); Bonnes and Wessel (2011) This model is known to exhibit a first-order phase transition.Wu (1982); Herrmann (1979) Thus we anticipate that the VBS-to-W transition is first order. In Fig.8 we show finite-size scaling of the diagonal structure factor $S_{\text{\bf K}}$ across the VBS-to-W transition. Assuming $z=1$, in Fig.8(a) we examine the rescaled $S_{\text{\bf K}}L^{z+\eta-2}$ as a function of $J_{z2}/J_{\pm}$ at fixed $J_{z1}/J_{\pm}=4.5$, $\beta/L=1$, and $m_{z}=-1/6$. In Fig.8(b) we attempt at obtaining the the scaling function $F_{S}$ by plotting $S_{\text{\bf K}}L^{z+\eta-2}$ vs $(J_{z2}/J_{\pm}-J_{z2c}/J_{\pm})L^{1/\nu}$. In order to produce Fig.7 we have used $\eta=-0.30\pm 0.01$ and $\nu=0.50\pm 0.04$. Even though the data collapse may appear consistent with criticality, we argue that the absence of a clear crossing in Fig.7(a) and the unusual negative value of $\eta$ suggests that the transition is instead first order, in agreement with the expected three-dimensional three-state Potts model. Furthermore, we consider histograms of the order parameter $S_{\text{\bf K}}$, the total energy, and the kinetic energy at the critical point. We find that the histograms of the energy and order parameter are not bimodal for all system sizes accessible in our simulations (not shown). Instead, in Fig.9 we show histograms of the kinetic energy $K/J_{\pm}$ for system sizes $L=60$ (Fig.9(a)) and $L=66$ (Fig.9(b)) where we find that for the largest system we could simulate ($L=66$) the kinetic energy starts developing a two-peak histogram with a dominant beside a small, albeit statistically robust, peak suggesting the onset of phase coexistence. While this signal seems rather stable, definitive evidence in favor of a first-order phase transition requires verifying that the double-peak histograms remain stable for larger system sizes unavailable in our simulation setup.

VI Conclusion

Inspired by various molybdenum-based cluster magnets, Sheckelton et al. (2012); Mourigal et al. (2014); Sheckelton et al. (2014); Haraguchi et al. (2015); Chen et al. (2014, 2016) we have studied the phase diagram of a generic XXZ spin model on the anisotropic Kagome lattice using large-scale SSE quantum Monte Carlo simulations. We have found a remarkable tripartite entangled plaquette state surrounded by a valence-bond solid state and a ferromagnet state, and studied the respective transitions between these phases. We find that all the transitions toward the tripartite entangled plaquette state are conventional order-disorder transitions, either continuous (FM-W) or first order (VBS-W), which supports the idea that the W phase is a featureless symmetric state.

We have reexamined the FM-to-VBS transition using a canonical algorithm in the light of a recently proposed scaling analysis with two length scales.Shao et al. (2016) While our data seem consistent with that scenario, the large error bars associated with fitting numerical data to the predictions of the two-length scaling analysis prevented us from drawing a firm conclusion about the applicability of such scenario to our data. Furthermore, the scaling of the structure factor, which shows drifting crossings between system sizes, is perhaps consistent with the idea that this transition is first order as suggested in Ref.Isakov et al., 2006.

As we described in Sec. II, the XXZ spin model in our work describes the charge sector physics of the cluster magnet, and, as such, one expects a one-to-one mapping between the phases in our generic phase diagram (Fig.1) and the phases of the extended Hubbard model used in their description. In particular, the FM phase in our phase diagram corresponds to the Fermi liquid phase, the tripartite entangled plaquette state W corresponds to the cluster Mott insulator with one electron localized on every strong triangle, and the VBS phase corresponds to the plaquette charge ordered state.

Besides the relevance to cluster magnets, our results directly apply to other areas outside of condensed matter. In particular, strongly correlated bosonic atoms can be loaded onto highly tunable trimerized optical Kagome lattices.Santos et al. (2004); Damski et al. (2005); Jo et al. (2012) Such systems are realizations of the XXZ model discussed here via the mapping between spin-1/2 and hardcore boson models.Matsuda and Matsubara (1957)

Acknowledgements.

We would like to thank G. Baskaran, Y. Qi, S. Sachdev, M. Stoudenmire, Y. Wan, and W. Witczak-Krempa for enlightening discussions. We thank F. Hanington for a careful reading of the manuscript. We are especially indebted to A. Sandvik for discussions and suggestions related to the scaling forms with two divergent length scales. This research was supported by NSERC of Canada, the Perimeter Institute for Theoretical Physics, and the John Templeton Foundation. G.C. acknowledges support from the Thousand-Youth-Talent program of People’s Republic of China. R.G.M. acknowledges support from a Canada Research Chair. Research at Perimeter Institute is supported through Industry Canada and by the Province of Ontario through the Ministry of Research & Innovation. Numerical simulations were carried out on the Shared Hierarchical Academic Research Computing Network (SHARCNET).

appendix

In this section we briefly detail the procedure to measure chiral-chiral correlation functions in SSE quantum Monte Carlo simulations. Measurements of chiral-chiral correlation functions defined as $\langle E_{t}E_{t^{\prime}}\rangle$, where $E_{t}={\mathbf{S}_{i_{1}}}\left({\mathbf{S}_{i_{2}}}\times{\mathbf{S}_{i_{3}}}\right)$ ($i_{1},i_{2},i_{3}\in t$), are possible within SSE whenever the basic Hamiltonian breakup contains the triangles defined by the different indices ($i_{1},i_{2},i_{3}\in t$) in the correlators. Since the Hamiltonian breakup we have used in our study uses the underlying corner-sharing trianglesMelko (2007) defining the Kagome lattice, chiral-chiral correlations defined over those triangles can be naturally obtained within our simulations efficiently. By expanding the vector products in the correlator $\langle E_{t}E_{t^{\prime}}\rangle$, the SSE measurements are simplified to a combination of terms the form $\langle S^{z}_{i}S^{+}_{j}S^{-}_{k}S^{z}_{i^{\prime}}S^{+}_{j^{\prime}}S^{-}_{k^{\prime}}\rangle$ where $i,j,k\in t$ and $i^{\prime},j^{\prime},k^{\prime}\in t^{\prime}$. Terms of the form above can be measured through the estimator

[TABLE]

where the sum runs for ordered sub-sequences and the operators $S^{+}_{j}S^{-}_{k}$ appear at imaginary-time slice $t$ followed by a $S^{+}_{j^{\prime}}S^{-}_{k^{\prime}}$ at $t_{\text{next}}$ in the operator sequence $S_{n}$, and $n$ is the expansion order.Sandvik (1999); Melko (2007) We have implemented the estimator in Eq.6 for each of the terms appearing in the chiral-chiral correlation function. To test the validity of our approach, in Fig.10 we benchmark SSE estimates of one term, $\langle S^{z}_{i}S^{+}_{j}S^{-}_{k}S^{z}_{l}S^{+}_{m}S^{-}_{n}\rangle$, on a small $V=2\times 2\times 3$ cluster with periodic boundary conditions against exact diagonalization (ED) calculations. The sites $i,j,k$ and $l,m,n$ are depicted in the inset. The agreement of our SSE calculations with the ED results validates our approach.