Frustrated Magnetism of Dipolar Molecules on a Square Optical Lattice: Prediction of a Quantum Paramagnetic Ground State

TL;DR

This paper predicts that dipolar molecules on a square optical lattice can host a quantum paramagnetic ground state, offering a new route to realize quantum spin liquids without complex lattice geometries.

Contribution

It demonstrates that tilting dipoles in a dipolar Heisenberg model induces frustration leading to a quantum spin liquid phase on a square lattice, which is a novel approach.

Findings

Identification of a quantum paramagnetic region between ordered phases.

Consistent phase diagrams from multiple theoretical methods.

Confirmation of a quantum spin liquid state using tensor network techniques.

Abstract

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$\acute{\textrm{e}}$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.

Equations

$$H_{XXZ}=\frac{J}{2}\sum _{i\ne j}f({\mathbf{r}}_i-{\mathbf{r}}_j)(S_i^x{S}_j^x+S_i^y{S}_j^y+\eta S_i^z{S}_j^z).$$

$$f(\mathbf{r})=[1-3(\hat{r}\cdot \hat{d}{)}^2](a/r{)}^3,$$

$$H_d=\frac{J}{2}\sum _{i\ne j}f({\mathbf{r}}_i-{\mathbf{r}}_j){\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$${\mathbf{S}}_{\mathbf{r}}/S=\hat{x}\cos (\mathbf{Q}\cdot \mathbf{r})+\hat{y}\sin (\mathbf{Q}\cdot \mathbf{r}),$$

$$2{\mathcal{H}}_{\text{cl}}/N{S}^2$$

$$S^{+}=a^{†}b,$$

$$S^{-}=b^{†}a,$$

$$S_z=\frac{1}{2}(a^{†}a-b^{†}b).$$

$$\frac{1}{2}(a^{†}a+b^{†}b)=S.$$

$$A_{ij}=\frac{1}{2}(a_i{b}_j-b_i{a}_j),$$

$$B_{ij}=\frac{1}{2}(a_i^{†}{a}_j+b_i^{†}{b}_j).$$

$${\mathbf{S}}_i\cdot {\mathbf{S}}_j=:B_{ij}^{†}{B}_{ij}:-A_{ij}^{†}{A}_{ij},$$

$${\mathbf{S}}_i\cdot {\mathbf{S}}_j\simeq [{\beta }_{ij}{B}_{ij}-{\alpha }_{ij}{A}_{ij}+h.c.]-|{\beta }_{ij}{|}^2+|{\alpha }_{ij}{|}^2,$$

$${\alpha }_{ij}=⟨A_{ij}^{†}⟩,$$

$${\beta }_{ij}=⟨B_{ij}^{†}⟩,$$

$${\mathcal{H}}_{MF}$$

$$a_i={\tilde{a}}_i+x_i,x_i=⟨a_i⟩,$$

$$b_i={\tilde{b}}_i+y_i,y_i=⟨b_i⟩,$$

$${\mathcal{H}}_{MF}$$

$${\omega }_k=\sqrt{(\lambda +{\beta }_k{)}^2-|{\alpha }_k{|}^2}.$$

$${\alpha }_k=\frac{1}{2}\sum _{\delta }{J}_{\delta }\sin (\mathbf{k}\cdot \mathbf{\delta }){\alpha }_{\delta },$$

$${\beta }_k=\frac{1}{2}\sum _{\delta }{J}_{\delta }\cos (\mathbf{k}\cdot \mathbf{\delta }){\beta }_{\delta }.$$

$$S+\frac{1}{2}=\frac{1}{2N}\sum _k\frac{\lambda +{\beta }_k}{{\omega }_k}+m,$$

$${\alpha }_{\delta }=\frac{1}{2N}\sum _k\frac{{\alpha }_k}{{\omega }_k}\sin (\mathbf{k}\cdot \mathbf{\delta })+m\sin (\frac{\mathbf{Q}}{2}\cdot \mathbf{\delta }),$$

$${\beta }_{\delta }=\frac{1}{2N}\sum _k\frac{\lambda +{\beta }_k}{{\omega }_k}\cos (\mathbf{k}\cdot \mathbf{\delta })+m\cos (\frac{\mathbf{Q}}{2}\cdot \mathbf{\delta }),$$

$$\lambda +{\beta }_{\frac{\mathbf{Q}}{2}}={\alpha }_{\frac{\mathbf{Q}}{2}}.$$

$$E_{MF}=\sum _k{\omega }_k-N\lambda (2S+1)-\frac{N}{2}\sum _{\delta }{J}_{\delta }(|{\beta }_{\delta }{|}^2-|{\alpha }_{\delta }{|}^2).$$

$$E_{\text{cl}}=\frac{1}{2}N{S}^2\sum _{\delta }{J}_{\delta }\cos (\mathbf{Q}\cdot \mathbf{\delta }),$$

$$S^{-}$$

$$S^{+}$$

$$S_z$$

$${\mathcal{H}}_{SW}={\mathcal{H}}_{\text{cl}}+{\mathcal{H}}^{(2)}+{\mathcal{H}}^{(4)},$$

$${\mathcal{H}}^{(2)}/S$$

$$4{\mathcal{H}}^{(4)}$$

$$⟨a_i{a}_i⟩=⟨a_i{a}_{i^{\prime }}⟩=⟨a_i^{†}{b}_j⟩=0.$$

$$f_0\equiv ⟨a_i^{†}{a}_i⟩,$$

Full text

Understanding highly entangled quantum matter remains a challenging goal of condensed matter physics Savary and Balents (2017). 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 Savary and Balents (2017); Zhou et al. (2017); Diep (2013). 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 Wen (2002). While the existence of spin liquids has been firmly established in a number of exactly solvable models, e.g., the toric code Kitaev (2003) or the honeycomb Kitaev model Kitaev (2006), 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 Yan et al. (2011); Iqbal et al. (2013); Jiang et al. (2012); Gong et al. (2014); Wang et al. (2016). An unambiguous experimental identification of quantum spin liquids in solid state materials also seems elusive Savary and Balents (2017). 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 de Paz et al. (2013) and polar molecules  Yan et al. (2013); Hazzard et al. (2014) 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 Gorshkov et al. (2011a); Yao et al. (2015). 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 Verstraete and Cirac (2004); Jordan et al. (2008); Corboz et al. (2010a); Levin and Nave (2007); Xie et al. (2012). 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,

[TABLE]

Here $i$ and $j$ label the lattice sites, $\mathbf{S}_{i}=(S^{x}_{i},S^{y}_{i},S^{z}_{i})$ are the spin (or pseudospin) operators at site $i$, and $\eta$ is the exchange anisotropy. The key new feature here is that the coupling between the two spins depends on their relative position $\mathbf{r}=\mathbf{r}_{i}-\mathbf{r}_{j}$ and the external field (dipole) direction $\hat{d}$

[TABLE]

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 $\eta=1$, $H_{XXZ}$ reduces to the dipolar Heisenberg model

[TABLE]

and for $\eta=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. de Paz et al. (2013), the spin dynamics of a gas of 52Cr atoms in optical lattices was observed. Each Cr atom carries a magnetic moment of $7\mu_{B}$ and hyperfine spin $S=3$. An external magnetic field is used to align the magnetic dipoles in the direction of $\hat{d}$. Such a dipolar gas of Cr in a deep lattice is shown to be described by $H_{XXZ}$ with $J=-\mu_{0}(g\mu_{B})^{2}/4\pi a^{3}<0$ and $\eta=-2$ de Paz et al. (2013). 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 40K87Rb confined in optical lattices with negligible tunneling provide another way to realize $H_{XXZ}$ with $S=1/2$ and tunable $J$ and $\eta$ Yan et al. (2013). Each molecule carries an electric dipole moment $\mathbf{d}$ and undergoes rotation with angular momentum $\mathbf{J}$ [see Fig. 1(a)]. Here, the pseudospin 1/2 refers to two rotational states of the molecule labeled by $|j,m\rangle$, where $j$ is the quantum number of the rotational angular momentum $\mathbf{J}$ 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. Yan et al. (2013); Gorshkov et al. (2011b); Yao et al. (2015). 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\rangle=|0,0\rangle$ and $|1,0\rangle$ as the pseudospin down and up respectively, Refs. Gorshkov et al. (2011b); Yao et al. (2015) showed that the system is described by the effective Hamiltonian $H_{XXZ}$ with $J=D_{t}^{2}/2\pi\epsilon_{0}a^{3}>0$ and $\eta=(D_{1}-D_{0})^{2}/2D_{t}^{2}>0$. Here the dipole matrix element $D_{t}=\langle 1,0|d^{0}|0,0\rangle$, $D_{1}=\langle 1,0|d^{0}|1,0\rangle$, $D_{0}=\langle 0,0|d^{0}|0,0\rangle$, and $d^{0}$ together with $d^{\pm}$ form the vector dipole operator in the spherical basis Gorshkov et al. (2011b); Yao et al. (2015).

The anisotropy $\eta$ increases monotonically with $E$. As shown in Ref. Yao et al. (2015), when $E\simeq 1.7B/|\mathbf{d}|$ with $B$ the energy splitting of the two pseudospin states, $\eta=1$, and one arrives at the dipolar Heisenberg model $H_{d}$. In the KRb experiment Yan et al. (2013) carried out at zero field and cubic lattice, $\eta\rightarrow 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 Yan et al. (2013) and modeled theoretically in Ref. Hazzard et al. (2014). Recently Yao et al. Yao et al. (2015) considered general $\eta$ 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, $\eta=1$ and $\theta=0$, in which $\theta$ is defined by $\hat{d}\cdot\hat{x}=\sin\theta\cos\phi$ with $\hat{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 $\eta$.

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 $\theta$, 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 $\theta$ values (not around $\theta=0$ as in Ref. Yao et al. (2015)) between the N$\acute{\textrm{e}}$el and the stripe orders. Thus, it differs qualitatively from the spin liquids studied in Ref. Yao et al. (2015). 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 $\hat{d}$ is tuned as well as its connection to frustrated quantum spin models Balents (2010); Diep (2013), let us consider the leading exchange couplings between the nearest neighbors, $J_{x}=Jf(a\hat{x})$ and $J_{y}=Jf(a\hat{y})$, and the next nearest neighbors, $J_{d}=Jf(a\hat{x}+a\hat{y})$ and $J^{\prime}_{d}=Jf(a\hat{x}-a\hat{y})$ [Fig. 1(b)]. Their relative magnitudes and signs depend sensitively on the dipole tilting angle $\theta$ and $\phi$. One example is shown in Fig. 1(b) for fixed $\phi=25^{\circ}$. At small $\theta$, $J_{x}\sim J_{y}$ dominates because it is about three times that of $J_{d}\sim J^{\prime}_{d}$. The situation is reminiscent of the $J_{1}$-$J_{2}$ model in the regime of the N$\acute{\textrm{e}}$el order. As $\theta$ is increased, $J_{d}$ and $J^{\prime}_{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 $\theta\simeq 40^{\circ}$, $J_{x}$ and $J_{d}$ vanish while $J^{\prime}_{d}\sim 0.4J_{y}$. The model can be viewed as coupled Heisenberg chains. For even larger $\theta$, $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, $\theta=\phi=0$, where $J_{2}/J_{1}=1/2\sqrt{2}\approx 0.35$ and the system is N$\acute{\textrm{e}}$el ordered.

The degree of frustration can be measured by the “spin gap” $\Delta$, the energy difference between the ground and the first excited state, from exact diagonalization of $H_{d}$ for a $4\times 4$ lattice not . For example, we observe a pronounced peak in $\Delta$ around $\theta\sim 28^{\circ}$ for $\phi=25^{\circ}$, which indicates strong frustration and points to a gapped, spin disordered ground state Sindzingre (2004). For fixed $\phi=35^{\circ}$, the spin structure factor shows a clear peak at $(\pi,\pi)$ for $\theta\sim 15^{\circ}$ for the N$\acute{\textrm{e}}$el order, a peak at $(0,\pi)$ for $\theta\sim 50^{\circ}$ for the stripe order, but no well defined peaks around $\theta\sim 35^{\circ}$, consistent with the argument above.

Spin-wave and Schwinger-boson theory.—First, we obtain a coarse phase diagram of $H_{d}$ on the $(\theta,\phi)$ 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 Luttinger and Tisza (1946). $H_{d}$ is of the form $\sum_{ij}J_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}$ with hard spin constraint $\mathbf{S}_{i}=S$ and $J_{ij}$ only depends on $\mathbf{r}_{i}-\mathbf{r}_{j}$. A theorem states that the classical ground state is a planar spin spiral, $\mathbf{S}_{\mathbf{r}}/S=\hat{x}\cos(\mathbf{Q}\cdot\mathbf{r})+\hat{y}\sin(\mathbf{Q}\cdot\mathbf{r})$ with an ordering wave vector $\mathbf{Q}=(Q_{x},Q_{y})$ Kaplan and Menyuk (2007). The classical phase diagram not consists of three phases. The first is the N$\acute{\textrm{e}}$el order corresponding to $\mathbf{Q}=(\pi,\pi)$ for small $\theta$. The second is the stripe phase with $\mathbf{Q}=(0,\pi)$ for large $\theta$ but not too large $\phi$. These two spin orders are collinear. The third, spiral phase fills the rest of the phase diagram, for large $\theta$ and $\phi$, where $\mathbf{Q}$ 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  Takahashi (1989); Chu and Shen (1991); Cui (1992). 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 not . The result is summarized in Fig. 2(a). We find that the phase boundary of the N$\acute{\textrm{e}}$el (stripe) phase moves towards smaller (larger) $\theta$ values, opening up an intermediate region in between where the magnetization vanishes. The spiral phase also recedes to higher $\phi$ 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$ Arovas and Auerbach (1988); Read and Sachdev (1991). It is a well tested method capable of describing both magnetically ordered and spin liquid states of frustrated spin models Mezio et al. (2013); Merino et al. (2014); Lee et al. (2014); Yang and Wang (2016). The resulting phase diagram is shown in Fig. 2(b). Here, each magnetic order corresponds to condensation of bosons at a certain wave vector $\mathbf{Q}$. Within a finite strip region labeled by QP between the N$\acute{\textrm{e}}$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 Verstraete and Cirac (2004); Jordan et al. (2008); Corboz et al. (2010a) has recently emerged as an accurate and unbiased algorithm for solving two dimensional frustrated quantum spin models Wang et al. (2016); Zou et al. (2016); Jiang et al. (2016); Liao et al. (2017). In this approach, the ground state many-body wave function $|\Psi\rangle$ is constructed from a network of tensors $T_{i}$ defined on lattice site $i$: $|\Psi\rangle=\mbox{tr}\prod_{i}T_{i}$, where tr 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\times 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, $|\psi^{\prime}\rangle=\exp(-\tau H)|\psi\rangle$, until convergence is achieved. We adopt the simple update scheme Jiang et al. (2008) based on singular value decomposition. By using the Trotter-Suzuki formula $\exp(-\tau H)\approx\prod_{i=1}^{4}\exp(-\tau H_{i})+O(\tau^{2})$, each iteration of projection for one plaquette can be done using $\exp(-\tau 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^{\prime}_{d}$ terms (See Refs. Corboz et al. (2010b); Corboz and Mila (2013); Wang et al. (2016); not ).

The expectation value of a local operator $O_{j}$ at site $j$, $\langle O_{j}\rangle=\langle\Psi|O_{j}|\Psi\rangle/\langle\Psi|\Psi\rangle$, can be computed by tensor contraction, $\langle O_{j}\rangle=\mbox{tr}(\mathcal{O}_{j}\prod_{i\neq j}\mathcal{T}_{i})/\mbox{tr}\prod_{i}\mathcal{T}_{i}$ where $\mathcal{T}_{i}=T_{i}^{\dagger}T_{i}$ and $\mathcal{O}_{j}=T_{j}^{\dagger}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 Levin and Nave (2007); Xie et al. (2012). In this way, we calculate the order parameters such as magnetization $M=\sqrt{\langle S_{x}\rangle^{2}+\langle S_{y}\rangle^{2}+\langle S_{z}\rangle^{2}}$ not .

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 $\phi=15^{\circ}$. By extrapolating the results to infinite $D$, we can determine the phase boundary of the N$\acute{\textrm{e}}$el and stripe phases. Repeating the procedure for different $\phi$ 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 $\theta$ and persisting from $\phi=0$ up to $\phi=20^{\circ}$, 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$\acute{\textrm{e}}$el phase, there is a sudden drop of $M$. Note that the spiral phase, in general, is incompatible with the $L\times L$ cluster choice, even for large $L$. So we refrain from carrying out the tensor network ansatz beyond $\phi=20^{\circ}$. 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 not . Finally, we point out that the quantum paramagnetic phase is a robust feature of the dipolar $XXZ$ model. It persists when $\eta$ is tuned away from the Heisenberg limit, e.g., down to $\eta=0.5$ not .

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 Gong et al. (2014) 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 Yang and Wang (2016). 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 C4 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 Jiang and Ran (2015); Mei et al. (2017) and Lanczos iteration Huang et al. (2016) 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 $(\theta,\phi)$ 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 Parsons et al. (2015); Cheuk et al. (2015); Parsons et al. (2016); Boll et al. (2016); Cheuk et al. (2016). Thus, it seems possible to probe the spin order or spin correlations of $H_{d}$ and related models in future experiments.

I Classical Phase Diagram

The classical ground state of a translationally invariant spin model $\sum_{ij}J_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}$ on a Bravais lattice can be obtained by minimizing the energy within the planar helix ansatz

[TABLE]

where $\hat{x}$ and $\hat{y}$ form an orthonormal basis and $\mathbf{Q}=(Q_{x},Q_{y})$ is the ordering wavevector Luttinger and Tisza (1946). This variational ansatz satisfies the hard spin constraint $|\mathbf{S}_{\mathbf{r}}|=S$. The classical energy of the dipolar Heisenberg model depends on $\mathbf{Q}$ via

[TABLE]

The set of wavevectors $\mathbf{Q}$ minimizing the classical energy will be denoted as $\{\mathbf{Q}\}$. We compare the energies of the incommensurate spiral $\{\mathbf{Q}_{I}\}$, the stripe $\mathbf{Q}_{s}=(0,\pi)$ and the N$\acute{\textrm{e}}$el $\mathbf{Q}_{n}=(\pi,\pi)$ order. The result is the classical phase diagram shown in Fig. S1. 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 $\phi=45^{\circ}$ (see the inset of Fig. S1), $Q_{x}=Q_{y}\equiv Q$, where $Q$ changes continuously from $\pi$ on the phase boundary between Neel and spiral phase to $\pi/4$ at the upper right corner of the $(\theta,\phi)$ diagram.

II Exact Diagonalization

We calculate the “spin gap” $\Delta$, the energy defference between the ground and the first excited state from exact diagonalization of the dipolar Heisenberg Hamiltonian for a $4\times 4$ lattice. Fig. S2 shows $\Delta$ as functions of $\theta$ for different $\phi$. For $\phi\leq 30^{\circ}$, the pronounced peak at each line indicates strong frustration. For $\phi=35^{\circ}$, The appearance of the second peak corresponds to the transition from the stripe phase to the incommensurate spiral phase as $\theta$ increases. For $\phi\geq 40^{\circ}$, the disappearance of the first peak indicates the transition from the Neel phase to the spiral phase.

III 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

[TABLE]

with the constraint

[TABLE]

For the square lattice, introduce the antiferromagnetic ($A$) and ferromagnetic ($B$) bond operators

[TABLE]

In terms of the bond operators, the spin exchange term becomes

[TABLE]

where $::$ means normal order of bosonic operators. Note that $B$ and $A$ are related by operator identity $:B^{\dagger}_{ij}B_{ij}:+A^{\dagger}_{ij}A_{ij}=S^{2}$. We adopt the rotational invariant formulation of SBMF and perform mean field decoupling for both $A$ and $B$,

[TABLE]

where

[TABLE]

and $\langle...\rangle$ 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. (S4) is only enforced on average by introducing the Lagrange multiplier $\lambda$. The mean field Hamiltonian then takes the form

[TABLE]

Long range magnetic order corresponds to condensation of the $a$ and/or $b$ bosons. To treat the condensate fraction, we decompose each operator into

[TABLE]

where $x_{i}$ and $y_{i}$ are c-numbers describing the condensate, while operators $\tilde{a}_{i}$ and $\tilde{b}_{i}$ annihilate excitations over the condensate. We assume $\alpha_{ij}=\alpha_{\delta}$ with $\boldsymbol{\delta}=\mathbf{r}_{j}-\mathbf{r}_{i}$ and similarly for $\beta_{ij}$. Namely they only depends on $\boldsymbol{\delta}$ and not on $\mathbf{r}_{i}$. For our model, it is sufficient to keep $\boldsymbol{\delta}=\pm\hat{x}$, $\pm\hat{y}$, $\pm\hat{x}\pm\hat{y}$, i.e. the nn and nnn couplings. Fourier transform to $\mathbf{k}$ space, e.g. $\tilde{a}_{i}\rightarrow\tilde{a}_{k}$, $\mathscr{H}_{MF}$ becomes a quadratic form of operators $\tilde{a}_{k}$, $\tilde{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 $\pm\mathbf{Q}/2$. It is then diagonalized by a standard Bogoliubov transformation,

[TABLE]

Here $N$ is the number of lattice sites, $c_{k}$ and $d_{k}$ are the eigenmodes of spin excitations with dispersion

[TABLE]

and

[TABLE]

We adopt the sprial ansatz $x_{i}=\sqrt{2m}\cos(\frac{\mathbf{Q}}{2}\cdot\mathbf{r}_{i})$, $y_{i}=\sqrt{2m}\sin(\frac{\mathbf{Q}}{2}\cdot\mathbf{r}_{i})$. Then $x_{Q/2}=\sqrt{Nm/2}$, $y_{Q/2}=-ix_{Q/2}$.

Minimizing the SBMF ground energy with respect to the variational parameters $\{\lambda,\alpha_{\delta},\beta_{\delta},x_{\frac{\mathbf{Q}}{2}},y_{\frac{\mathbf{Q}}{2}}\}$ leads to the self-consistency equations,

[TABLE]

The last equation is equivalent to the requirement that $\mathbf{Q}$ is chosen to be the minimum of $\omega_{k}$. And the SBMF ground state energy simplifies to

[TABLE]

In the large $S$ limit, we have $m\simeq S$, $\alpha_{\delta}=m\sin(\frac{\mathbf{Q}}{2}\cdot\boldsymbol{\delta})$, $\beta_{\delta}=m\cos(\frac{\mathbf{Q}}{2}\cdot\boldsymbol{\delta})$, $\lambda=-\sum_{\delta}J_{\delta}\cos(\mathbf{Q}\cdot\boldsymbol{\delta})$, and

[TABLE]

which agrees with the classical result as expected.

IV Modified Spin Wave Theory

We represent the spin operator using Holstein-Primakoff (HP) bosons,

[TABLE]

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,

[TABLE]

where the classical part $\mathscr{H}_{\textrm{cl}}$ is given previously in Eq. S2, the quadratic part $\mathscr{H}^{(2)}$ is

[TABLE]

and the quartic part $\mathscr{H}^{(4)}$ is

[TABLE]

For the Neel phase, the expectation values of many operator pairs vanish, e.g.,

[TABLE]

We define the following nonzero averages of boson operators describing the quantum fluctuations of spins

[TABLE]

and apply self-consistent mean field decoupling of the quartic terms in $\mathscr{H}^{(4)}$

[TABLE]

where $0\leq\lambda_{1},\lambda_{2}\leq 1$ are parameters determined by minimizing the ground state energy. The magnitudes of $\lambda_{1}$ or $\lambda_{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 $\{\lambda_{1},\lambda_{2}\}$. The energy and the staggered magnetization are given by

[TABLE]

where

[TABLE]

and

[TABLE]

with

[TABLE]

and

[TABLE]

The self-consistency equations are

[TABLE]

The criterion for Neel order is a finite $\langle S_{z}\rangle_{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 $\mathbf{\lambda}_{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 $\langle S_{z}\rangle_{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. S8, the expectation values of bond operators are used in SBMF while the quartic terms in Eq. IV are described by the variational parameters of ${\lambda_{1},\lambda_{2}}$ for the modified spin wave theory.

V Tensor Network Ansatz

V.1 Simple Update

We choose a $L\times L$ unit cell (i.e. $L\times L$ local tensors) with different virtual bond dimension $D=2,4,6$ to form the initial tensor network state $|\Psi\rangle$ and set the time interval $\tau=0.005J^{-1}$ for imaginary time evolution iterations for local tensors $|\psi\rangle$,

[TABLE]

until convergence is achieved.

Taking $L=2$ as an example and using the Trotter-Suzuki formula Jiang et al. (2008); Corboz et al. (2010), we can express the projection operator as

[TABLE]

where

[TABLE]

This means that each iteration of projection can be done using $\exp(-\tau 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. S3).

V.2 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. S4(a)),

[TABLE]

where $O_{i}$ is an operator.

The expectation value of $O_{i}$, $\langle O_{i}\rangle=\langle\Psi|O_{i}|\Psi\rangle/\langle\Psi|\Psi\rangle$, can then be obtained by

[TABLE]

in the thermodynamic limit by using Tensor Renormalization Group (TRG) method Levin and Nave (2007); Xie et al. (2012), where tr stands for contraction of neighboring tensors. Taking the denominator $\mbox{tr}\prod_{i}\mathcal{T}_{i}$ as an example (the numerator can be coarse-grained with the same procedure because the local operator $O$ has the same structure with $\mathcal{T}$). As shown in Fig. S4(b), for each step, one can decompose each $\mathcal{T}$ to two $\mathcal{S}$ via singular value decomposition,

[TABLE]

The truncation bond dimension of coarse-graining (dimension of the third leg of $\mathcal{S}$) is set as $\chi$. The new local tensor $\tilde{\mathcal{T}}$ with the same structure as $\mathcal{T}$ can be constructed from contracting the inner legs of four $\mathcal{S}$,

[TABLE]

Using $\tilde{\mathcal{T}}$ as the starting tensors, these steps are repeated until $\mbox{tr}\prod_{i}\mathcal{T}_{i}$ is converged. In our TRG calculation, the truncation bond dimension is fixed as $\chi=8$ to make $D$ the only tuning parameter of the whole procedure.

V.3 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. S5 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.

V.4 Average Energy and Extrapolation of $M$

The simple update and coarse-graining TRG steps are repeated until the average energy $E$ (Fig. S6) 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$ Jiang et al. (2012); Wang et al. (2016). One example at $\phi=15^{\circ}$ is shown in Fig. S7. Suppression of the magnetization to zero as $D\rightarrow\infty$ suggests a quantum paramagnetic region.

V.5 Results for Finite Anisotropy

We apply the same extrapolations of $M$ for different anisotropy $\eta$ (Fig. S8), which shows that the quantum paramagnetic region persists away from the Heisenberg limit $\eta=1$. Specifically, for $\eta<1$, the quantum paramagnetic region remains robust for a large region, e.g., down to $\eta=0.5$. While for $\eta>1$, long range order is preferred when $\eta$ is increased to $\eta=1.1$. This seems to suggest that the Heisenberg limit is close to the upper limit of the quantum paramagnetic region.