Quantum dimer models emerging from large-spin ultracold atoms

TL;DR

This paper proposes a method using cold atoms in optical lattices to simulate quantum dimer models, enabling the exploration of exotic phases like spin liquids and resonating bonds that are difficult to study in solid-state systems.

Contribution

It introduces an experimental protocol with large-spin ultracold atoms and derives an effective Hamiltonian similar to quantum dimer models, facilitating the study of complex quantum phases.

Findings

Effective Hamiltonian similar to quantum dimer models

Numerical analysis of ground states and phase tuning

Protocol for direct measurement of dimer resonances

Abstract

We propose an experimental protocol for using cold atoms to create and probe quantum dimer models, thereby exploring the Pauling-Anderson vision of a macroscopic collection of resonating bonds. This process can allow the study of exotic crystalline phases, fractionalization, topological spin liquids, and the relationship between resonating dimers and superconductivity subjects which have been challenging to address in solid-state experiments. Our key technical development is considering the action of an off-resonant photoassociation laser on large spin atoms localized at the sites of a deep optical lattice. The resulting superexchange interaction favors nearest-neighbor singlets. We derive an effective Hamiltonian in terms of these dimer degrees of freedom, finding that it is similar to well-known quantum dimer models, which boast a rich variety of valence bond crystal and spin liquid phases. We numerically study the ground state, explain how to tune the parameters, and develop a protocol to directly measure the dimers and their resonating patterns.

Tables

Table 1. Table 1: List of ring exchange amplitudes and bond interactions obtained from Eq. ( 4 ), for different lattice geometries. a There is some debate in the literature about the phases of dimer models on a square lattice Banerjee  et al. ( 2014 ) .

Lattice geometry	$\frac{t}{J^2/U_0}$	$\frac{V}{J^2/U_0}$	$\frac{t^{\prime }}{J^2/U_0}$	Ground state
				at large $f$
Square lattice	$\frac{4}{2f+1}$	$\frac{8}{{\left(2f+1\right)}^2}$	$\frac{8}{{\left(2f+1\right)}^2}$	columnara
Triangular lattice	$\frac{8\left(f+1\right)}{{\left(2f+1\right)}^2}$	$\frac{4}{{\left(2f+1\right)}^2}$	$\frac{8}{{\left(2f+1\right)}^2}$	$\sqrt{12}\times \sqrt{12}$
Cubic lattice	$\frac{4}{2f+1}$	$\frac{8}{{\left(2f+1\right)}^2}$	$\frac{8}{{\left(2f+1\right)}^2}$
Honeycomb lattice	$\frac{12}{{\left(2f+1\right)}^2}$	$𝒪\left(1/f^4\right)$	$𝒪\left(1/f^4\right)$	plaquette
Kagome lattice	$\frac{12}{{\left(2f+1\right)}^2}$	$𝒪\left(1/f^4\right)$	$𝒪\left(1/f^4\right)$

Equations

$$\hat{H}=-J\sum _{⟨ij⟩}\sum _{m=-f}^f{\hat{b}}_{i,m}^{†}{\hat{b}}_{j,m}^{}+\sum _F{U}_F\sum _{M=-F}^F{\hat{A}}_{ii}^{F,M†}{\hat{A}}_{ii}^{F,M},$$

$$\begin{array}[]{rcr}\hat{A}_{ij}^{FM\dagger}&=\frac{1}{\sqrt{2}}\sum_{m}C^{F,M}_{m,M-m}\hat{b}_{i,m}^{\dagger}\hat{b}_{i,M-m}^{\dagger},&{\rm when}\ i=j,\\ &=\sum_{m}C^{F,M}_{m,M-m}\hat{b}_{i,m}^{\dagger}\hat{b}_{j,M-m}^{\dagger},&{\rm when}\ i\neq j.\end{array}$$

$${\hat{H}}_{\mathrm{eff}}$$

$${\hat{H}}_{\mathrm{eff}}\approx -\frac{2J^2}{U_0}\sum _{⟨ij⟩}{\hat{A}}_{ij}^{00†}{\hat{A}}_{ij}^{00}.$$

$$U_F^{\mathrm{Fesh}}={\alpha }_F\frac{{\Omega }^2}{\delta +i\Gamma /2},$$

$${\alpha }_F=\frac{(2i+1)(2f+1)-F(F+1)}{2(2i+1{)}^2}.$$

$$|a⟩=\prod _{(i,j)\in a}{\hat{A}}_{ij}^{†}|0⟩.$$

$${\hat{A}}_{ij}^{†}{\hat{A}}_{ij}^{}|a⟩$$

$$\hat{H}|a⟩=\sum _b|b⟩H_{ba}.$$

$${\hat{H}}_{\mathrm{QDM}}=$$

$$-t^{\prime }\left(|\widetilde{=\mid }⟩⟨\widetilde{\mid =}|+|\widetilde{\mid =}⟩⟨\widetilde{=\mid }|\right)+\cdots$$

$${\hat{H}}_{\mathrm{eff}}=-\frac{2J^2}{U_0}\sum _{⟨ij⟩}{\hat{A}}_{ij}^{00†}{\hat{A}}_{ij}^{00}-\frac{2J^2}{U_2}\sum _{⟨ij⟩M}{\hat{A}}_{ij}^{2M†}{\hat{A}}_{ij}^{2M},$$

$${\hat{H}}_{\mathrm{Fesh}}=$$

$$+\Omega \left(e^{-i\omega t}{|\mathrm{mol}⟩}_{m{m}^{\prime }}{⟨\mathrm{at}|}_{m{m}^{\prime }}+h.c\right),$$

$${|\mathrm{at}⟩}_{m{m}^{\prime }}=\frac{|\uparrow m⟩\otimes |\downarrow m^{\prime }⟩-|\downarrow m⟩\otimes |\uparrow m^{\prime }⟩}{\sqrt{2}}$$

$${\widehat{H^{\prime }}}_{\mathrm{Fesh}}=\frac{{\Omega }^2}{\delta +i\Gamma /2}\sum _{m{m}^{\prime }}{|\mathrm{at}⟩}_{m{m}^{\prime }}{⟨\mathrm{at}|}_{m{m}^{\prime }}$$

$${\widehat{H^{\prime }}}_{\mathrm{Fesh}}^f$$

$$=$$

$${\alpha }_F^f$$

$${\alpha }_F^f=|⟨F(f_1(i_1{s}_1)f_2(i_2{s}_2))|F(I(i_1{i}_2)S(s_1{s}_2))⟩{|}^2.$$

$${\alpha }_F^f=⟨F,M|P_{S=0}|F,M⟩,$$

$$|f{m}_f⟩$$

$$|FM⟩$$

$$|s\sigma ,im⟩$$

$$+\sqrt{\frac{f+1/2+\sigma m}{2i+1}}|f=i+\frac{1}{2},m_f=m+\frac{\sigma }{2}⟩,$$

$${\alpha }_F$$

$${\lambda }_m$$

$$\sum _n{\left(C_{n,-n}^{f+f\to F}\right)}^2$$

$$\sum _n{n}^2{\left(C_{n,-n}^{f+f\to F}\right)}^2$$

$$\sum _n(f-n)(f+n+1)$$

$$=⟨F,0|{\hat{\vec{f}}}_1\cdot {\hat{\vec{f}}}_2-{\hat{f}}_1^z{\hat{f}}_2^z|F,0⟩.$$

$${\alpha }_F^f=\frac{(2i+1)(2f+1)-F(F+1)}{2(2i+1{)}^2}.$$

$$S_{ab}^{\{i_1,i_2,\cdots i_L\}}=⟨A_{i_1{i}_2}{A}_{i_3,i_4}\cdots A_{i_{L-1}{i}_L}{A}_{i_2,i_3}^{†}{A}_{i_4,i_5}^{†}\cdots A_{i_L,i_1}^{†}⟩$$

$$A_{ij}=\frac{1}{\sqrt{2f+1}}\sum _m(-1{)}^{f-m}{b}_{im}{b}_{j-m}$$

$$S_{ab}=(2f+1{)}^{N_{\text{loops}}}{\left(\frac{1}{\sqrt{2f+1}}\right)}^{L_{\text{loops}}},$$

Full text

I Introduction

Quantum dimer models—which describe the dynamics of close-packed hard-core dimers on a lattice—have received continued attention since their original proposal by Rokhsar and Kivelson in 1988 Rokhsar and Kivelson (1988). Several factors have motivated these studies, including connections to Pauling and Anderson resonating valence bonds Pauling (1949); Anderson (1973), Anderson’s theory of high-$T_{c}$ superconductivity Anderson (1987); Lee et al. (2006), the appearance of quantum critical points Sachdev (1989); Senthil et al. (2004); Vishwanath et al. (2004), topological order and fractionalized excitations Kivelson et al. (1987); Wen (1991); Moessner et al. (2001a); Senthil and Fisher (2001); Fradkin et al. (2004), their mapping to lattice gauge theories Moessner et al. (2000); Moessner and Sondhi (2001); Misguich et al. (2002), and their potential applications in quantum computation Ioffe et al. (2002); Kitaev (2003). This rich variety of physics emerges due to the interplay between quantum fluctuations, hard-core constraints, and the lattice geometry of these systems. However, there are relatively few experimental realizations of dimer models. In this article, we show that experiments using an atomic gas trapped in an optical lattice can realize and probe dramatic dimer resonances in a range of quantum dimer models.

The spin physics emerging from atoms in optical lattices can be qualitatively different from those in electron systems. The atoms typically have spin greater than $1/2$, and the natural exchange processes lead to couplings which are more complicated than a simple Heisenberg model. As we demonstrate, these processes can be tuned to favor singlet dimers.

Our technique works for both fermionic and bosonic atoms with vanishing electronic orbital angular momentum ($l=0$) and relatively large hyperfine spin $f$, but with relatively weak dipole-dipole interactions. Alkalis such as 7Li or 23Na are potential candidates. By tuning the lattice depth and trapping potential one drives the system into a Mott insulating state with one atom per site Jaksch and Zoller (2005)—effectively yielding immobile spins on each site which interact via a virtual superexchange process Anderson (1950). We propose manipulating these superexchange interactions by optically coupling pairs of atoms to an excited molecular state which has $L=1$ and $S=0$. When tuned sufficiently off-resonance, this optical coupling favors the formation of nearest-neighbor hyperfine spin singlets, which we refer to as dimers. In the large-$f$ limit the dimers are monogamous and orthogonal, e.g., a state where site $i$ forms a singlet with site $j$ is orthogonal to one in which $i$ forms a singlet with $k\neq j$. The resulting theory has the form of a quantum dimer model and, depending on lattice geometry and scattering length parameters, has the potential to realize dimer crystals and dimer liquid (resonating valence bond) ground states. At smaller $f$ the dimers are not orthogonal, but nonetheless the dimer configurations span the low-energy subspace. We show how to work with this nonorthogonal basis, and derive an effective dimer model.

We numerically find the ground state of our system for small lattices. When $f$ is large we find strong dimer crystal correlations indicative of the columnar dimer state on the square lattice and the $\sqrt{12}\times\sqrt{12}$ plaquette phase on the triangular lattice Ralko et al. (2005); Leung et al. (1996) (see Sec. V). At small $f$ the results are more ambiguous, and may point towards a spin liquid or a translationally invariant symmetry broken state (such as the nematic state predicted in Ref. Rutkowski and Lawler (2016)). We explain how to further tune parameters to explore phase space—a useful requirement for the search for a spin liquid. We additionally propose a protocol to detect the dimer correlations which are central to many of these states.

To further characterize our model, we perform a large $f$ expansion, and find that as $f\to\infty$ it reduces to a special case of the Rokhsar and Kivelson model Rokhsar and Kivelson (1988). On the 2D square and 3D cubic lattice, there is some contention about the ground states of that model Banerjee et al. (2014); Sachdev (1989); Leung et al. (1996), an issue which an experimental realization of our proposal could resolve.

There have been previous proposals to observe related physics in cold atom experiments, including crystallized dimer phases Demler and Zhou (2002); Yip (2003); Imambekov et al. (2003); Zhou and Semenoff (2006); Eckert et al. (2007), resonating plaquette phases Chen et al. (2005); Xu and Wu (2008); Szirmai and Lewenstein (2011), and dimer liquid phases Hermele et al. (2009); Chan and Duan (2012); Sinkovicz et al. (2013); Song and Hermele (2013); Rutkowski and Lawler (2016)—all hallmarks of the quantum dimer model. However, these studies were generally based on different mechanisms, and did not exploit mappings of their systems onto quantum dimer models. Additionally, Ising models can be implemented in cold atoms, and such models may be mapped onto dimer models Kasteleyn (1963).

This paper is organized as follows. In Sec. II we present the system we study, and its microscopic Hamiltonian. In Sec. III we describe our proposal to tune the interactions. In Sec. IV we present the effective model describing our system. In Sec. V we numerically find the eigenstates of this model, and describe their properties. In Sec. VI we explore the large $f$ limit, mapping our system onto more traditional dimer models in Sec. VI.1, showing how to tune parameters in Sec. VI.3, and describing the phases of this model in Sec. VI.2. In Sec. VII we propose a method to observe these phases.

II Microscopic model

II.1 Setup

To produce a quantum dimer model, we begin with a tight-binding Hamiltonian for atoms in an optical lattice, which includes spin-dependent interactions. There will be a hopping term, where an atom with hyperfine spin projection $m$ moves between neighboring sites. There will also be an on-site two-particle interaction term. In the presence of rotational symmetry, these interactions can be decomposed into different angular momenta channels $F$ Stamper-Kurn and Ueda (2013). Thus, in complete generality we write

[TABLE]

where $i$ runs over all lattice sites, and $\langle ij\rangle$ runs over all distinct nearest-neighbor pairs. Due to particle statistics, we sum over only even values of $F$, up to a maximum value $F=2f$ for bosons, and $F=2f-1$ for fermions. The $\hat{b}^{\dagger}_{i,m}(\hat{b}_{i,m})$ operators create (annihilate) an atom at lattice site $i$ with hyperfine spin $f$ and spin projection $m$, while the $\hat{A}_{ij}^{F,M\dagger}(\hat{A}_{ij}^{F,M})$ operators create (annihilate) a pair of atoms on sites $i$ and $j$ in total angular momentum state $F$ with total projection $M$. These operators may be defined via the relation

[TABLE]

Here, $C^{F,M}_{m,m^{\prime}}=C^{f+f\to F}_{mm^{\prime}}=\langle f,m;f,m^{\prime}|F,M\rangle$ are Clebsch-Gordan coefficients, and the factor of $\sqrt{2}$ is chosen so that $\langle\hat{A}_{ij}^{FM}\hat{A}_{ij}^{FM\dagger}\rangle=1$ in the vacuum state. The kinetic energy term in Eq. (1)—parameterized by the positive constant $J$—models the tunneling of atoms between neighboring lattice sites. The parameters $U_{F}$ encode the local spin-dependent interactions. While typically one expects that the scattering in different spin channels to be of similar magnitude, in Sec. III we argue that one can engineer an optical Feshbach resonance so that the interactions are significantly weaker in the hyperfine singlet channel than all others: $U_{F\neq 0}\gg U_{0}>0.$

We refer readers to the review article by Stamper-Kurn and Ueda Stamper-Kurn and Ueda (2013), for further background on Eq. (1).

II.2 Effective nearest-neighbor inearest-neighbornteraction

In the limit where the interactions are strong compared to the hopping ($U_{F}\gg J$), and there is exactly one particle per site, this system should form a Mott Insulator. Super-exchange will lead to a magnetic coupling between neighboring sites. In particular, we let $P$ be the projector into the space with one particle per site, and define $H_{0}=PHP$, $\Lambda=(1-P)HP$, $\Lambda^{\dagger}=PH(1-P)$, and $H_{1}=(1-P)H(1-P)$. Here $H_{0}=0$, and $\Lambda\propto J$ is considered small. We consider an eigenstate $\psi$, with $\psi_{0}=P\psi$, and $\phi=(1-P)\psi$. The Schrodinger equation $H\psi=E\psi$ can be projected into the space with one particle per site, and into the complementary space to give $H_{0}\psi_{0}+\Lambda^{\dagger}\phi=E\psi_{0}$ and $H_{1}\phi+\Lambda\psi_{0}=E\phi$. To lowest order in $J$, the second equation yields $\phi=-H_{1}^{-1}\Lambda\psi_{0}+{\cal O}(J)$, and hence $(H_{0}-\Lambda^{\dagger}H_{1}^{-1}\Lambda)\psi_{0}=E\psi_{0}+{\cal O}(J^{3})$, which yields the effective Hamiltonian $H_{\rm eff}=H_{0}-\Lambda^{\dagger}H^{-1}\Lambda$, or explicitly

[TABLE]

Under the condition $U_{F\neq 0}\gg U_{0}>0$, we can neglect all but the $F=0$ term to find

[TABLE]

On bipartite lattices this model is an example of a $SU(N)$ antiferromagnet model Read and Sachdev (1989); Harada et al. (2003). We show in Sec. IV that our model [Eq. (4)] can be mapped onto a dimer model.

For typical parameters (lattice depth $V_{x}=V_{y}=10E_{R},V_{z}=30E_{R}$, wavelength $\lambda=1064\ {\rm nm}$, and scattering length $a_{0}=20\ {\rm Bohr}$), the superexchange coefficient is $\frac{2J^{2}}{U_{0}}=200\ {\rm Hz}$. This scale is large compared to neglected physics such as off-site dipole interactions ($\sim 0.5{\rm Hz}$ for alkali atoms). One may also worry about tensor light shifts from the lattice or photoassociation beams. These will be minimal if the laser detunings are larger than hyperfine splitting Grimm et al. (2000). If any residual light shifts remain, they can be canceled by adding additional fields (e.g., as in Ref. Gerbier et al. (2006)).

III Tuning the interactions via an Optical Feshbach Resonance

We propose inducing an optical Feshbach resonance Fedichev et al. (1996); Bohn and Julienne (1997); Ciuryło et al. (2005); Nicholson et al. (2015); Miller et al. (1993); Abraham et al. (1995, 1996); Junker et al. (2008); Hamley et al. (2009); Gerton et al. (2001); Stwalley and Wang (1999); Jelassi et al. (2006a, b); Kemmann et al. (2004); Wang et al. (1996, 2000); Pichler et al. (2003); Fatemi et al. (2000); Enomoto et al. (2008); Yamazaki et al. (2010); Theis et al. (2004); Yan et al. (2013); Blatt et al. (2011); Kim et al. (2016) between pairs of atoms by shining a laser tuned near a transition to an excited molecular state, labeled by orbital angular momentum $L=1$, electronic spin $S=0$, and total electronic angular momentum $J_{a}=1$. These are good quantum numbers in molecules formed from lighter elements such as Lithium or Sodium, where spin-orbit coupling is relatively weak [Hund’s case (b)] Brown and Carrington (2003). For example, the laser can be tuned to couple the atoms to ${}^{1}\Sigma_{u/g}$ molecular states, as in Refs. Abraham et al. (1996) and Fatemi et al. (2000). The nuclear angular momentum is not important as long as the detuning of the laser is large compared to the hyperfine splitting. In the cold collision limit the rotational angular momentum of the nuclei vanishes, $R=0$. As in Ref. Sundar and Mueller (2016), second-order perturbation theory then gives a contribution to $U_{F}=U_{F}^{\rm bg}+U_{F}^{\rm Fesh}$ of

[TABLE]

with $U_{F}^{\rm bg}$ encoding the background scattering, including any influence of on-site dipole-dipole interactions. The matrix element $\Omega^{2}$ is proportional to the intensity of the laser. The detuning $|\delta|$ must be taken much larger than the linewidth $\Gamma$, so molecular decay can be neglected Hamley et al. (2009). This limit is reasonable, as in typical experiments $\delta\sim$GHz and $\Gamma\sim$100 MHz Hamley et al. (2009). The line must be chosen judiciously so that unwanted transitions are avoided.

The coefficient $\alpha_{F}$ is the square of the overlap between the electronic spin singlet ($S=0$) and the hyperfine state with total spin $F$ Messiah (1961): It is related to the Wigner $9j$ symbols. In Appendix A we show how to calculate this matrix element for the most relevant case of alkali atoms, which have $s=1/2$ and nuclear spin $i$. We find

[TABLE]

Assuming that $U_{F}^{\rm bg}$ depends only weakly on $F$, one can then choose the laser intensity and detuning so that $U_{F\neq 0}\gg U_{0}\gg\left|\frac{\Gamma}{\delta}U_{F}^{\rm Fesh}\right|$, which then yields Eq. (4). For Lithium, the laser intensity required to achieve this limit is only a few ${\rm W/cm}^{2}$.

IV From Singlet Coverings to Dimer Models

The Hamiltonian in Eq. (4) appears to count nearest-neighbor singlet bonds. One might therefore expect that the ground state would be formed by creating some pattern of nearest-neighbor singlets, which we will describe as a “singlet covering.” For example, the first image in Fig. 1 illustrates one possible singlet covering of six sites that are laid out in a rectangle: $\left|a\right\rangle=\hat{A}^{\dagger}_{1,2}\hat{A}^{\dagger}_{4,5}\hat{A}^{\dagger}_{3,6}\left|0\right\rangle$, where $\left|0\right\rangle$ is the vacuum state with no particles. The label $a=\{(1,2),(4,5),(3,6)\}$ is the set of all bonds. In general,

[TABLE]

This definition works even for coverings which involve longer-range bonds (such as the third image in Fig. 1). The operator $\hat{A}^{\dagger}_{i,j}$ adds exactly one particle each to sites $i$ and $j$. We are working in the sector with exactly one particle per site, and therefore we require each site to appear in one and only one of the bonds. Furthermore, by the standard rules of adding angular momentum, a particle cannot be in a singlet bond with more than one other particle.

As explained by Rokhsar and Kivelson Rokhsar and Kivelson (1988) in the context of spin-1/2 electrons, the singlet coverings are not eigenstates of the Hamiltonian, but they are closed under the action of Eq. (4), and the ground state is a superposition of such coverings. In particular, a single term in the Hamiltonian maps one single covering into another:

[TABLE]

where the notation $\left|(i,j):a\right\rangle$ denotes a state where sites $i$ and $j$ are paired together into a singlet, the original partners of $i$ and $j$ in $\left|a\right\rangle$ are paired together into a singlet, and all the other bonds in $\left|a\right\rangle$ are left unchanged. An example of a singlet covering $\left|a\right\rangle$ and a few of the related states $\left|(i,j):a\right\rangle$ are illustrated in Fig. 1. Note, the labeling is not unique: One goes from $\left|a\right\rangle$ in Fig. 1 to $\left|(1,4):a\right\rangle=\left|(2,5):a\right\rangle$ by either acting with $\hat{A}_{14}^{\dagger}\hat{A}_{14}^{\phantom{\dagger}}$ or $\hat{A}_{25}^{\dagger}\hat{A}_{25}^{\phantom{\dagger}}$. As illustrated by the right-most figure, the nearest-neighbor bond operators acting on a state with nearest-neighbor bonds can generate configurations with longer-ranged bonds.

Thus, the Hamiltonian is a map on the space of singlet coverings. Somewhat complicating the analysis, however, is the fact that the singlet coverings are not orthogonal. In fact, they are not even linearly independent. Nonetheless, as detailed below, it is straightforward to work with these states. In Appendix B we review the more traditional approach of orthogonalizing the states, which is somewhat more involved. Working with either basis gives equivalent results.

We consider a general state $|\psi\rangle=\sum_{a}\psi_{a}|a\rangle$. Equation (8) allows us to write

[TABLE]

Clearly if $\sum_{a}H_{ba}\psi_{a}=E\psi_{b}$, then $\hat{H}|\psi\rangle=E|\psi\rangle$, and the eigenstates of the matrix $H_{ba}$ either yield eigenstates of $\hat{H}$, or are null-vectors. The latter have eigenvalue 0. We will solely be concerned with states with negative energy, and hence will not encounter any of these null-states. Because of the nonorthogonality, $H_{ba}$ is not a Hermitian matrix—but it is self-adjoint with respect to the natural inner product: $\langle\phi|\psi\rangle=\sum_{ab}\phi_{a}^{*}\psi_{b}S_{ab}$, with $S_{ab}=\langle a|b\rangle$.

In Sec. V, we consider a small system, and enumerate all dimer coverings. Although exponentially large in the system size, this is a much smaller Hilbert space than a spin model on the same lattice. We then numerically calculate $H_{ab}$, and find its eigenstates.

If we formally set $f=\infty$, then the dimer coverings are orthogonal, and become eigenstates of the Hamiltonian. The energies of these states are negative, and proportional to the number of nearest-neighbor dimers. Thus the ground state manifold is highly degenerate, consisting of all nearest-neighbor coverings. In Sec. VI we derive a systematic expansion in $1/(2f+1)$, and find that the leading terms break this degeneracy and stabilize various dimer crystal or plaquette phases. The structure of this expansion corresponds to Rokhsar and Kivelson’s model Rokhsar and Kivelson (1988).

In the limit of small $f$, we anticipate the dimer crystal order to disappear. In particular, Rutkowski and Lawler Rutkowski and Lawler (2016), carried out a variational study of the Hamiltonian in Eq. (3), and argued that for $f<3$, a translationally invariant nematically ordered state will be found. Another, even more exciting possibility is a spin liquid—which could either occur as an intermediate phase, or at the phase boundary.

Related physics is seen in studies of anisotropic 2D lattices of coupled spin-1 objects Harada et al. (2007); Grover and Senthil (2007). In those studies, the researchers finds regions with nematic order, and others with dimer crystal order. Surprisingly, there appears to be a direct second-order phase transition between these phases: Within the Landau paradigm such a direct transition would require fine-tuning. Moreover there is evidence that this transition displays “deconfined quantum criticality,” where the transition is described by an emergent gauge theory, and is a spin liquid Senthil et al. (2004). By analogy, one might expect that our model would display a similar critical point as $f$ is changed. We have, however, not yet verified this conjecture.

V Numerical Results

We numerically diagonalize the matrix $H_{ab}$ in Eq. (9). To visualize the ground state, we calculate the correlation functions $C_{ij}^{kl}=\langle\hat{A}_{ij}^{\dagger}\hat{A}_{ij}^{\phantom{\dagger}}\hat{A}_{kl}^{\dagger}\hat{A}_{kl}^{\phantom{\dagger}}\rangle$. As we argue in Sec. VII this is an experimental observable. It corresponds to the probability of simultaneously having a dimer on bonds $(i,j)$ and $(k,l)$. In a state $|a\rangle$ consisting of a single dimer covering, the expectation value $\langle a|\hat{A}_{ij}^{\dagger}\hat{A}_{ij}^{\phantom{\dagger}}|a\rangle$ is equal to $1$ if $(i,j)\in a$ and $1/(2f+1)^{2}$ otherwise. Similarly, $\langle a|\hat{A}_{ij}^{\dagger}\hat{A}_{ij}^{\phantom{\dagger}}\hat{A}_{kl}^{\dagger}\hat{A}_{kl}^{\phantom{\dagger}}|a\rangle$ is equal to $1$ if both $(i,j)\in a$ and $(k,l)\in a$, and is otherwise suppressed by factors of $1/(2f+1)^{2}$.

Figure 2 shows the calculated correlations between a horizontal bond $(i,j)$ and a bond $(k,l)$ in the ground state on a square lattice for two different values of $f$, where the bond $(i,j)$ is fixed in the lower left corner, and the bond $(k,l)$ is varied. The fixed bond is colored red, while the other is colored based upon the strength of the correlations. A distinctive “ladder” pattern can be observed at $f=100$. This bond configuration is characteristic of the “columnar state.” Our large $f$ expansion in Sec. VI indeed confirms that the columnar state is expected to be the ground state at large $f$. In a finite size sample, the idealized columnar state is a quantum superposition of four symmetry related dimer crystals, as shown in Fig. 3. Within this cartoon, and taking $f\to\infty$, all correlations $C_{ij}^{kl}=\langle\hat{A}_{ij}^{\dagger}\hat{A}_{ij}\hat{A}_{kl}^{\dagger}\hat{A}_{kl}\rangle$ will be either 0 or 0.25, depending on the two bonds $(i,j)$ and $(k,l)$ are both found in the same configuration. The reference value should be $C^{ij}_{ij}=0.25$. The correlations in Fig. 2(a) share this same pattern, but the contrast is somewhat weaker than in the idealized picture. Such “quantum fluctuations” are due to the fact that the quantum state has weight on configurations other than those given by this cartoon.

For $f=3$ [Fig. 2b], the correlations are somewhat more ambiguous. The pattern includes short-range columnar order, but it is unclear if there is long-range order.

Figure 4 shows the correlations on a triangular lattice for the same two values of $f$. For this lattice, we anticipate the ground state to be the $\sqrt{12}\times\sqrt{12}$ phase Ralko et al. (2005), and so we take our system to have the shape of a unit cell in the $\sqrt{12}\times\sqrt{12}$ phase. As illustrated in Fig. 5(b), the expected unit cell consists of 12 sites, which resonate between two plaquette configurations. Figure 5(a) further illustrates that in each plaquette the spins are expected to resonate between two different dimer configurations. The correlations corresponding to this ansatz are shown in Figs. 5(c) and  5(d), for two different reference bonds. The thicker (blue) bonds in the idealized model have $C=0.25$, while thinner (green) bonds have $C=0.125$, as labeled. The interpretation is that when one expands out the superposition in Fig. 5(b), 1/4 of the terms will simultaneously have bonds at a given red and blue position, and 1/8 of the terms will simultaneously have bonds at a given red and green position. The red bond is thicker in Fig. 5(d) than 5(c), as it appears in more terms of the superposition.

For large $f$, the pattern of bonds in Fig. 4 are nearly identical to what one expects from the $\sqrt{12}\times\sqrt{12}$ phase. The deviations are of comparable size to those seen in the square lattice. The $f=3$ pattern shares some of the same symmetries, but one observes significant differences, which can be interpreted as spatial broadening. For example, in the lower image of Fig. 4, the central three bonds have relatively substantial weight, while no such weight is found in the cartoon of the $\sqrt{12}\times\sqrt{12}$ phase. Similar discrepancies were seen in the numerical studies of the models in Ref. Ralko et al. (2005). Since our exact diagonalization approach is only able to capture a single unit cell, we cannot say anything about long-range order from this calculation.

We also use these numerical results to investigate the stability of our system against small perturbations, such as off-site dipole-dipole interactions and magnetic field noise. Using parameters appropriate for Lithium, we find that the amplitude for off-site dipole interactions to create an excitation is only $\sim 0.008\ {\rm Hz}$, which is small compared to the lowest energy plaquette-flip excitation $\sim 0.9J^{2}/U_{0}\sim 90\ {\rm Hz}$. Our system is however sensitive to small perturbations in the magnetic field, which will break dimers to align spins in its direction. In order to remain in the ground state, the magnetic field in the experiment should be less than $30\mu$G.

VI Large $f$ limit

VI.1 Mapping onto the Rokhsar-Kivelson model

We now consider the large $f$ limit, and show that our model maps onto the classic Rokhsar-Kivelson dimer model.

For notational convenience it is useful to introduce a fictitious Hilbert space with orthonormal basis states labeled by the singlet coverings $|\tilde{a}\rangle$. This allows us to use the familiar language of bras and kets in manipulating $H_{ab}$. For example, we write $\tilde{H}=\sum_{ab}H_{ab}|\tilde{a}\rangle\langle\tilde{b}|$. An eigenstate $|\tilde{\psi}\rangle$ of $\tilde{H}$ can be mapped into the physical space via $|\psi\rangle=\hat{P}|\tilde{\psi}\rangle$, where $\hat{P}=\sum_{a}|a\rangle\langle\tilde{a}|$. All eigenvalues of $\tilde{H}$ are also eigenvalues of $H$, with the caveat that zero energy eigenvectors may be unphysical as quantum states.

The diagonal elements of $\tilde{H}$, $H_{aa}=(-2J^{2}/U_{0})N_{a}$ count the number of nearest-neighbor singlets in $|a\rangle$. The low-energy space is then spanned by nearest-neighbor singlet coverings. We will elimate the other modes to derive an effective Hamiltonian which acts only in this low-energy space. The key point is that the off-diagonal matrix elements are of order $(2f+1)^{-1}$, and are small in the limit $f\to\infty$.

The leading order term in the effective Hamiltonian comes from the parts of $\hat{H}$ which directly take one between nearest-neighbor singlet coverings. For example, the $\hat{A}_{14}^{\dagger}\hat{A}_{14}$ or $\hat{A}_{25}^{\dagger}\hat{A}_{25}$ terms acting on the state $|a\rangle$ in Fig. 1. Terms of this form take two parallel vertical nearest-neighbor bonds, and replaces them with horizontal bonds (or vice versa), and can be represented as $\tilde{H}_{0}=-t\left(\left|\tilde{=}\right\rangle\left\langle\tilde{\shortparallel}\right|+\left|\tilde{\shortparallel}\right\rangle\left\langle\tilde{=}\right|\right)$ with $t=(2J^{2}/U_{0})(2/(2f+1))$. The factor of 2 comes from the fact that the same term is generated if one acts on either of the new bonds.

Acting on any of the other nearest-neighbor bonds will introduce a long-range bond. For example, the $\hat{A}_{23}^{\dagger}\hat{A}_{23}$ term acting on the state $|a\rangle$ in Fig. 1, yields the state with long range bonds on the far right. There are $3N/2-2(N_{=}+N_{\shortparallel})$ such terms, each of which will contribute to the effective Hamiltonian in second-order perturbation theory. Here $N_{=}$ and $N_{\shortparallel}$ count the number of plaquettes with two horizontal or vertical bonds. As illustrated in Fig. 6(a), for each of these terms, there are two ways to return to the initial state. Up to an additive constant, one thus finds a contribution to the effective Hamiltonian of $\tilde{H}_{1a}=V\left(\left|\tilde{=}\right\rangle\left\langle\tilde{=}\right|+\left|\tilde{\shortparallel}\right\rangle\left\langle\tilde{\shortparallel}\right|\right)$ with $V=2\times\lambda\times(2\lambda)/\epsilon=(8J^{2}/U_{0})/(2f+1)^{2}$. Here one factor of 2 comes from bond counting, one $\lambda=(2J^{2}/U_{0})/(2f+1)$ is for the forward matrix element, the $2\lambda$ is for the backward matrix element, and $\epsilon=(2J^{2}/U_{0})$ is the energy denominator.

There are also second order processes which, as illustrated in Fig. 6(b), rotates a set of three bonds, and can be written as $\hat{H}_{1b}=-t^{\prime}\left(\left|\widetilde{=\shortmid}\right\rangle\left\langle\widetilde{\shortmid=}\right|+\rm{h.c.}\right)$, with $t^{\prime}=(2\lambda)^{2}/\epsilon_{0}=(8J^{2}/U_{0})/(2f+1)^{2}$. Here one factor of $2$ comes from the existence of two possible intermediate states (each of which can be reached in a single way). The second factor of $2$ comes from the two ways to reach the rotated configuration from each intermediate state.

This reasoning can be continued to generate terms involving longer and longer ring exchanges. Any term which appears at $m^{\rm th}$ order scales as $(2f+1)^{-m}$. Thus, unlike the spin-1/2 electronic case, larger ring exchange terms are strongly suppressed for large $f$.

This same argument goes through on any lattice. On a square, cubic, or triangular lattice, the effective Hamiltonian is of the form

[TABLE]

This defines the matrix elements $H_{ab}$. Table 1 lists the parameters $t,t^{\prime}$, and $V$ for different lattice geometries. On the triangular lattice one interprets parallel nearest-neighbor bonds as those that are on opposite sides of a rhombus made from two triangular units. The effective dimer model on the honeycomb and Kagome lattices have similar terms, but the smallest kinetic term involves three-bond loops, and therefore $t$ has an amplitude of $\mathcal{O}\left(2f+1\right)^{-2}$. Similarly, on those lattices the potential term $V\sim\mathcal{O}\left(2f+1\right)^{-4}$ penalizes parallel bonds on alternate sides of a hexagon, while $t^{\prime}\sim\mathcal{O}\left(2f+1\right)^{-4}$ involves a ring with bonds extending over two hexagons.

VI.2 Phases

The effective model for our system from Eq. (VI.1) has a rich phase diagram, which has been well explored along $t^{\prime}=0$ in a number of geometries Moessner and Raman (2011); Sachdev (1989); Moessner et al. (2001a); Moessner and Sondhi (2001); Moessner et al. (2001b); Moessner and Sondhi (2002); Misguich et al. (2002); Ralko et al. (2005); Zeng and Elser (1995); Leung et al. (1996); Syljuåsen (2006); Ralko et al. (2008); Banerjee et al. (2014). For 2D bipartite lattices with $t^{\prime}=0$, one finds only valence bond solid phases, except for the Rokhsar-Kivelson point $V=t$. On 3D and nonbipartite 2D lattices, dimer liquids may be found for nonvanishing ranges of $t/V$. The phase diagram at finite $t^{\prime}$ is less explored Nakata and Totsuka (2011).

The valence bond solid phases described in the literature fall into four types: columnar, plaquette, mixed, and staggered. The columnar phase is built from vertical columns of horizontal parallel bonds, or vice versa. In the plaquette phase, dimer bonds resonate between different configurations inside a multi-site unit cell. For example, on a square lattice, the plaquette phase has a unit cell with four lattice sites; two parallel bonds resonate between horizontal and vertical configurations inside a plaquette. The plaquette phases on a triangular lattice have larger unit cells. The mixed phase is a hybrid between the columnar and plaquette phases, which is best described in terms of the symmetries it breaks Banerjee et al. (2014). The staggered phase has no flippable plaquettes ($=$ or $\shortparallel$). The columnar phase is favored at large negative $V$, and the staggered phase at large positive $V$.

As $f\to\infty$, the dominant coupling constant in the effective model, Eq. (VI.1), is $t$. On a square lattice, this generally is believed to lead to a columnar phase, though there is some contention Banerjee et al. (2014). (Experiments may be able to resolve these issues.) On the triangular lattice, as $f\rightarrow\infty$, we expect to see a plaquette phase, called the $\sqrt{12}\times\sqrt{12}$ phase, which has a $12$-site unit cell, and quantum resonances that extend throughout the cell Ralko et al. (2005). Observing these resonances is part of Pauling and Anderson’s vision of quantum resonances that manifest throughout a macroscopic system Anderson (1973). The analysis in Sec. V confirms that at $f=100$ these orders appear to be present. The smaller $f$ data is more ambiguous, and could point towards a spin liquid or some other phase.

VI.3 Tuning the parameters in the large $f$ dimer model

As presented, a given experimental realization yields a unique dimer model: Aside from the overall energy scale, all parameters are determined by the spin $f$ and the lattice geometry. We can gain the ability to tune the parameters by modifying the detuning and coupling strength of the optical Feshbach resonance. For example, consider the case $U_{F\neq 0,2}\gg U_{2}>U_{0}\gg J$. Then the effective Hamiltonian at second order in the tunneling strength is

[TABLE]

where $\hat{A}_{ij}^{2M\dagger}=\sum_{m}C_{m,M-m}^{2M}\hat{b}_{i,m}\hat{b}_{j,M-m}$ creates a neighboring atom pair with total spin $F=2$ and azimuthal spin $M$. In Appendix C we use our large-$f$ perturbation techniques to obtain a dimer model of the form of Eq. (VI.1). We find that to leading order $t$ and $t^{\prime}$ are independent of $U_{2}$, but $V$ depends on $U_{2}$. By tuning $U_{2}/U_{0}$ via the Feshbach laser intensity and detuning, one can control the relative size of $V$.

VII Detection

To probe the valence bond solid order and observe the resonating dimers in the plaquette and spin liquid phase, we propose measuring the dimer-dimer correlation function $\langle\hat{A}_{ij}^{\dagger}\hat{A}_{ij}^{\phantom{\dagger}}\hat{A}_{kl}^{\dagger}\hat{A}_{kl}^{\phantom{\dagger}}\rangle$. Similar correlation functions have been used to characterize order in quantum dimer models Misguich et al. (2002); Ralko et al. (2005); Zeng and Elser (1995); Leung et al. (1996). We provide an experimental protocol to image these correlations. Furthermore, in Sec. V we numerically calculated these correlations in our system for both large and moderate values of $f$.

To image the dimer bonds we propose shining a weak near-resonant photoassociation laser on the system, tuned near a molecular state with angular momenta $L=1$ and $S=0$. In our system, when virtual hopping brings two atoms forming an $S=0,L=0$ dimer onto the same lattice site, the near-resonant light drives these atoms into the molecular state. The excited molecule has a short lifetime and so those atoms are lost from the trap.

After driving this photoassociation, one would use a quantum gas microscope to image the location of all remaining atoms Bakr et al. (2009); Jones et al. (2006). All adjacent pairs of empty sites in the image were likely occupied by atoms entangled in dimers. In this way, a fraction of the dimers in the system can be imaged. Quantitative dimer-dimer correlations can be extracted by analyzing data from multiple realizations of this imaging process, and can be used to identify the phase. Similar techniques have been used in the past to probe atomic correlations Partridge et al. (2005).

One formal way to model this process is to take $U_{0}\to U_{0}+i\Gamma/2$, where $\Gamma$ quantifies the photoassociation rate. We thus see that the Hamiltonian in Eq. (4) gains an imaginary term which removes a pair of neighboring particles. The probability that (after a fixed time) atoms at neighboring sites $i$ and $j$ are missing will be proportional to $\langle\hat{A}_{ij}^{\dagger}\hat{A}_{ij}^{\phantom{\dagger}}\rangle$. The probability that there are also missing atoms at neighboring sites $k$ and $\ell$ will then be proportional to $\langle\hat{A}_{ij}^{\dagger}\hat{A}_{ij}^{\phantom{\dagger}}\hat{A}_{kl}^{\dagger}\hat{A}_{kl}^{\phantom{\dagger}}\rangle$.

We emphasize that the ability to directly image the valence-bond correlations is one of the greatest strengths of using cold atoms to explore dimer models. This imaging will allow unambiguous identification of the various valence-bond ordered phases. Spin liquid phases will be characterized by the absence of long-range valence bond order. The experimental systems are much larger than those we can model numerically.

VIII Summary

In summary, we propose experimental protocols to produce quantum dimer models and detect both static and resonating patterns of dimer configurations. In particular, we show that appropriately tuned off-resonant photoassociation light modifies the interactions in a gas of cold atoms, yielding a low-energy Hilbert space spanned by short-range dimers. By expanding in powers of $(2f+1)^{-1}$ we develop an effective dimer model Hamiltonian, and discuss its phase diagram. We find that a number of valence bond solid and plaquette phases are readily produced, and suggest techniques which are suited to searching for even more exotic states such as topological spin liquids. We demonstrate that by combining photoassociation with quantum gas microscopy one can directly detect the dimers and the dimer-dimer correlations, thereby probing the defining features of these phases. We numerically calculate the dimer correlations, finding that on triangular lattices one will be able to image an intricate pattern of resonating bonds, extending over a 12-site unit cell.

Quantum dimer models have been highly influential in developing an understanding of how geometric constraints lead to new emergent physics Moessner and Sondhi (2002); Moessner and Raman (2011), and they have been used as a theoretical foundation for attempting to understand phenomena ranging from high temperature superconductivity to exotic antiferromagnets Savary and Balents (2017). A direct experimental realization of dimer models is key to validating and refining these ideas.

Acknowledgment

This material is based upon work supported by the NSF Grants No. PHY-1508300, No. PHY-1806357, and No. PHY11-25915. This work was supported in part by the Data Analysis and Visualization Cyberinfrastructure funded by NSF under Grant No. OCI-0959097. B.S. acknowledges useful discussions with Randall Hulet, Anna Marchant and Jacob Fry. B.S. and T.C.R. contributed equally to this work.

Appendix A Deriving the effective interaction due to an optical Feshbach resonance

Here we derive the effective interaction induced by our proposed optical Feshbach resonance. We closely follow the argument in our previous work, Ref. Sundar and Mueller (2016). A laser is tuned near a transition to a molecule state with well-defined electronic spin $S=0$, and well-defined electronic angular momentum $L=1$. Keeping only the relevant degrees of freedom, and neglecting any coupling to the nuclear degrees of freedom, we model the photoassociation as

[TABLE]

where the electronic singlet state is

[TABLE]

Here, $\uparrow/\downarrow$ represent the spin projection $s_{z}$ of the spin-1/2 electrons, while $m,m^{\prime}$ are the spin projections of the nuclei. Due to hyperfine interactions, $\left|\rm at\right\rangle_{mm^{\prime}}$ is not an eigenstate of the atomic Hamiltonian.

The energy of the molecule is $E$, and we have included an imaginary part, $\Gamma$ to model its finite lifetime. In principle, the molecular energy should have some dependence on the nuclear spin projections, but these will play no role as long as the detuning of the laser is large compared to the hyperfine splitting. Thus, we ignore them. We further assume that the atoms are both in the same spatial mode of a single site of the optical lattice, and therefore drop spatial indices. The coupling $\Omega$ will depend on the shape of this mode. The laser frequency $\omega=E+\delta$ is detuned from the atom-molecule transition by $\delta$. When the detuning $\delta$ is large compared to $\Omega$, we use second-order perturbation theory to eliminate the molecule and obtain in a rotating frame

[TABLE]

If the incoming and outgoing atoms are restricted to being in a single hyperfine manifold ($f=i\pm 1/2$), then symmetry implies that this expression can be replaced by

[TABLE]

where $|F,M\rangle_{f}$ is the two-particle state with total hyperfine spin $F$ and total spin projection $M$, built from two particles in the manifold with hyperfine spin $f$. The operator $P_{f}$ projects into the space where each atom has spin projection $f$. The $SU(2)$ symmetry implies that the coefficients $\alpha_{F}^{f}$ do not depend on $M$. In the main text we do not explicitly write the $f$ labels. Equating these expressions for $\hat{H}^{\prime}_{\rm Fesh}$ gives

[TABLE]

where the state $\left|\rm at\right\rangle_{m(M-m)}$ is given by Eq. (13) with $m^{\prime}=M-m$. Alternatively, this can be written as the square overlap of two states: the first is formed by combining $i_{1}$ and $s_{1}$ into $f_{1}$, $i_{2}$ and $s_{2}$ into $f_{2}$, then $f_{1}$ and $f_{2}$ into $F$. The second is formed by combining $s_{1}$ and $s_{2}$ into $S$, $i_{1}$ and $i_{2}$ into $I$, then $S$ and $I$ into $F$. The nine angular momenta $s_{1},s_{2},i_{1},i_{2},f_{1},f_{2},S,I,F$ can be combined into a Wigner 9-j symbol DLMF . The most natural notation for this construction involves recursively noting how each angular momentum is constructed. For example $f_{1}(i_{1}s_{1})$ indicates that $f_{1}$ is built from $i_{1}$ and $s_{1}$. In this notation,

[TABLE]

Here we have a particularly simple case where $s_{1}=1/2,s_{2}=1/2,i_{1}=i,i_{2}=i,f_{1}=f,f_{2}=f,S=0,I=F$. A third representation of the coefficient is the expectation value

[TABLE]

where $P_{S=0}=(1/4)-{\bf S_{1}\cdot S_{2}}$ is the projector into the space where ${\bf S}={\bf S_{1}+S_{2}}=0$.

There are several ways to evaluate $\alpha_{F}$. The simplest is to note that the condition $S=0$ reduces Eq. (18) to a 6-j symbol—which is tabulated in Ref. DLMF for the case $s_{1}=s_{2}=1/2$. The second is to directly evaluate Eq. (17). Presumably there is also an approach based upon Eq. (19). In the remainder of this section we outline the second method, based on Eq. (17).

We first note that since the result is independent of $M$, we can set $M=0$. We then find a common basis for each set of states, using the Clebsch-Gordon coefficients, defined by

[TABLE]

Using tabulated expressions for $S=1/2$, we can invert this relationship to arrive at

[TABLE]

where $\sigma=+1(-1)$ corresponds to $\uparrow(\downarrow)$, and as before $m$ is the nuclear spin projection. Substituting this result into Eq. (13) and combining it with Eq. (21) yields

[TABLE]

We then derive a series of sum rules: First, we express $A=\sum_{n}\left(C^{f+f\to F}_{n,-n}\right)^{2}$, as $A=\sum_{n}\langle F,0|f,n;f-n\rangle\langle f,n;f-n|F,0\rangle$. This sum contains a resolution of the identity in the sector with $M=0$, and hence $A=\langle F,0|F,0\rangle=1$. Second, by the same reasoning $B=\sum_{n}n^{2}\left(C^{f+f\to F}_{n,-n}\right)^{2}$ $=\sum_{n}\langle F,0|f,n;f,{-n}\rangle n^{2}\langle f,n;f,-n|F,0\rangle$ which can be identified as the expectation value $B=-\langle F,0|\hat{f}_{1}^{z}\hat{f}_{2}^{z}|F,0\rangle$. Finally, $C=\sum_{n}(f-n)(f+n+1)\left(C^{f+f\to F}_{n,-n}\right)^{2}$ $=\langle F,0|\hat{f}_{1}^{+}\hat{f}_{2}^{-}|F,0\rangle$. Using the symmetry between the two spins, we can simplify this to $C=\langle F,0|\hat{\vec{f}}_{1}\cdot\hat{\vec{f}}_{2}-\hat{f}_{1}^{z}\hat{f}_{2}^{z}|F,0\rangle$. The resulting three identities are

[TABLE]

These sum rules, plus the expression $\langle F,0|\hat{F}^{2}|F,0\rangle=F(F+1)=2f(f+1)+2\langle F,0|\hat{\vec{f}}_{1}\cdot\hat{\vec{f}}_{2}|F,0\rangle$, allow us to write Eq. (A) as

[TABLE]

For our purposes, the most important feature of this expression is that it is monotonic in $F$.

Appendix B Construction of an orthogonal basis

In the main text we work with a nonorthogonal basis, while the more traditional approach involves orthogonalizing the basis states, as originally developed for spin-$1/2$ systems in Ref. Rokhsar and Kivelson (1988). Here we follow a similar procedure to construct the orthogonal basis of dimer states for cold atoms for large spin $f$, perturbatively in $(2f+1)^{-1}$. The singlet coverings contain both short- and long-ranged singlet bonds.

We first note that the singlet operators commute with one-another unless they share a site. Hence the overlap between two singlet coverings $S_{ab}=\langle a|b\rangle$ factors into expectation values of loops: sets of connected sites $\{i_{1},i_{2},\cdots i_{L}\}$. We can always label these loops so that $(i_{2j+1},i_{2j+2})\in a$ and $(i_{2j+2},i_{2j+3})\in b$ for $j=0,1,\cdots L/2-1$. We will also have $(i_{L},i_{1})\in b$. The contribution to $S_{ab}$ from such a loop is

[TABLE]

where the expectation value is in the vacuum state with no particles. The full $S_{ab}$ is the product of the contribution from all such loops. To graphically generate this set of loops, one simply takes the set of all bonds which are in only one of $a$ and $b$, but not the other. Figure 7 shows two examples of such a graphical construction of the overlap matrix elements for a six-site system.

Using the expression for the Clebsch-Gordan coefficients, and assuming $i\neq j$, Eq. (2) becomes

[TABLE]

We substitute Eq. (27) into Eq. (26), and use Wick’s theorem to evaluate the vacuum expectation value. There is only one nonzero contraction, as there is only one creation operator and only one annihilation operator acting on each site. Once the $m$ of a single site is set, all others are fixed. There are $2f+1$ choices for $m$, and each term contributes equally. Hence $S_{ab}^{\{i_{1},i_{2},\cdots i_{L}\}}=(2f+1)^{1-L/2}$.

The full expression for $S_{ab}$ is just the product of the contribution from each loop, and hence

[TABLE]

where $N_{\text{loops}}$ is the total number of closed loops formed by the dimers not common to $\left|a\right\rangle$ and $\left|b\right\rangle$, while $L_{\text{loops}}$ is the total number of sites involved in all loops.

For large $f$, we expand $S_{ab}$ in powers of $(2f+1)^{-1}$ to obtain

[TABLE]

Here, $\Box_{ab}=1$ if $\left|a\right\rangle$ and $\left|b\right\rangle$ differ by a four-site loop in their transition graph, and is zero otherwise. The sites making up the loops do not need to be nearest neighbors. The symbol $\Box^{(2)}_{ab}=1$ if $\left|a\right\rangle$ and $\left|b\right\rangle$ differ by either a single six-site loop, or two distinct four-site loops in their transition graph, and is zero otherwise.

We can now construct orthogonal “dimer states” via

[TABLE]

The expansion in Eq. (29) formally leads to

[TABLE]

Figure 8 shows this construction for a small system of six sites. Although this expansion is commonly used in the literature Rokhsar and Kivelson (1988), it is at best formal. Intuitively one expects that for a given $f$, any given orthogonal dimer state $|\bar{a}\rangle$ will differ from the singlet covering $|a\rangle$ by a finite density of loops. This intuition is reflected in the fact that subsequent terms in Eq. (31) contain ever higher factors of the volume of space, and the limit $f\to\infty$ does not commute with the thermodynamic limit. For example, the diagonal element of third term in parentheses, $\sum_{c}\Box_{ac}\Box_{ca}$, is proportional to the total number of four-site loops that can be constructed, and scales as $N_{b}^{2}$, where $N_{b}$ is the number of bonds in $a$.

Despite its formal nature, we note that one can use this expansion to derive an effective dimer model Hamiltonian. We omit the details as it is a lengthy argument, and the result is the same as we found in Sec. IV.

Appendix C Beyond Singlets

In this appendix we analyze the case of Eq. (3) to the case when two of the terms are significant, namely,

[TABLE]

where as in the main text, $\hat{A}_{ij}^{FM\dagger}$ creates a pair on sites $i$ and $j$ with total spin $F$ and spin projection $M$. We will assume that $U_{0}\ll U_{2}$, so that the first term is large compared to the second, but that these two terms are large compared to all others.

The space of singlet coverings is not closed under Eq. (32), and we must enlarge our Hilbert space to include coverings with both spin-2 and spin-0 dimers—the former of which carry a quantum number $M$. For example, given two sites we have a six-dimensional Hilbert space, spanned by the singlet dimer, and the five spin-2 dimers. Given four sites, our Hilbert space is spanned by $3\times 6^{2}$ states—corresponding to the three different ways to pair up the four sites, and the six different flavors of each dimer. As in the purely singlet case, these states are not orthogonal. That is, the state |\raisebox{-4.55244pt}{\includegraphics[width=15.00002pt]{ahoriz1234.pdf}}\rangle=\hat{A}_{12}^{FM\dagger}\hat{A}_{34}^{F^{\prime}M^{\prime}\dagger}|{\rm vac}\rangle is not orthogonal to the state |\raisebox{-4.55244pt}{\includegraphics[width=15.00002pt]{avert1234.pdf}}\rangle=\hat{A}_{14}^{F^{\prime\prime}M^{\prime\prime}\dagger}\hat{A}_{23}^{F^{\prime\prime\prime}M^{\prime\prime\prime}\dagger}|{\rm vac}\rangle, regardless of angular momenta and projections, which are denoted by the different styles of lines joining the sites. Different flavor bonds on the same sites, however, are orthogonal: |\raisebox{-0.56905pt}{\includegraphics[width=15.00002pt]{a12.pdf}}\rangle=\hat{A}_{12}^{00\dagger}|{\rm vac}\rangle is orthogonal to |\raisebox{-0.56905pt}{\includegraphics[width=15.00002pt]{a12alt.pdf}}\rangle=\hat{A}_{12}^{2M\dagger}|{\rm vac}\rangle.

The low-energy space is spanned by nearest-neighbor singlets. As described in the main text, acting on states of this form with $\hat{A}_{ij}^{00\dagger}\hat{A}_{ij}^{00}$ can either move us in this space, or generate longer-range singlets. We need to calculate how operators of the form $\hat{A}_{ij}^{2M\dagger}\hat{A}_{ij}^{2M}$ act on these states. Let $a$ describe the singlet covering. If $(i,j)\in a$, then $\hat{A}_{ij}^{2M\dagger}\hat{A}_{ij}^{2M}|a\rangle=0$. If $(i,j)\not\in a$ then the action of $\hat{A}_{ij}^{2M\dagger}\hat{A}_{ij}^{2M}$ will involve the sites $i,j$ and their partners $k,l$. No other bonds matter, so we consider the action on |\raisebox{-4.55244pt}{\includegraphics[width=15.00002pt]{aijklstraight.pdf}}\rangle=\hat{A}_{ik}^{00\dagger}\hat{A}_{jl}^{00\dagger}|{\rm vac}\rangle. The notation does not imply any spatial relationship between the sites—just that they are connected. We then calculate |\psi_{M}\rangle=\hat{A}_{ij}^{2M\dagger}\hat{A}_{ij}^{2M}|\raisebox{-4.55244pt}{\includegraphics[width=15.00002pt]{aijklstraight.pdf}}\rangle as

[TABLE]

where we have used $C^{00}_{n,-n}=(2f+1)^{-1/2}(-1)^{f-n}$. The Hamiltonian will always generate a superposition of all different $M$, $|\psi\rangle=\sum_{M}|\psi_{M}\rangle.$ Returning to the space containing only singlets requires acting with either $\hat{A}_{ik}^{00\dagger}\hat{A}_{ik}^{00}$ or $\hat{A}_{jl}^{00\dagger}\hat{A}_{jl}^{00}$. We therefore calculate $|\phi\rangle=$ $\hat{A}_{ik}^{00\dagger}\hat{A}_{ik}^{00}|\psi\rangle$ $=\hat{A}_{jl}^{00\dagger}\hat{A}_{jl}^{00}|\psi\rangle$ as

[TABLE]

Below we show that $\sum_{M}\left(C_{m,M-m}^{2M}\right)^{2}=5/(2f+1)$, which then gives

[TABLE]

This last Clebsch-Gordan identity is a special case of the more general result $V_{Ffm}=\sum_{M}(C^{FM}_{m,M-M})^{2}=(2F+1)/(2f+1)$, which is proven by writing

[TABLE]

where $\hat{X}={\rm Tr}_{2}\sum_{M}|F,M\rangle\langle F,M|$ is the trace over the second spin of the projector into the space of fixed $F$. This operator clearly transforms as a singlet under rotation, and hence $V_{Ffm}$ must be independent of $m$. We therefore sum over $m$, and divide by $2f+1$ to find

[TABLE]

Having established the action of the various terms in the Hamiltonian on dimer coverings, we perturbatively eliminate the coverings, which either contain longer range singlet bonds, or any $F=2$ bonds. This process gives us an effective model which only involves nearest-neighbor singlets. The contributions from long-range singlets are identical to those derived in Sec. IV. Below we show that the leading contributions from the $F=2$ bonds renormalize $V$, while leaving $t$ and $t^{\prime}$ unchanged.

Let $a$ be a nearest-neighbor singlet covering, and consider nearest neighbors $i$ and $j$ such that $a$ does not contain the bond connecting them: i.e., $(i,j)\not\in a$. We will separately consider the case where the partners of $i$ and $j$ are also nearest neighbors, and the case where they are not. The first circumstance corresponds to parallel bonds. In that case, acting with $\hat{A}_{ij}^{2M\dagger}\hat{A}_{ij}^{2M\dagger}$ yields a state with two-fewer nearest-neighbor singlet bonds, but two extra nearest-neighbor $F=2$ bonds, and hence an excited state with energy $4J^{2}/U_{0}-4J^{2}/U_{2}$. In the second case, one instead finds an intermediate state with two-fewer nearest-neighbor singlets, but only one extra nearest-neighbor $f=2$ bonds. The other $F=2$ bond is long-ranged. Thus the second order process in which one returns to the initial state will have different coefficients for parallel and nonparallel bonds, hence shifting $V$. There will also be an unimportant constant energy shift to all states. In particular, the change in $V$ will be

[TABLE]

where $\lambda^{\prime}=(2J^{2}/U_{2})(2f+1)^{-1}$ is the forward amplitude, $\bar{\lambda}=5(2J^{2}/U_{0})(2f+1)^{-2}$ is the backward amplitude. The factors of $2$ account for the multiplicity of processes: There are two ways to produce a given intermediate state, and two ways back. The energy denominators $\epsilon=4J^{2}/U_{0}-4J^{2}/U_{2}$ and $\epsilon^{\prime}=(4J^{2}/U_{0}-2J^{2}/U_{2})$ are the energy denominators associated with parallel and nonparallel bonds. We note that this shift is of $\mathcal{O}(f^{-3})$ and so the effect of sub-dominant scattering channels remain negligible in the large $f$ limit, furthering the validity of our description in terms of quantum dimers.