Frustrated Magnetism of Dipolar Molecules on a Square Optical Lattice: Prediction of a Quantum Paramagnetic Ground State
Haiyuan Zou Erhai Zhao W. Vincent Liu







    𝐫𝐒𝐐Wilczek Quantum Center, School of Physics and Astronomy and T. D. Lee Institute,  Shanghai Jiao Tong University, Shanghai 200240, ChinaDepartment of Physics and Astronomy, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, USADepartment of Physics and Astronomy, George Mason University, Fairfax, Virginia 22030, USAWilczek Quantum Center, School of Physics and Astronomy and T. D. Lee Institute,  Shanghai Jiao Tong University, Shanghai 200240, ChinaDepartment of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260Center for Cold Atom Physics, Chinese Academy of Sciences, Wuhan 430071, China
        Motivated by the experimental realization of quantum spin models of polar molecule KRb in optical lattices,
        we analyze the spin 1/2 dipolar Heisenberg model with competing anisotropic, long-range exchange interactions.
        We show that, by tilting the orientation of dipoles using an external electric field, the dipolar spin system on square lattice 
        comes close to a maximally frustrated region similar, but not identical, to that of the J_1-J_2 model.
        This provides a simple yet powerful route to potentially realize a quantum spin liquid without the need for a triangular or kagome lattice. 
        The ground state phase diagrams obtained from Schwinger-boson and spin-wave theories consistently show
        a spin disordered region between the Néel, stripe, and spiral phase. 
        The existence of a finite quantum paramagnetic region is further confirmed by an unbiased variational ansatz
        based on tensor network states and a tensor renormalization group.
        Frustrated Magnetism of Dipolar Molecules on a Square Optical Lattice: Prediction of a Quantum Paramagnetic Ground State
            W. Vincent Liu
            Received: date / Accepted: date
        ========================================================================================================================
        
        
        
        Understanding highly entangled quantum matter remains a challenging goal of condensed matter physics <cit.>.
        One paradigmatic example is quantum spin liquids in frustrated spin systems which defy any conventional long range order characterized by broken symmetry at zero temperature <cit.>. Instead, the ground state features long-range entanglement and nonlocal excitations.
        
        Spin liquids are also fertile ground for studying quantum phases described by gauge field theories and topological order <cit.>.
        
        While the existence of spin liquids has been firmly established in a number of exactly solvable models, e.g., the toric code <cit.>
        or the honeycomb Kitaev model <cit.>, the nature of the ground states for many frustrated spin models, e.g., the Heisenberg model on kagome 
        lattices or the J_1-J_2 model on square lattices, still remains controversial despite the great theoretical progress in recent years <cit.>.  
        
        An unambiguous experimental identification of quantum spin liquids in solid state materials also seems elusive <cit.>. 
        It is, then, important to explore new physical systems that can cleanly realize
        well-defined spin models which have potential spin liquid ground states.
        
        Recent breakthrough experiments on magnetic atoms <cit.> and polar molecules  <cit.> confined in deep optical lattices introduced 
        a new class of lattice spin models with competing exchange interactions that are long-ranged and anisotropic. 
        The resulting spin Hamiltonians, such as the dipolar XXZ
        and dipolar Heisenberg models, are highly tunable by the external fields that couple to the magnetic and electric dipoles <cit.>. 
        
        Here, we show that these models on square lattices
        feature strong exchange (not geometric) frustration and a quantum paramagnetic ground state for intermediate dipole tilting angles. This claim is consistently supported by physical arguments, two independent semiclassical analytical methods, and full numerical calculation based on tensor network ansatz <cit.>.
        
        Our key insight is that spin liquids may arise naturally from the system of tilted, interacting dipoles on square lattices, without the
        requirement of peculiar (e.g., triangular or kagome) lattices or exotic (e.g., Kitaev or ring-exchange) interactions. 
        
         The dipolar XXZ and Heisenberg model.—First, we define the dipolar XXZ model on a square optical lattice,
         
        H_XXZ=J/2∑_i≠ j f(𝐫_i-𝐫_j) (S^x_i S^x_j+S^y_i S^y_j+η S^z_i S^z_j).
        
        Here i and j label the lattice sites, 𝐒_i=(S^x_i,S^y_i,S^z_i) are the spin (or pseudospin) operators at site i, and η is the exchange anisotropy.
        The key new feature here is that the coupling between the two spins depends on their relative position 𝐫=𝐫_i-𝐫_j
        and the external field (dipole) direction d̂
        f(𝐫)=[1-3(r̂·d̂)^2](a/r)^3,
        
        with a the lattice constant [Fig. 1(a)]. This geometric factor, characteristic of the dipole-dipole interaction, dictates that spin interactions are long-ranged and anisotropic. For the special case of η=1, H_XXZ reduces to the dipolar Heisenberg model
        
        H_d=J/2∑_i≠ j f(𝐫_i-𝐫_j) 𝐒_i ·𝐒_j,
        
        and for η=0, it reduces to the dipolar XY model, H_XY.
        
        Spin models of the form of H_XXZ have been realized experimentally in two settings. In Ref. <cit.>, the spin dynamics of a gas of ^52Cr atoms in optical lattices was observed. Each Cr atom carries a magnetic moment of 7μ_B and hyperfine spin S=3. An external magnetic field is used to align the magnetic dipoles in the direction of d̂. Such a dipolar gas of Cr in a deep lattice is shown to be described by H_XXZ with J=-μ_0(gμ_B)^2/4π a^3<0 and η=-2<cit.>. 
        
        Note that J induced by the dipolar interaction is, contrary to the superexchange, independent of the tunneling, and it can be set as the unit of energy. 
         
        
        Polar molecules such as ^40K^87Rb confined in optical lattices with negligible tunneling provide another way to realize H_XXZ with S=1/2 and tunable J and η<cit.>. Each molecule carries an electric dipole moment 𝐝 and undergoes rotation with angular momentum 𝐉 [see Fig. 1(a)].  
        Here, the pseudospin 1/2 refers to two rotational states of the molecule labeled by |j,m⟩, where j is the quantum number of the rotational angular momentum 𝐉 and m is its projection onto the quantization axis, chosen as the direction of the external electric field E. More details can be found in Ref. <cit.>.
        
        The dipole-dipole interaction projected onto the sub-Hilbert space of the pseudospins then takes the form of a spin Hamiltonian, where the spin flips correspond to transitions between the rotational states.
        
        For example, by choosing |j,m⟩=|0,0⟩ and |1,0⟩ as the pseudospin down and up respectively, Refs. <cit.> showed that the system is described by the effective Hamiltonian H_XXZ with
        J=D_t ^2/2πϵ_0 a^3>0 and η=(D_1-D_0)^2/2D_t^2>0. 
        Here the dipole matrix element D_t=⟨1,0|d^0|0,0⟩, D_1=⟨1,0|d^0|1,0⟩, D_0=⟨0,0|d^0|0,0⟩, and d^0 together with d^± form the vector dipole operator in the spherical basis <cit.>. 
        
        
        The anisotropy η increases monotonically with E. As shown in Ref. <cit.>, 
        when E≃ 1.7B/|𝐝| with B the energy splitting of the two pseudospin states, η=1, and one arrives at the dipolar Heisenberg model H_d. In the KRb experiment <cit.> carried out at zero field and cubic lattice, η→ 0, the dipolar XY model H_XY was realized with  J on the order of 100 Hz. Despite the low filling factor and high entropy, coherent spin dynamics was observed via Ramsey spectroscopy <cit.> and modeled theoretically in Ref. <cit.>. Recently Yao  et al. <cit.> considered general η and worked out the phase diagram of H_XXZ on the Kagome and triangular lattice using Density Matrix Renormalization Group (DMRG).. For both lattices, they found evidence for quantum spin liquid centering around the Heisenberg limit, η=1 and θ=0, in which θ is defined by d̂·x̂=sinθcosϕ with x̂ representing a base vector of the square lattice. Thus the physics is connected to a geometrically frustrated Heisenberg model on both lattices, with additional longer range interactions and anisotropy η.   
        
        In this Letter, we study the phases of H_d on a square lattice as the dipoles are tilted towards the lattice plane [see Fig. 1(a)] for S=1/2 and J>0. We show that strong frustration occurs at intermediate dipole tilting angle θ, 
        leading to a quantum paramagnetic ground state. We emphasize that, here, the frustration is not imposed by the lattice geometry, but instead, is due to the competition between the exchange interactions, analogous to the J_1-J_2 model. Relatedly, the quantum paramagnetic phase appears at intermediate θ values (not around θ=0 as in Ref. <cit.>) between the Néel and the stripe orders. Thus, it differs qualitatively from the spin liquids studied in Ref. <cit.>. We will also employ different methods to solve the dipolar quantum spin models. 
        
        Competing exchanges for tilted dipoles.—To appreciate the possible phases of H_d as d̂ is tuned as well as its connection to frustrated quantum spin models <cit.>, let us consider the leading exchange couplings between the nearest neighbors, J_x=Jf(ax̂) and J_y=Jf(aŷ), and the next nearest neighbors, J_d=Jf(ax̂+aŷ) and J'_d=Jf(ax̂-aŷ) [Fig. 1(b)]. Their relative magnitudes and signs depend sensitively on the dipole tilting angle θ and ϕ. One example is shown in Fig. 1(b) for fixed ϕ=25^∘. At small θ, J_x∼ J_y dominates because it is about three times that of J_d∼ J'_d. The situation is reminiscent of the J_1-J_2 model in the regime of the Néel order.  As θ is increased, J_d and J'_d grow relative to J_x and J_y. The system becomes more frustrated due to the increased competition of the exchanges. This is the most interesting parameter region. Around θ≃ 40^∘, J_x and J_d vanish while J'_d∼ 0.4 J_y. The model can be viewed as coupled Heisenberg chains. For even larger θ, J_x and J_d switch signs to become ferromagnetic, and the stripe order is expected. Clearly, the physics of H_d is much richer than the J_1-J_2 model.  In fact, the two models only overlap at one single point, θ=ϕ=0, where J_2/J_1=1/2√(2)≈ 0.35 and the system is Néel ordered. 
        
        The degree of frustration can be measured by the “spin gap" Δ, the energy difference between the ground and the first excited state, from exact diagonalization of H_d for a 4× 4 lattice <cit.>. For example, we observe a pronounced peak in Δ around θ∼ 28^∘ for ϕ=25^∘, which indicates strong frustration and points to a gapped, spin disordered ground state <cit.>. For fixed ϕ=35^∘, the spin structure factor shows a clear peak at (π,π) for θ∼ 15^∘ for the Néel order, a peak at (0,π) for θ∼ 50^∘ for the stripe order, but no well defined peaks around θ∼ 35^∘, consistent with the argument above. 
        
        
        
        
        Spin-wave and Schwinger-boson theory.—First, we obtain a coarse phase diagram of H_d on the (θ,ϕ) plane using two widely adopted analytical methods in frustrated quantum magnetism. This will help identify the interesting regions for the more expensive tensor network calculations to focus on. The starting point is the classical solution of H_d by the Luttinger-Tisza method <cit.>. H_d is of the form ∑_ijJ_ij𝐒_i·𝐒_j with hard spin constraint 𝐒_i=S and J_ij only depends on 𝐫_i-𝐫_j. A theorem states that the classical ground state is a planar spin spiral, _/S=x̂cos(·)+ŷsin(·) with an ordering wave vector = (Q_x,Q_y)<cit.>. The classical phase diagram <cit.> consists of three phases. The first is the Néel order corresponding to =(π,π) for small θ. The second is the stripe phase with =(0,π) for large θ but not too large ϕ. These two spin orders are collinear. The third, spiral phase fills the rest of the phase diagram, for large θ and ϕ, where  varies continuously and, in general, is incommensurate with the lattice.
        
        Beyond the classical limit, quantum fluctuations will suppress the magnetic order and shift the phase boundary. These effects can be described qualitatively by modified spin wave theory  <cit.>.  In the Holstein-Primakoff representation, we expand H_d in a series of 1/S and keep up to the quartic order of bosonic operators, i.e., we take into account the interactions between the linear spin waves. The bosonic Hamiltonian is solved by self-consistent mean field theory <cit.>. The result is summarized in Fig. 2(a). We find that the phase boundary of the Néel (stripe) phase moves towards smaller (larger) θ values, opening up an intermediate region in between where the magnetization vanishes. The spiral phase also recedes to higher ϕ values. We label this quantum paramagnetic region with QP. This is precisely the region where the various exchanges compete and the system is most frustrated. 
        
        Alternatively, we can take into account quantum fluctuations by the rotationally invariant Schwinger boson mean field theory which is nonperturbative in S<cit.>. It is a well tested method capable of describing both magnetically ordered and spin liquid states of frustrated spin models <cit.>. 
        The resulting phase diagram is shown in Fig. 2(b). Here, each magnetic order corresponds to condensation of bosons at a certain wave vector . Within a finite strip region labeled by QP between the Néel and stripe phase, the condensation fraction vanishes and the spin excitations are gapped, corresponding to a quantum paramagnetic phase. The fact that two different approximations agree on the existence of QP indicates that it must be a robust feature of the model H_d.
        
        Phase diagram from a tensor network ansatz.—A variational ansatz based on tensor network states <cit.> has recently emerged
        as an accurate and unbiased algorithm for solving two dimensional frustrated quantum spin models <cit.>. In this approach, the ground state many-body wave function |Ψ⟩ is constructed from a network of tensors T_i defined on lattice site i: |Ψ⟩=∏_i T_i, where  stands for contraction of neighboring tensors. Each tensor T_i has four virtual legs (indices), each with bond dimension D designed to build up the quantum entanglement between lattice sites, and one physical leg representing the spin. We choose a L× L cluster as the unit cell with periodic boundary conditions. The algorithm starts with L^2 random tensors, and imaginary time evolution is used to update the local tensors, |ψ'⟩=exp(-τ H)|ψ⟩, until convergence is achieved. We adopt the simple update scheme <cit.> based on singular value decomposition. By using the Trotter-Suzuki formula exp(-τ H)≈∏_i=1^4exp(-τ H_i)+O(τ^2), each iteration of projection for one plaquette can be done using exp(-τ H_i)(i=1,2,3,4) in four separate steps, in which each step evolves three sites (a right triangle) in one plaquette with H_i contains only three terms of the Hamiltonian. For example, H_1,2 contains J_x, J_y, and J_d terms and H_3,4 contains J_x, J_y and J'_d terms (See Refs. <cit.>).
        
        The expectation value of a local operator O_j at site j, ⟨ O_j⟩=⟨Ψ|O_j|Ψ⟩/⟨Ψ|Ψ⟩, can be computed by tensor contraction, ⟨ O_j⟩= (𝒪_j∏_i≠ j𝒯_i)/∏_i 𝒯_i where 𝒯_i=T_i^† T_i and 𝒪_j=T_j^† O_j T_j. We evaluate it using an iterative, real space coarse-graining procedure known as the tensor renormalization group which enables one to reach the thermodynamic limit <cit.>. In this way, we calculate the order parameters such as magnetization M=√(⟨ S_x⟩^2+⟨ S_y⟩^2+⟨ S_z⟩^2) <cit.>.
        
        
         
        
        With increasing D, quantum fluctuations beyond spin wave or Schwinger boson analysis are taken into account. The suppression of M is illustrated in Fig. 3(a) for different D values at fixed ϕ=15^∘. By extrapolating the results to infinite D, we can determine the phase boundary of the Néel and stripe phases. Repeating the procedure for different ϕ values, we obtain the phase diagram Fig. 3(b). It firmly establishes the existence of a finite quantum paramagnetic region (in red), about one degree wide in θ and persisting from ϕ=0 up to ϕ=20^∘, where the magnetization is completely suppressed to zero. The paramagnetic phase is narrower than the prediction of the Schwinger boson mean field theory which tends to overestimate the spin disordered region. Inside the Néel phase, there is a sudden drop of M.
        
        Note that the spiral phase, in general, is incompatible with the L× L cluster choice, even for large L. So we refrain from carrying out the tensor network ansatz beyond ϕ=20^∘. On the other hand, our numerics indicates that the phase boundary presented in Fig. 3(b) is not expected to depend sensitively on L as it varies <cit.>. Finally, we point out that the quantum paramagnetic phase is a robust feature of the dipolar XXZ model. It persists when η is tuned away from the Heisenberg limit, e.g., down to η=0.5 <cit.>.
        
        
        It is challenging to pin down the precise nature of the paramagnetic phase found here in the dipolar Heisenberg model. Similar difficulties also arise for the J_1-J_2 model where the latest DMRG result <cit.> suggests that the paramagnetic region may consist of a subregion with a plaquette valence bond solid (VBS) order and a second, spin liquid or quantum critical region. Possible spin liquid states for the J_1-J_2 model on square lattices have been classified within the framework of the Schwinger boson mean field theory <cit.>. Yet it remains unclear which one is realized in the ground state. It is possible that the QP region of H_d may contain some VBS order. Unlike the J_1-J_2 model, the C_4 rotation symmetry is broken in H_d as soon as the dipoles are tilted, which may disfavor the plaquette VBS. Because of the limitation of the cluster size, we could not accurately compute the dimer correlation functions.
        Future numerical work with larger L and D is required to shed light on this open issue.
        The new formulation of symmetric tensor networks <cit.> and Lanczos iteration <cit.> seems promising to detect the possible topological order and accessing the excitation spectrum.
        
        In summary, we presented consistent evidence that a quantum paramagnetic phase emerges from the simple physical system of interacting, tilted dipoles confined on square optical lattices. Our analysis of the dipolar Heisenberg model for general (θ,ϕ) adds a new dimension to frustrated quantum magnetism. It allows the exploration of potential spin liquids beyond the J_1-J_2 model which has not been realized cleanly so far. For KRb, J is about 100 Hz, or 5 nK, similar to the superexchange scale t^2/U of the Fermi Hubbard model recently studied using quantum gas microscope <cit.>. Thus, it seems possible to probe the spin order or spin correlations of H_d and related models in future experiments.
        
         
        We thank Ying-Jer Kao, Bo Liu, Jaime Merino, and Ling Wang for helpful discussions. This work is supported by the U.S. AFOSR Grant No. FA9550-16-1-0006 (H.Z., E.Z., and W.V.L.), NSF PHY-1205504 (E.Z.), and ARO Grant No. W911NF-11-1-0230, the Chinese National Science Foundation through the Overseas Scholar Collaborative Program (Grant No. 11429402) sponsored by Peking University and another grant (No.11227803), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB21010100) (W.V.L.). Part of the simulations were done at the supercomputing resources provided by the University of Pittsburgh Center for Simulation and Modeling.
        
        
        
        
        
        
        
        
        
        
        50
        fxundefined [1]
         ifx#1
        fnum [1]
          #1firstoftwo
         secondoftwo
         
        fx [1]
          #1firstoftwo
         secondoftwo
         
        
        
        
        
        
        noop [0]secondoftwo
        
        ref[1]@startlink#1@href
        href[1]#1@endlink
        anitize@url [0] `
        12 `$12
          `&12 `#12 `1̂2 `_12 `%12
        startlink[1]
        endlink[0]
        
        rl [1]href #1
        
        
        
        
        
        
        
        
        
        
        
        
        @bib@innerbibempty
        
         [Savary and Balents(2017)]newRev2
          authorauthorL. Savary and authorL. Balents, http://stacks.iop.org/0034-4885/80/i=1/a=016502journaljournalReports on Progress in Physicsvolume80, pages016502 (year2017)NoStop
         [Zhou et al.(2017)Zhou,
          Kanoda, and Ng]newRev
          authorauthorY. Zhou, authorK. Kanoda,  andauthorT.-K. Ng,  10.1103/RevModPhys.89.025003journaljournalRev. Mod. Phys. volume89, pages025003 (year2017)NoStop
         [Diep(2013)]diep2013frustrated
          authorauthorH. Diep, @noop titleFrustrated spin
          systems (publisherWorld Scientific, year2013)NoStop
         [Wen(2002)]PhysRevB.65.165113
          authorauthorX.-G. Wen,  10.1103/PhysRevB.65.165113journaljournalPhys. Rev. B volume65,pages165113 (year2002)NoStop
         [Kitaev(2003)]Kitaev20032
          authorauthorA. Kitaev,  https://doi.org/10.1016/S0003-4916(02)00018-0journaljournalAnnals of Physics volume303, pages2  (year2003)NoStop
         [Kitaev(2006)]Kitaev20062
          authorauthorA. Kitaev,  https://doi.org/10.1016/j.aop.2005.10.005journaljournalAnnals of Physics volume321, pages2  (year2006)NoStop
         [Yan et al.(2011)Yan,
          Huse, and White]Yan1173
          authorauthorS. Yan, authorD. A. Huse,   and authorS. R. White,  10.1126/science.1201080journaljournalScience volume332, pages1173
          (year2011)NoStop
         [Iqbal et al.(2013)Iqbal,
          Becca, Sorella, and Poilblanc]Gapless2013
          authorauthorY. Iqbal, authorF. Becca,
          authorS. Sorella,  andauthorD. Poilblanc,  10.1103/PhysRevB.87.060405journaljournalPhys. Rev. B volume87, pages060405 (year2013)NoStop
         [Jiang et al.(2012)Jiang,
          Yao, and Balents]J1J2Jiang
          authorauthorH.-C. Jiang, authorH. Yao,  andauthorL. Balents,  10.1103/PhysRevB.86.024424journaljournalPhys. Rev. B volume86, pages024424 (year2012)NoStop
         [Gong et al.(2014)Gong,
          Zhu, Sheng, Motrunich, andFisher]J1J2Plaq
          authorauthorS.-S. Gong, authorW. Zhu, authorD. N. Sheng, authorO. I. Motrunich,  and authorM. P. A. Fisher, 
          10.1103/PhysRevLett.113.027201journaljournalPhys. Rev. Lett. volume113, pages027201 (year2014)NoStop
         [Wang et al.(2016)Wang,
          Gu, Verstraete, and Wen]tnj1j2
          authorauthorL. Wang, authorZ.-C. Gu,
          authorF. Verstraete,  andauthorX.-G. Wen,  10.1103/PhysRevB.94.075143journaljournalPhys. Rev. B volume94, pages075143 (year2016)NoStop
         [de Paz et al.(2013)de Paz,
          Sharma, Chotia, Maréchal,
          Huckans, Pedri, Santos,
          Gorceix, Vernac, and Laburthe-Tolra]PhysRevLett.111.185305
          authorauthorA. de Paz, authorA. Sharma,
          authorA. Chotia, authorE. Maréchal, authorJ. H. Huckans, authorP. Pedri, authorL. Santos, authorO. Gorceix, authorL. Vernac,  and authorB. Laburthe-Tolra, 
          10.1103/PhysRevLett.111.185305journaljournalPhys. Rev. Lett. volume111, pages185305 (year2013)NoStop
         [Yan et al.(2013)Yan,
          Moses, Gadway, Covey,
          Hazzard, Rey, Jin, andYe]Yan:2013xe
          authorauthorB. Yan, authorS. A. Moses,
          authorB. Gadway, authorJ. P. Covey, authorK. R. A. Hazzard, authorA. M. Rey, authorD. S. Jin,  and authorJ. Ye, http://dx.doi.org/10.1038/nature12483journaljournalNature volume501, pages521 (year2013)NoStop
         [Hazzard et al.(2014)Hazzard, Gadway, Foss-Feig, Yan, Moses, Covey, Yao,
          Lukin, Ye, Jin, andRey]PhysRevLett.113.195302
          authorauthorK. R. A.Hazzard, authorB. Gadway, authorM. Foss-Feig,
          authorB. Yan, authorS. A. Moses, authorJ. P. Covey, authorN. Y. Yao, authorM. D. Lukin, authorJ. Ye, authorD. S.Jin,  and authorA. M.Rey,  10.1103/PhysRevLett.113.195302journaljournalPhys. Rev. Lett. volume113, pages195302 (year2014)NoStop
         [Gorshkov et al.(2011a)Gorshkov, Manmana, Chen, Ye, Demler,
          Lukin, and Rey]PhysRevLett.107.115301
          authorauthorA. V. Gorshkov, authorS. R. Manmana, authorG. Chen,
          authorJ. Ye, authorE. Demler, authorM. D. Lukin,  and authorA. M. Rey, 
          10.1103/PhysRevLett.107.115301journaljournalPhys. Rev. Lett. volume107, pages115301 (year2011a)NoStop
         [Yao et al.(2015)Yao,
          Zaletel, Stamper-Kurn, and Vishwanath]DSL2015
          authorauthorN. Y. Yao, authorM. P. Zaletel,
          authorD. M. Stamper-Kurn,   and authorA. Vishwanath,@noop  (year2015), http://arxiv.org/abs/arXiv:1510.06403arXiv:1510.06403NoStop
         [Verstraete and Cirac(2004)]PEPS1
          authorauthorF. Verstraete and authorJ. I. Cirac, @noop  (year2004), http://arxiv.org/abs/arXiv:cond-mat/0407066arXiv:cond-mat/0407066NoStop
         [Jordan et al.(2008)Jordan,
          Orús, Vidal, Verstraete,  and Cirac]iPEPS
          authorauthorJ. Jordan, authorR. Orús,
          authorG. Vidal, authorF. Verstraete,  and authorJ. I. Cirac, 
          10.1103/PhysRevLett.101.250602journaljournalPhys. Rev. Lett. volume101, pages250602 (year2008)NoStop
         [Corboz et al.(2010a)Corboz, Orús,
          Bauer, and Vidal]iPEPS2
          authorauthorP. Corboz, authorR. Orús,
          authorB. Bauer,  and authorG. Vidal, 
          10.1103/PhysRevB.81.165104journaljournalPhys.
          Rev. B volume81, pages165104
          (year2010a)NoStop
         [Levin and Nave(2007)]LevinTRG
          authorauthorM. Levin and authorC. P. Nave,  10.1103/PhysRevLett.99.120601journaljournalPhys. Rev. Lett. volume99, pages120601 (year2007)NoStop
         [Xie et al.(2012)Xie,
          Chen, Qin, Zhu, Yang, and Xiang]HOTRG1
          authorauthorZ. Y. Xie, authorJ. Chen, authorM. P. Qin, authorJ. W. Zhu, authorL. P. Yang,  and authorT. Xiang,  10.1103/PhysRevB.86.045139journaljournalPhys. Rev. B volume86, pages045139 (year2012)NoStop
         [Gorshkov et al.(2011b)Gorshkov, Manmana, Chen, Demler, Lukin, and Rey]PhysRevA.84.033619
          authorauthorA. V. Gorshkov, authorS. R. Manmana, authorG. Chen,
          authorE. Demler, authorM. D. Lukin,  and authorA. M. Rey, 
          10.1103/PhysRevA.84.033619journaljournalPhys.
          Rev. A volume84, pages033619
          (year2011b)NoStop
         [Balents(2010)]Balents:2010kq
          authorauthorL. Balents, http://dx.doi.org/10.1038/nature08917journaljournalNature volume464,pages199 (year2010)NoStop
         [not()]noteS
          @noop noteSee Supplemental Material for a detailed
          descriptionNoStop
         [Sindzingre(2004)]refId0
          authorauthorP. Sindzingre,  10.1103/PhysRevB.69.094418journaljournalPhys. Rev. B volume69, pages094418 (year2004)NoStop
         [Luttinger and Tisza(1946)]LuttingerTisza
          authorauthorJ. M. Luttinger and authorL. Tisza,  10.1103/PhysRev.70.954journaljournalPhys. Rev. volume70,pages954 (year1946)NoStop
         [Kaplan and Menyuk(2007)]kaplan2007spin
          authorauthorT. A. Kaplan and authorN. Menyuk,  10.1080/14786430601080229journaljournalPhilosophical Magazine volume87, pages3711 (year2007)NoStop
         [Takahashi(1989)]MSW
          authorauthorM. Takahashi,  10.1103/PhysRevB.40.2494journaljournalPhys. Rev. B volume40, pages2494 (year1989)NoStop
         [Chu and Shen(1991)]SelfCSW
          authorauthorD. Chu and authorJ.-l. Shen, 10.1103/PhysRevB.44.4689journaljournalPhys. Rev. B volume44, pages4689 (year1991)NoStop
         [Cui(1992)]SpinW1
          authorauthorS.-M. Cui, http://stacks.iop.org/0953-8984/4/i=28/a=003journaljournalJournal of Physics: Condensed Mattervolume4, pagesL389 (year1992)NoStop
         [Arovas and Auerbach(1988)]PhysRevB.38.316
          authorauthorD. P. Arovas and authorA. Auerbach,  10.1103/PhysRevB.38.316journaljournalPhys. Rev. B volume38,pages316 (year1988)NoStop
         [Read and Sachdev(1991)]PhysRevLett.66.1773
          authorauthorN. Read and authorS. Sachdev, 10.1103/PhysRevLett.66.1773journaljournalPhys. Rev. Lett. volume66, pages1773 (year1991)NoStop
         [Mezio et al.(2013)Mezio,
          Sposetti, Manuel, and Trumper]mezio2013broken
          authorauthorA. Mezio, authorC. N. Sposetti,
          authorL. O. Manuel,  andauthorA. E. Trumper, http://stacks.iop.org/0953-8984/25/i=46/a=465602journaljournalJournal of Physics: Condensed Matter volume25, pages465602 (year2013)NoStop
         [Merino et al.(2014)Merino,
          Holt, and Powell]merino2014spin
          authorauthorJ. Merino, authorM. Holt,  andauthorB. J. Powell,  10.1103/PhysRevB.89.245112journaljournalPhys. Rev. B volume89, pages245112 (year2014)NoStop
         [Lee et al.(2014)Lee,
          Jeong, Hwang, and Kim]lee2014emergent
          authorauthorS. Lee, authorJ.-S. Jeong,
          authorK. Hwang,  and authorY. B. Kim, 
          10.1103/PhysRevB.90.134425journaljournalPhys.
          Rev. B volume90, pages134425
          (year2014)NoStop
         [Yang and Wang(2016)]SchwingerJ1J2
          authorauthorX. Yang and authorF. Wang, 10.1103/PhysRevB.94.035160journaljournalPhys. Rev. B volume94, pages035160 (year2016)NoStop
         [Zou et al.(2016)Zou,
          Liu, Zhao, and Liu]compass120
          authorauthorH. Zou, authorB. Liu, authorE. Zhao,  and authorW. V. Liu, http://stacks.iop.org/1367-2630/18/i=5/a=053040journaljournalNew Journal of Physics volume18, pages053040 (year2016)NoStop
         [Jiang et al.(2016)Jiang,
          Kim, Han, and Ran]HKagome1
          authorauthorS. Jiang, authorP. Kim, authorJ. H. Han,  and authorY. Ran, @noop 
          (year2016), http://arxiv.org/abs/arXiv:1610.02024arXiv:1610.02024NoStop
         [Liao et al.(2017)Liao,
          Xie, Chen, Liu, Xie, Huang, Normand, and Xiang]HKagome2
          authorauthorH. J. Liao, authorZ. Y. Xie,
          authorJ. Chen, authorZ. Y. Liu, authorH. D. Xie, authorR. Z.Huang, authorB. Normand,  and authorT. Xiang,  10.1103/PhysRevLett.118.137202journaljournalPhys. Rev. Lett. volume118, pages137202 (year2017)NoStop
         [Jiang et al.(2008)Jiang,
          Weng, and Xiang]simpleup
          authorauthorH. C. Jiang, authorZ. Y. Weng,   and authorT. Xiang,  10.1103/PhysRevLett.101.090603journaljournalPhys. Rev. Lett. volume101, pages090603 (year2008)NoStop
         [Corboz et al.(2010b)Corboz, Jordan,  and Vidal]nnnsimple
          authorauthorP. Corboz, authorJ. Jordan,   and authorG. Vidal,  10.1103/PhysRevB.82.245119journaljournalPhys. Rev. B volume82, pages245119 (year2010b)NoStop
         [Corboz and Mila(2013)]tnSSmodel
          authorauthorP. Corboz and authorF. Mila, 10.1103/PhysRevB.87.115144journaljournalPhys. Rev. B volume87, pages115144 (year2013)NoStop
         [Jiang and Ran(2015)]SymmTN
          authorauthorS. Jiang and authorY. Ran, 10.1103/PhysRevB.92.104414journaljournalPhys. Rev. B volume92, pages104414 (year2015)NoStop
         [Mei et al.(2017)Mei,
          Chen, He, and Wen]HKagomeZ2SL
          authorauthorJ.-W. Mei, authorJ.-Y. Chen,
          authorH. He,  and authorX.-G. Wen, 
          10.1103/PhysRevB.95.235107journaljournalPhys.
          Rev. B volume95, pages235107
          (year2017)NoStop
         [Huang et al.(2016)Huang,
          Liao, Liu, Xie, Xie, Zhao, Chen, and Xiang]LanczosTN
          authorauthorR.-Z. Huang, authorH.-J. Liao,
          authorZ.-Y. Liu, authorH.-D. Xie, authorZ.-Y. Xie, authorH.-H. Zhao, authorJ. Chen,  and authorT. Xiang, @noop  (year2016),http://arxiv.org/abs/arXiv:1611.09574arXiv:1611.09574NoStop
         [Parsons et al.(2015)Parsons, Huber, Mazurenko, Chiu, Setiawan, Wooley-Brown, Blatt, and Greiner]Parsons1253
          authorauthorM. F. Parsons, authorF. Huber,
          authorA. Mazurenko, authorC. S. Chiu, authorW. Setiawan, authorK. Wooley-Brown, authorS. Blatt,  and authorM. Greiner, 
          10.1103/PhysRevLett.114.213002journaljournalPhys. Rev. Lett. volume114, pages213002 (year2015)NoStop
         [Cheuk et al.(2015)Cheuk,
          Nichols, Okan, Gersdorf,
          Ramasesh, Bakr, Lompe, andZwierlein]Cheuk1260
          authorauthorL. W. Cheuk, authorM. A. Nichols,
          authorM. Okan, authorT. Gersdorf, authorV. V. Ramasesh, authorW. S. Bakr, authorT. Lompe,  and authorM. W. Zwierlein, 
          10.1103/PhysRevLett.114.193001journaljournalPhys. Rev. Lett. volume114, pages193001 (year2015)NoStop
         [Parsons et al.(2016)Parsons, Mazurenko, Chiu, Ji, Greif, and Greiner]Parsons2016
          authorauthorM. F. Parsons, authorA. Mazurenko,
          authorC. S. Chiu, authorG. Ji, authorD. Greif,  and authorM. Greiner, 
          10.1126/science.aag1430journaljournalScience volume353, pages1253
          (year2016)NoStop
         [Boll et al.(2016)Boll,
          Hilker, Salomon, Omran,
          Nespolo, Pollet, Bloch, andGross]Boll1257
          authorauthorM. Boll, authorT. A. Hilker,
          authorG. Salomon, authorA. Omran, authorJ. Nespolo, authorL. Pollet, authorI. Bloch,  and authorC. Gross,  10.1126/science.aag1635journaljournalScience volume353, pages1257 (year2016)NoStop
         [Cheuk et al.(2016)Cheuk,
          Nichols, Lawrence, Okan,
          Zhang, Khatami, Trivedi,
          Paiva, Rigol, and Zwierlein]Cheuk2016
          authorauthorL. W. Cheuk, authorM. A. Nichols,
          authorK. R. Lawrence, authorM. Okan, authorH. Zhang, authorE. Khatami, authorN. Trivedi, authorT. Paiva, authorM. Rigol,  and authorM. W.Zwierlein,  10.1126/science.aag3349journaljournalScience volume353, pages1260 (year2016)NoStop
        
        
         Supplemental Materials for “Frustrated magnetism of dipolar molecules on a square optical lattice: Prediction of a quantum paramagnetic ground state"
        
        Haiyuan Zou, Erhai Zhao and W. Vincent Liu
        
        
        § CLASSICAL PHASE DIAGRAM
        
         The classical ground state of a translationally invariant spin model ∑_ijJ_ij_i·_j
         on a Bravais lattice can be obtained by minimizing the energy within the planar helix ansatz
        _/S=x̂cos(·)+ŷsin(·),
        
        where x̂ and ŷ form an orthonormal basis and = (Q_x,Q_y) is the ordering wavevector <cit.>.
        
        This variational ansatz satisfies the hard spin constraint |_|=S. The classical energy of the dipolar 
        Heisenberg model depends on  via
        
        2ℋ_cl/N S^2   =   J_xcos(Q_x)+J_ycos(Q_y)
           +   J_d'cos(Q_x+Q_y)+J_dcos(Q_x-Q_y).
        
         The set of wavevectors  minimizing the classical energy will be denoted as {}.  
         We compare the energies of the incommensurate spiral {_I}, the stripe _s = (0, π) and the Néel _n = (π, π) order. The result is the classical phase diagram shown in Fig. <ref>.
        Upon crossing the phase boundary, e.g., from the Neel (or the stripe) phase to the incommensurate spiral phase, the wavevector Q varies continuously. 
        
        For example, in the special case of ϕ=45^∘ (see the inset of Fig. <ref>), Q_x=Q_y≡ Q, where Q changes continuously from π on the phase boundary between Neel and spiral phase to π/4 at the upper right corner of the (θ, ϕ) diagram.
        
        
        § EXACT DIAGONALIZATION
        
        We calculate the “spin gap" Δ, the energy defference between the ground and the first excited state from exact diagonalization of the dipolar Heisenberg Hamiltonian for a 4× 4 lattice. Fig. <ref> shows Δ as functions of θ for different ϕ. For ϕ≤ 30^∘, the pronounced peak at each line indicates strong frustration. For ϕ=35^∘, The appearance of the second peak corresponds to the transition from the stripe phase to the incommensurate spiral phase as θ increases. For ϕ≥ 40^∘, the disappearance of the first peak indicates the transition from the Neel phase to the spiral phase. 
        
        
        § SCHWINGER BOSON THEORY
        
        We outline the Schwinger boson mean field theory (SBMF) of the dipolar Heisenberg model. The starting point is the bosonic representation of the spin operators
        
        S^+=a^† b,
        
        S^-=b^† a,
        
        S_z=1/2(a^† a- b^† b).
        
        with the constraint
        1/2(a^† a + b^† b)=S. 
        For the square lattice, introduce the antiferromagnetic (A) and ferromagnetic (B) bond operators
        
        A_ij=1/2(a_i b_j - b_i a_j),
        
        B_ij=1/2(a^†_i a_j + b^†_i b_j).
         
        In terms of the bond operators, the spin exchange term becomes
        _i ·_j=:B^†_ij B_ij:-A^†_ij A_ij,
        
        where :: means normal order of bosonic operators. Note that B and A are related by operator identity :B^†_ij B_ij:+A^†_ij A_ij=S^2. 
        We adopt the rotational invariant formulation of SBMF and perform mean field decoupling for both A and B,
        _i ·_j≃ [β_ij B_ij-α_ijA_ij + h.c. ]-|β_ij|^2+|α_ij|^2,
         
         where
         α_ij=⟨ A^†_ij⟩,
        β_ij=⟨ B^†_ij⟩,
         
        and ⟨ ... ⟩ denotes the ground state expectation value. This is known to perform better in describing the phases of frustrated spin systems compared to antatz that only keep either A or B. Also, within SBMF, the constraint Eq. (<ref>) is only enforced on average by introducing the Lagrange multiplier λ. The mean field Hamiltonian then takes the form
         ℋ_MF   =   1/2∑_i≠ jJ_ij [β_ij B_ij-α_ijA_ij+ h.c.- |β_ij|^2+|α_ij|^2 ] 
           +   λ∑_i[a^†_i a_i +
        b^†_i b_i-2S].
        
        
        Long range magnetic order corresponds to condensation of the a and/or b bosons. To treat the condensate fraction, we decompose each operator into
        
        a_i = ã_i + x_i,   x_i=⟨ a_i ⟩ ,
        
        b_i = b̃_i +y_i,   y_i=⟨ b_i ⟩ ,
        
        where x_i and y_i are c-numbers describing the condensate, while operators ã_i and b̃_i annihilate excitations over the condensate. We assume α_ij=α_δ with δ=_j-_i and similarly for β_ij. Namely they only depends on δ and not on _i. For our model, it is sufficient to keep δ=±x̂, ±ŷ, ±x̂±ŷ, i.e. the nn and nnn couplings.
        Fourier transform to 𝐤 space, e.g. ã_i→ã_k,  ℋ_MF becomes a quadratic form of operators ã_k, b̃_k and c-numbers x_k, y_k. In accordance with the classical analysis, we assume x_k and y_k are nonzero only at a pair of wave vector ±𝐐/2. It is then diagonalized by a standard Bogoliubov transformation, 
          ℋ_MF   =   ∑_k [c^†_k c_k + d^†_k d_k +1]ω_k +  N/2∑_δ J_δ (|α_δ|^2-|β_δ|^2) 
           +   ∑_𝐤=±𝐐/2 [β_k (|x_k|^2+|y_k|^2) +(ix_ky_-kα^*_k + h.c.)]
           -   Nλ (2S +1) +λ∑_𝐤=±𝐐/2 (|x_k|^2+|y_k|^2).
         
        Here N is the number of lattice sites, c_k and d_k are the eigenmodes of spin excitations with dispersion
        ω_k=√((λ+β_k)^2-|α_k|^2).
        
        and 
        α_k=1/2∑_δ J_δsin (𝐤·δ) α_δ,
        β_k=1/2∑_δ J_δcos (𝐤·δ) β_δ.
        
        We adopt the sprial ansatz x_i=√(2m)cos(𝐐/2·_i), y_i=√(2m)sin(𝐐/2·_i). Then x_Q/2=√(Nm/2), y_Q/2=-ix_Q/2. 
        
        Minimizing the SBMF ground energy with respect to the variational parameters {λ, α_δ, β_δ,x_𝐐/2,y_𝐐/2} leads to the self-consistency equations,
        
        S+1/2=1/2N∑_k λ+β_k/ω_k+m, 
        α_δ=1/2N∑_k α_k/ω_ksin (𝐤·δ)+msin (𝐐/2·δ),
        β_δ=1/2N∑_k λ+β_k/ω_kcos (𝐤·δ)+mcos (𝐐/2·δ),
        λ+β_𝐐/2 = α_𝐐/2.
        
        The last equation is equivalent to the requirement that 𝐐 is chosen to be the minimum of ω_k.
        And the SBMF ground state energy simplifies to
        
        E_MF=∑_k ω_k -Nλ (2S +1) -N/2∑_δ J_δ (|β_δ|^2-|α_δ|^2).
        
        
        
        In the large S limit, we have m≃ S, α_δ= m sin (𝐐/2·δ), β_δ= m cos (𝐐/2·δ), λ=-∑_δ J_δcos ( 𝐐·δ), and 
        
        E_cl=1/2NS^2∑ _δ J_δcos ( 𝐐·δ),
        
        which agrees with the classical result as expected.
        
        
        
        
        § MODIFIED SPIN WAVE THEORY
        
        
        We represent the spin operator using Holstein-Primakoff (HP) bosons,
        
        S^-   =   a^†√(2S-a^† a),
        
        S^+   =   √(2S-a^† a) a,
        
        S_z   =   S-a^† a.
        
        Proper number of boson operators are introduced for the two-sublattice case (Neel phase) and the four-sublattice case (stripe phase). Take the two-sublattice for example, a_i (or b_j) are Bose annihilation operators on the A (or B) sublattice. 
        The dipolar Heisenberg Hamiltonian can then be expanded in series of boson operators, 
        ℋ_SW=ℋ_cl+ℋ^(2)+ℋ^(4),
        
        where the classical part ℋ_cl is given previously in Eq. <ref>, the quadratic part ℋ^(2) is
        ℋ^(2)/S   =   J_x∑_nn_x(a^†_i a_i+b^†_j b_j-a^†_i b^†_j-a_i b_j) 
           +   J_y∑_nn_y(a^†_i a_i+b^†_j b_j-a^†_i b^†_j-a_i b_j) 
           +   J_d'∑_nnn_1(a^†_i a_i'+b^†_j b_j'-a^†_i a_i-b^†_j b_j) 
           +   J_d∑_nnn_2(a^†_i a_i'+b^†_j b_j'-a^†_i a_i-b^†_j b_j),
        
        and the quartic part ℋ^(4) is 
        
        4ℋ^(4)   =   J_x∑_nn_x(a^†_i a_ia_ib_j+a_jb^†_j b_jb_j-2a^†_ia_ib^†_jb_j) 
           +   J_y∑_nn_y(a^†_i a_ia_ib_j+a_jb^†_j b_jb_j-2a^†_ia_ib^†_jb_j) 
           +   J_d'∑_nnn_1(a^†_ia^†_i'a_ia_i'-a^†_i a^†_i'a_i'a_i') 
           +   J_d'∑_nnn_1(b^†_jb^†_j'b_jb_j'-b^†_j b^†_j'b_j'b_j') 
           +   J_d∑_nnn_2(a^†_ia^†_i'a_ia_i'-a^†_i a^†_i'a_i'a_i') 
           +   J_d∑_nnn_2(b^†_jb^†_j'b_jb_j'-b^†_j b^†_j'b_j'b_j')+H.C..
        
        
        For the Neel phase, the expectation values of many operator pairs vanish, e.g., 
        ⟨ a_ia_i⟩=⟨ a_ia_i'⟩=⟨ a^†_ib_j⟩=0.
        
        We define the following nonzero averages of boson operators describing the quantum fluctuations of spins
        
        f_0≡⟨ a^†_ia_i⟩, 
        
        g_1≡⟨ a_ib_j⟩, 
        
        f_2≡⟨ a^†_ia_j⟩,
        
        and apply self-consistent mean field decoupling of the quartic terms in ℋ^(4)  
        a^†_ia_ia_ib_j   =   f_0 a_ib_j+g_1 a^†_ia_i-f_0 g_1, 
        
        a^†_ia^†_i'a_i'a_i'   =   f_0 a^†_ia_i'+f_2 a^†_i'a_i'-f_0 f_2, 
        
        a^†_ia_ib^†_jb_j   =   (1-λ_1)(f_0 b^†_jb_j+f_0 a^†_ia_i-f_0^2) 
           +   λ_1(g_1a_ib_j+g_1a^†_ib^†_j-g_1^2), 
        
        a^†_ia_ia^†_i'a_i'   =   (1-λ_2)(f_0 a^†_i'a_i'+f_0 a^†_ia_i-f_0^2) 
           +   λ_2(f_2a^†_ia_i'+f_2a^†_i'a_i-f_2^2),
        
        where 0≤λ_1,λ_2≤ 1 are parameters determined by minimizing the ground state energy. The magnitudes of λ_1 or λ_2 describe the competition between the diagonal and off-diagonal terms of spin deviation operators. 
        
        After Bogoliubov transformation to diagonalize the resulting Hamiltonian, the self-consistent equations can be solved by minimizing the ground state energy E_0 with respect to the variational parameters {λ_1, λ_2}. 
        The energy and the staggered magnetization are given by
        
        E_0   =   ℋ_cl+E_1+∑_𝐤ϵ_k, 
        ⟨ S_z⟩   =   S-f_0,
        
        where
              E_1=
              J_x[(1-λ_1)(f_0^2-f_0)+λ_1g_1^2+(1/2-f_0)g_1+S] 
           +   J_y[(1-λ_1)(f_0^2-f_0)+λ_1g_1^2+(1/2-f_0)g_1+S] 
           -   J_d'[(1-λ_2)(f_0^2-f_0)+λ_2f_2^2+(1/2-f_0)f_2+S] 
           -   J_d[(1-λ_2)(f_0^2-f_0)+λ_2f_2^2+(1/2-f_0)f_2+S], 
        
        and
        ϵ_k=√(h_k^2-Δ_k^2),
        
        with 
        
        h_k   =   J_x[S-(1-λ_1)f_0+1/2g_1] 
           +   J_y[S-(1-λ_1)f_0+1/2g_1] 
           -   J_d'[S-(1-λ_2)f_0+1/2f_2)] 
           -   J_d[S-(1-λ_2)f_0+1/2f_2)] 
           +   J_d'[cos(k_x+k_y)(S-f_0/2+λ_2f_2)] 
           +   J_d[cos(k_x-k_y)(S-f_0/2+λ_2f_2)],
        
        and
        Δ_k   =   J_xcos(k_x)(f_0/2-S-λ g_1) 
           +   J_ycos(k_y)(f_0/2-S-λ g_1).
        
        
        The self-consistency equations are
        
        f_0   =   1/N∑_𝐤(h_k/ϵ_k-1), 
        
        g_1   =   -1/N∑_𝐤Δ_k/2ϵ_k(cos k_x+cos k_y), 
        
        f_2   =   1/N∑_𝐤1/ϵ_kcos k_xcos k_y.
        
        The criterion for Neel order is a finite ⟨ S_z⟩_Neel.  
        
        For the stripe case, a similar procedure can be applied except that four types of boson operators should be introduced. Correspondingly, three variational parameters λ_i are needed due to the difference between x and y directions. Using a similar self-consistent mean-field approximation, the boundary of stripe phase can be determined. The criterion for stripe phase is a finite ⟨ S_z⟩_stripe and real, positive-definite spin deviation operators.  
        
        The mean field phase diagrams in Fig. 2 of the main text obtained by two different methods give us the same qualitatively picture but  different areas of the spin disordered region. This is not surprising, since different spin representations and mean field decoupling schemes are used. For example, in Eq. <ref>, the expectation values of bond operators are used in SBMF while the quartic terms in Eq. <ref> are described by the variational parameters of λ_1, λ_2 for the modified spin wave theory. 
        
        
        
        
        § TENSOR NETWORK ANSATZ
        
        
         §.§ Simple Update
        
        We choose a L× L unit cell (i.e. L× L local tensors) with different virtual bond dimension D=2,4,6 to form the initial tensor network state |Ψ⟩ and set the time interval τ=0.005J^-1 for imaginary time evolution iterations for local tensors |ψ⟩,
        
        |ψ'⟩=exp(-τ H)|ψ⟩,
        
        until convergence is achieved. 
         
        Taking L=2 as an example and using the Trotter-Suzuki formula <cit.>, we can express the projection operator as
        exp(-τ H)≈∏_i=1^4exp(-τ H_i)+O(τ^2),
        
        where
        
        H_1   =   J_x𝐒_A ·𝐒_B+J_y𝐒_B ·𝐒_C+2J_d𝐒_A ·𝐒_C,
        
        H_2   =   J_x𝐒_C ·𝐒_D+J_y𝐒_A ·𝐒_D+2J_d𝐒_A ·𝐒_C,
        
        H_3   =   J_x𝐒_ C·𝐒_D+J_y𝐒_B ·𝐒_C+2J'_d𝐒_B ·𝐒_D,
        
        H_4   =   J_x𝐒_A ·𝐒_B+J_y𝐒_A ·𝐒_D+2J'_d𝐒_B ·𝐒_D.
        
        This means that each iteration of projection can be done using exp(-τ H_i)(i=1,2,3,4) in four separate steps for one plaquette. While in each step three out of four tensors are evolved (Fig. <ref>). 
        
        
        
         §.§ Tensor Renormalization Group
        
        Starting from the converged local tensors T_i obtained from the simple update, one can construct new two-dimensional local tensors (Fig. <ref>(a)),
        𝒯_i   =   T_i^† T_i,
        𝒪_i   =   T_i^† O_i T_i,
        
        where O_i is an operator.
         
        
        The expectation value of  O_i, ⟨ O_i⟩=⟨Ψ|O_i|Ψ⟩/⟨Ψ|Ψ⟩, can then be obtained by 
        ⟨ O_i⟩= (𝒪_i∏_j≠ i𝒯_j)/∏_j 𝒯_j,
        
        in the thermodynamic limit by using Tensor Renormalization Group (TRG) method <cit.>, where tr stands for contraction of neighboring tensors.
        Taking the denominator ∏_i 𝒯_i as an example (the numerator can be coarse-grained with the same procedure because the local operator O has the same structure with 𝒯). As shown in Fig. <ref>(b), for each step, one can decompose each 𝒯 to two 𝒮 via singular value decomposition, 
        𝒯_i_1,a,k_1,b≈∑_ĩ𝒮_i_1,k_1,ĩ𝒮_a,b,ĩ.
        
        The truncation bond dimension of coarse-graining (dimension of the third leg of 𝒮) is set as χ. The new local tensor 𝒯̃ with the same structure as 𝒯 can be constructed from contracting the inner legs of four 𝒮,
        𝒯̃_ĩ,j̃,k̃,l̃=∑_a,b,c,d𝒮_a,b,ĩ𝒮_d,c,j̃𝒮_a,c,k̃𝒮_d,b,l̃.
        
        Using 𝒯̃ as the starting tensors, these steps are repeated until ∏_i 𝒯_i is converged. In our TRG calculation, the truncation bond dimension is fixed as χ = 8 to make D the only tuning parameter of the whole procedure.
        
        
        
        
         §.§ Comparison of Different Unit Cell Sizes
         
        The phase boundary between the quantum paramagnetic phase and other long-range orders can be inferred from the disappearance of magnetic order parameters. It is crucial to determine whether the phase boundary depends sensitively on L, the size of the unit cell. To address this question, we calculate the magnetization M and the average energy E with different unit cell size (L=2,4,8) at fixed D. Fig. <ref> shows examples with D=4, from which we conclude that increasing L does not increase the accuracy significantly. Thus the scaling of L to larger value gives essentially the same result as L=2 and we can use L=2 to obtain the phase diagram in the main text. 
        
        
        
         §.§ Average Energy and Extrapolation of M
         
        The simple update and coarse-graining TRG steps are repeated until the average energy E (Fig. <ref>) is converged for given D. To obtain the phase boundary, we apply the finite-size extrapolations of M using second-order polynomial fit in 1/D to infinite D <cit.>. One example at ϕ=15^∘ is shown in Fig. <ref>. Suppression of the magnetization to zero as D→∞ suggests a quantum paramagnetic region. 
        
        
         §.§ Results for Finite Anisotropy
        
        
        We apply the same extrapolations of M for different anisotropy η (Fig. <ref>), which shows that the quantum paramagnetic region persists away from the Heisenberg limit η=1. Specifically, for η< 1, the quantum paramagnetic region remains robust for a large region, e.g., down to η=0.5. While for η>1, long range order is preferred when η is increased to η = 1.1. This seems to suggest that the Heisenberg limit is close to the upper limit of the quantum paramagnetic region.
        7
        fxundefined [1]
         ifx#1
        fnum [1]
          #1firstoftwo
         secondoftwo
         
        fx [1]
          #1firstoftwo
         secondoftwo
         
        
        
        
        
        
        noop [0]secondoftwo
        
        ref[1]@startlink#1@href
        href[1]#1@endlink
        anitize@url [0] `
        12 `$12
          `&12 `#12 `1̂2 `_12 `%12
        startlink[1]
        endlink[0]
        
        rl [1]href #1
        
        
        
        
        
        
        
        
        
        
        
        
        @bib@innerbibempty
        
         [Luttinger and Tisza(1946)]LuttingerTisza2
          authorauthorJ. M. Luttinger and authorL. Tisza,  10.1103/PhysRev.70.954journaljournalPhys. Rev. volume70,pages954 (year1946)NoStop
         [Jiang et al.(2008)Jiang,
          Weng, and Xiang]simpleup2
          authorauthorH. C. Jiang, authorZ. Y. Weng,   and authorT. Xiang,  10.1103/PhysRevLett.101.090603journaljournalPhys. Rev. Lett. volume101, pages090603 (year2008)NoStop
         [Corboz et al.(2010)Corboz,
          Jordan, and Vidal]nnnsimple2
          authorauthorP. Corboz, authorJ. Jordan,   and authorG. Vidal,  10.1103/PhysRevB.82.245119journaljournalPhys. Rev. B volume82, pages245119 (year2010)NoStop
         [Levin and Nave(2007)]LevinTRG2
          authorauthorM. Levin and authorC. P. Nave,  10.1103/PhysRevLett.99.120601journaljournalPhys. Rev. Lett. volume99, pages120601 (year2007)NoStop
         [Xie et al.(2012)Xie,
          Chen, Qin, Zhu, Yang, and Xiang]HOTRG12
          authorauthorZ. Y. Xie, authorJ. Chen, authorM. P. Qin, authorJ. W. Zhu, authorL. P. Yang,  and authorT. Xiang,  10.1103/PhysRevB.86.045139journaljournalPhys. Rev. B volume86, pages045139 (year2012)NoStop
         [Jiang et al.(2012)Jiang,
          Yao, and Balents]J1J2Jiang2
          authorauthorH.-C. Jiang, authorH. Yao,  andauthorL. Balents,  10.1103/PhysRevB.86.024424journaljournalPhys. Rev. B volume86, pages024424 (year2012)NoStop
         [Wang et al.(2016)Wang,
          Gu, Verstraete, and Wen]tnj1j22
          authorauthorL. Wang, authorZ.-C. Gu,
          authorF. Verstraete,  andauthorX.-G. Wen,  10.1103/PhysRevB.94.075143journaljournalPhys. Rev. B volume94, pages075143 (year2016)NoStop
