Bottom-up configuration-interaction emulations of ultracold fermions in entangled optical plaquettes: building blocks of unconventional superconductivity

TL;DR

This paper introduces a configuration-interaction method for emulating unconventional superconductivity using ultracold fermions in optical traps, analyzing entanglement, symmetry-breaking, and coupling effects in plaquette systems.

Contribution

The paper presents a novel microscopic CI methodology for simulating superconductivity phenomena in ultracold atoms, connecting to Hubbard models and exploring multi-plaquette interactions.

Findings

Single-occupancy RVB states only partially explain many-body solutions.

CI results map onto a Hubbard Hamiltonian, not a t-J model.

Symmetry-breaking and coupling effects significantly influence spectral properties.

Abstract

A microscopic configuration-interaction (CI) methodology is introduced to enable bottom-up Schroedinger-equation emulation of unconventional superconductivity in ultracold optical traps. We illustrate the method by exploring the properties of Lithium-6 atoms in a single square plaquette in the hole-pairing regime, and by analyzing the entanglement (symmetry-preserving) and disentanglement physics (via symmetry-breaking, associated with the separation of charge and spin density waves) of two coupled plaquettes in the same regime. The single-occupancy RVB states contribute only partially to the exact many-body solutions, and the CI results map onto a Hubbard Hamiltonian, but not onto the double-occupancy-excluding t-J one. For the double-plaquette case, effects brought about by breaking the symmetry between two weakly-interacting plaquettes, either by distorting, or by tilting and detuning, one of the plaquettes with respect to the other, as well as spectral changes caused by increased coupling between the two plaquettes, are explored.

Figures

Equations

$$V({\mathbf{r}}_i,{\mathbf{r}}_j)=\frac{\lambda }{{\sigma }^2\pi }{e}^{-({\mathbf{r}}_i-{\mathbf{r}}_j{)}^2/{\sigma }^2}.$$

$${\Delta }_p=2E_{\mathrm{GS}}(N=3)-[E_{\mathrm{GS}}(N=4)+E_{\mathrm{GS}}(N=2)],$$

$${\mathcal{H}}_{\mathrm{MB}}=\sum _{i=1}^{N}H(i)+\sum _{i=1}^N\sum _{j>i}^{N}V({\mathbf{r}}_i,{\mathbf{r}}_j),$$

$${\chi }_j(x,y)={\phi }_j(x,y)\alpha ,\text{mbox}if1\le j\le K,$$

$${\chi }_j(x,y)={\phi }_{j-K}(x,y)\beta ,\text{mbox}ifK<j\le 2K,$$

$${\Phi }_{N,q}^{\text{CI}}({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)=\sum _I{C}_I^q{\Psi }_I^N({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N),$$

$$\Psi^{N}_{I}=\frac{1}{\sqrt{N!}}\left|\begin{array}[]{ccc}\chi_{j_{1}}({\bf r}_{1})&\dots&\chi_{j_{N}}({\bf r}_{1})\\ \vdots&\ddots&\vdots\\ \chi_{j_{1}}({\bf r}_{N})&\dots&\chi_{j_{N}}({\bf r}_{N})\\ \end{array}\right|,$$

$${\mathcal{H}}_{\mathrm{MB}}{\Phi }_{N,q}^{\text{EXD}}=E_{N,q}^{\text{EXD}}{\Phi }_{N,q}^{\text{EXD}}$$

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }d{\mathbf{r}}_{1}d{\mathbf{r}}_2{\phi }_i^{\ast }({\mathbf{r}}_1){\phi }_j^{\ast }({\mathbf{r}}_2)V({\mathbf{r}}_1,{\mathbf{r}}_2){\phi }_k({\mathbf{r}}_1){\phi }_l({\mathbf{r}}_2),$$

$${\hat{\mathbf{S}}}^2{\Psi }_I^N=\left[(N_{\alpha }-N_{\beta }{)}^2/4+N/2+\sum _{i<j}{\varpi }_{ij}\right]{\Psi }_I^N,$$

$${\hat{\mathcal{P}}}_{xy}{\Phi }_{N,q}^{\text{CI}}({\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N)={\Phi }_{N,q}^{\text{CI}}(-{\mathbf{r}}_1,-{\mathbf{r}}_2,\dots ,-{\mathbf{r}}_N)$$

$$n(\mathbf{r})=⟨{\Phi }^{\text{CI}}|\sum _{i=1}^N\delta ({\mathbf{r}}_i-\mathbf{r})|{\Phi }^{\text{CI}}⟩.$$

$$\text{SPOL}(\mathbf{r})=⟨{\Phi }^{\text{CI}}|\sum _{i=1}^N\delta ({\mathbf{r}}_i-\mathbf{r})({\delta }_{\uparrow {\sigma }_i}-{\delta }_{\downarrow {\sigma }_i})|{\Phi }^{\text{CI}}⟩.$$

$$P(\mathbf{r}\sigma ,{\mathbf{r}}_0{\sigma }_0)=$$

$$⟨{\Phi }^{\text{CI}}|\sum _{i=1}^N\sum _{j\ne i}^N\delta ({\mathbf{r}}_i-\mathbf{r})\delta ({\mathbf{r}}_j-{\mathbf{r}}_0){\delta }_{\sigma {\sigma }_i}{\delta }_{{\sigma }_o{\sigma }_j}|{\Phi }^{\text{CI}}⟩,$$

$$P(\mathbf{r}\sigma ;{\mathbf{r}}_1$$

$$\begin{array}{rl}{\chi }_{00}=& \sqrt{\frac{1}{3}}\sin \theta |\uparrow \uparrow \downarrow \downarrow ⟩+\left(\frac{1}{2}\cos \theta -\sqrt{\frac{1}{12}}\sin \theta \right)|\uparrow \downarrow \uparrow \downarrow ⟩\\ & -\left(\frac{1}{2}\cos \theta +\sqrt{\frac{1}{12}}\sin \theta \right)|\uparrow \downarrow \downarrow \uparrow ⟩-\left(\frac{1}{2}\cos \theta +\sqrt{\frac{1}{12}}\sin \theta \right)|\downarrow \uparrow \uparrow \downarrow ⟩\\ & +\left(\frac{1}{2}\cos \theta -\sqrt{\frac{1}{12}\sin \theta }\right)|\downarrow \uparrow \downarrow \uparrow ⟩+\sqrt{\frac{1}{3}}\sin \theta |\downarrow \downarrow \uparrow \uparrow ⟩\end{array}$$

$$\sqrt{\frac{1}{3}}\sin \theta |\uparrow \uparrow \downarrow \downarrow ⟩;$$

$$\left(\frac{1}{2}\cos \theta -\sqrt{\frac{1}{12}}\sin \theta \right)|\uparrow \downarrow \uparrow \downarrow ⟩;$$

$$\left(\frac{1}{2}\cos \theta +\sqrt{\frac{1}{12}}\sin \theta \right)|\uparrow \downarrow \downarrow \uparrow ⟩.$$

$$\begin{array}{rl}& {\Pi }_{\uparrow \downarrow }(1,2)=\\ & {\left(\frac{1}{2}\cos \theta -\sqrt{\frac{1}{12}}\sin \theta \right)}^2+{\left(\frac{1}{2}\cos \theta +\sqrt{\frac{1}{12}}\sin \theta \right)}^2\end{array}$$

$$\begin{array}{rl}& \delta =({\Pi }_{\uparrow \downarrow }(1,2)-Vol(2){)}^2+\\ & ({\Pi }_{\uparrow \downarrow }(1,3)-Vol(3){)}^2+({\Pi }_{\uparrow \downarrow }(1,4)-Vol(4){)}^2,\end{array}$$

$$\begin{array}{rl}{\chi }_{d-RVB}=& \frac{1}{2\sqrt{3}}|\uparrow \uparrow \downarrow \downarrow ⟩-\frac{1}{\sqrt{3}}|\uparrow \downarrow \uparrow \downarrow ⟩+\frac{1}{2\sqrt{3}}|\uparrow \downarrow \downarrow \uparrow ⟩+\\ & \frac{1}{2\sqrt{3}}|\downarrow \uparrow \uparrow \downarrow ⟩-\frac{1}{\sqrt{3}}|\downarrow \uparrow \downarrow \uparrow ⟩+\frac{1}{2\sqrt{3}}|\downarrow \downarrow \uparrow \uparrow ⟩\end{array}$$

$$1-\left|\frac{{\theta }_{CI}-{\theta }_{d-RVB}}{{\theta }_{d-RVB}}\right|=99.98\%$$

$$H=-t\sum _{<i,j>,\sigma }(c_{i,\sigma }^{†}{c}_{j,\sigma }+c_{j,\sigma }^{†}{c}_{i,\sigma })+U\sum _{i=1}^N{n}_{i\uparrow }{n}_{i\downarrow }$$

$$\displaystyle\begin{split}\psi_{gs}&=0.455\big{(}|\uparrow,\downarrow,\uparrow,\downarrow\rangle+|\uparrow,\downarrow,\uparrow,\downarrow\rangle\big{)}\\ &+0.228\big{(}-|\uparrow,\uparrow,\downarrow,\downarrow\rangle-|\uparrow,\downarrow,\downarrow,\uparrow\rangle-|\downarrow,\uparrow,\uparrow,\downarrow\rangle-|\downarrow,\downarrow,\uparrow,\uparrow\rangle\big{)}\\ &+0.149\big{(}|\circ,\downarrow,\uparrow,\uparrow\downarrow\rangle-|\circ,\uparrow,\downarrow,\uparrow\downarrow\rangle\big{)}+0.149\big{(}|\uparrow,\downarrow,\circ,\uparrow\downarrow\rangle-|\downarrow,\uparrow,\circ,\downarrow\uparrow\rangle\big{)}\\ &+0.149\big{(}|\downarrow,\circ,\uparrow\downarrow,\uparrow\rangle-|\uparrow,\circ,\uparrow\downarrow,\downarrow\rangle\big{)}+0.149\big{(}|\uparrow,\downarrow,\uparrow\downarrow,\circ\rangle-|\downarrow,\uparrow,\uparrow\downarrow,\circ\rangle\big{)}\\ &+0.149\big{(}|\downarrow,\uparrow\downarrow,\circ,\uparrow\rangle-|\uparrow,\uparrow\downarrow,\circ,\downarrow\rangle\big{)}+0.149\big{(}|\circ,\uparrow\downarrow,\uparrow,\downarrow\rangle-|\circ,\uparrow\downarrow,\downarrow,\uparrow\rangle\big{)}\\ &+0.149\big{(}|\uparrow\downarrow,\circ,\downarrow,\uparrow\rangle-|\uparrow\downarrow,\circ,\uparrow,\downarrow\rangle\big{)}+0.149\big{(}|\uparrow\downarrow,\uparrow,\downarrow,\circ\rangle-|\uparrow\downarrow,\downarrow,\uparrow,\circ\rangle\big{)}\\ &+0.078\big{(}-|\circ,\circ,\uparrow\downarrow,\uparrow\downarrow\rangle+|\circ,\uparrow\downarrow,\uparrow\downarrow,\circ\rangle+|\uparrow\downarrow,\circ,\circ,\uparrow\downarrow\rangle-|\uparrow\downarrow,\uparrow\downarrow,\circ,\circ\rangle\big{)}\\ &\;\;\;\;\;=0.79\chi_{d-RV\!B}+0.38\chi_{do}\end{split}$$

$${\chi }_{s-RVB}=\frac{1}{2}\left(|\uparrow \uparrow \downarrow \downarrow ⟩+|\downarrow \downarrow \uparrow \uparrow ⟩-|\downarrow \uparrow \uparrow \downarrow ⟩-|\uparrow \downarrow \downarrow \uparrow ⟩\right)$$

$$1-\left|\frac{{\theta }_{CI}-{\theta }_{s-RVB}}{{\theta }_{s-RVB}}\right|=99.95\%$$

$$\displaystyle\begin{split}\psi_{2nd\;exc.}&=0.361\big{(}-|\uparrow,\uparrow,\downarrow,\downarrow\rangle+|\uparrow,\downarrow,\downarrow,\uparrow\rangle+|\downarrow,\uparrow,\uparrow,\downarrow\rangle-|\downarrow,\downarrow,\uparrow,\uparrow\rangle\big{)}\\ &+0.185\big{(}-|\uparrow\downarrow,\circ,\uparrow\downarrow,\circ\rangle-|\circ,\uparrow\downarrow,\circ,\uparrow\downarrow\rangle\big{)}\\ &+0.153\big{(}|\circ,\downarrow,\uparrow,\uparrow\downarrow\rangle-|\circ,\uparrow,\downarrow,\uparrow\downarrow\rangle\big{)}+0.153\big{(}-|\uparrow,\downarrow,\circ,\uparrow\downarrow\rangle+|\downarrow,\uparrow,\circ,\downarrow\uparrow\rangle\big{)}\\ &+0.153\big{(}|\downarrow,\circ,\uparrow\downarrow,\uparrow\rangle-|\uparrow,\circ,\uparrow\downarrow,\downarrow\rangle\big{)}+0.153\big{(}-|\uparrow,\downarrow,\uparrow\downarrow,\circ\rangle+|\downarrow,\uparrow,\uparrow\downarrow,\circ\rangle\big{)}\\ &+0.153\big{(}|\downarrow,\uparrow\downarrow,\circ,\uparrow\rangle-|\uparrow,\uparrow\downarrow,\circ,\downarrow\rangle\big{)}+0.153\big{(}-|\circ,\uparrow\downarrow,\uparrow,\downarrow\rangle+|\circ,\uparrow\downarrow,\downarrow,\uparrow\rangle\big{)}\\ &+0.153\big{(}|\uparrow\downarrow,\circ,\downarrow,\uparrow\rangle-|\uparrow\downarrow,\circ,\uparrow,\downarrow\rangle\big{)}+0.153\big{(}-|\uparrow\downarrow,\uparrow,\downarrow,\circ\rangle+|\uparrow\downarrow,\downarrow,\uparrow,\circ\rangle\big{)}\\ &+0.093\big{(}-|\circ,\circ,\uparrow\downarrow,\uparrow\downarrow\rangle-|\circ,\uparrow\downarrow,\uparrow\downarrow,\circ\rangle-|\uparrow\downarrow,\circ,\circ,\uparrow\downarrow\rangle-|\uparrow\downarrow,\uparrow\downarrow,\circ,\circ\rangle\big{)}\\ &\;\;\;\;\;=0.72\chi_{s-RV\!B}+0.48\chi_{do}\end{split}$$

$$v(x_i,y_i)=v_{dw}(x_i)+v_{dw}(y_i),$$

Full text

I Introduction

Rapid experimental advances in the creation of finite systems of ultracold atoms using few optical traps (bottom-up) joch12 ; joch15 ; joch15.2 ; kauf14 ; kauf15 or through the use of extended optical lattices (top-down) hule15 ; bloc16 ; grei16 ; zwie16 are promising approaches toward realization of Feynman’s vision of a quantum simulator feyn82 , capable of finding solutions to systems that are otherwise numerically and/or analytically intractable. However, apart from a few double-well investigations yann15 ; yann16 ; poll15 , to date model-independent microscopic studies of multiwell systems providing theoretical insights and guidance to experimental efforts are largely lacking.

Here we introduce a configuration-interaction (CI) Schrödinger-equation methodology yann07 ; yann09 ; yann15 ; yann16 ; szabobook for exploring finite plaquette systems assembled from individual optical traps; for a brief description of the CI method, see Appendix A.

These systems are fundamental stepping stones toward bottom-up realization of large scale checkerboard or square ulracold atom lattices which are promising candidates for emulating the physics underlying $d$-wave high-Tc superconductivity ande08 ; fuku08 ; scal96 ; tsai06 ; tsai08 in optical lattices hofs02 ; zoll06 ; rey09 ; pare08 ; nasc12 ; hule15 ; bloc16 ; grei16 ; zwie16 . The work described here, demonstrating the feasibility of such exact CI calculations for interatomic contact interactions, can be extended to electronic plaquettes, i.e., to quantum-dot-arrays, governed by long-range Coulomb interactions.

The plan of the paper is as follows:

In section II we explore first the properties of ultracold fermionic atoms (6Li) confined in a single square plaquette (4-sites) in the regime of hole pairing, and subsequently analyze the entanglement physics of two coupled plaquettes in the hole-pairing regime.

In section III, for the “hole-doped” coupled plaquettes (8-sites, six atoms), we analyze the wave function anatomy of the entangled yann15 ; yann16 ; poll15 (Schrödinger-cat) two almost-degenerate manifolds, comprising: (${\it A}$) the ground state (GS) and 1st excited (1EX) state, and (${\it B}$) the two higher excited states (2EX and 3EX). The almost-degenerate states have good, but opposite, parities. When symmetry-broken (SB) – either by superposing $(\pm)$ the degenerate pair in each manifold, or via offsetting the energies (tilting) of the two plaquettes – the SB states coming from ${\it A}$ are characterized by a particle (“charge”)-density modulation (“wave”), i.e, a CDW portraying the hole-paired, (4,2) or (2,4), components, whereas the SB states originating from ${\it B}$ remain in a (3,3) particle distribution (each of the plaquettes having an unpaired hole), but exhibit a spin-polarizarion density modulation (“wave”), i.e., a SDW. For the the double-plaquette case, we further explore effects brought about by breaking the symmetry between two weakly-coupled plaquettes, either by distorting one of the plaquettes, or by tilting and detuning one of the plaquettes with respect to the other. Spectral changes caused by increased tunnel coupling (e.g., by decreased inter-plaquette distance) between the two plaquettes are also considered.

Noteworthy is our finding that the GS $d$-wave resonating valence bond (RVB) state contributes only partially to the exact many-body wave function – i.e., double occupancies (referred also as doublons) need be included. Indeed, our microscopic results map properly onto a Hubbard Hamiltonian (including extended Hubbard models fali69 ; lueh15 depending upon the range of the experimental parameters), but not onto the double-occupancy-excluding $t$-$J$ model ande08 ; rice88 ; fuku08 . Our conclusions agree with recent bloc16 ; grei16 ; zwie16 observations of doublons in two-dimensional (2D) optical lattices.

The rest of the paper comprises the Summary (Sec. IV) and three Appendices, concerning a brief description of the CI method (Appendix A), the mathematical definitions of single-particle densities and two-body and $N$-body conditional probability distributions (Appendix B), and the specification of the correspondence between spin and RVB functions and CI many-body wave functions (Appendix C).

II The single plaquette

The short-range 2-body repulsion in the Hamiltonian is described by a Gaussian of width $\sigma$, i.e., by

[TABLE]

$\lambda$ is the strength parameter. Here and throughout the paper: $\sigma=\sqrt{2}l_{0}/10=0.1833$ $\mu$m, where the oscillator length $l_{0}^{2}=\hbar/(M_{{}^{6}{\rm Li}}\omega$), $M_{{}^{6}{\rm Li}}$ being the 6Li mass and $\hbar\omega=1$ kHz being the trapping frequency of the plaquette potential wells note1 . This form of interaction provides a good approximation for the atom-atom interactions and avoids the peculiarities of the delta function in two dimensions doga13 .

The potential surface of the 4-site (or 8-site) plaquette is constructed with the help of the two-center/smooth-neck oscillator Hamiltonian which was previously introduced in Refs. yann09 ; yann15 ; yann16 ; yann99 . The smooth-neck interwell barrier $V_{b}$ can be varied independently and is controlled by the parameter $\epsilon_{b}=V_{b}/V_{0}$, where $V_{0}$ is the intersection height of the two bare potential parabolas from neighboring sites, i.e., prior to inserting the smooth-neck potential contribution yann15 ; yann09 ; yann99 . For a graphical illustration of the two-center-oscillator (TCO) double-well potential, including the definitions of $V_{b}$ and $V_{0}$, see Fig. 1. Here and throughout the paper: the intersite distance in a single plaquette is $d_{w}=6$ $\mu$m and $\epsilon_{b}=0.5$ (yielding $V_{b}=1.34$ kHz), unless noted otherwise. The potential surfaces of the 4-site and 8-site plaquettes is constructed by combining such TCO potentials along the $x$ and $y$ directions (see Supplemental Material for details suppl ).

The CI single-plaquette low-energy spectrum as a function of the two-body repulsion strength $\lambda$, is displayed in Fig. 2(a), with the single-particle density (SPD) for the GS (see red star), or for the second excited state (blue star), being shown in Fig. 2(b); for the definition of the SPD in the CI approach, see Appendix B. A total-spin projection $S_{z}=0$ has been assumed. The CI wave functions depend both on the atom (continuous) position, ${\bf r}_{i}$, and the spin, $\sigma_{i}$, $i=1,2,\ldots$, variables. However, for plaquettes with sufficiently high interwell barriers, the microscopic structure of a CI state $|\Phi\rangle$ can be mapped onto the usual Hubbard Hilbert space involving the superposition of many primitive basis functions $\Omega_{i}\equiv|\ldots d\ldots\uparrow\ldots\downarrow\dots 0\ldots\rangle$; $d$ here denotes a doubly occupied site, $\uparrow$ and $\downarrow$ denote the site’s spin occupancy, and 0 denotes an empty site. Inset (I) in Fig. 2 describes the 4-site external potential. The coefficients $c_{i}$ of the primitives $\Omega_{i}$ in the expansion of $|\Phi\rangle$, can be extracted from the CI wave functions with the help of the spin-resolved conditional probability distributions (SR-CPDs) yann15 ; yann16 ; yann07.2 . Two SR-CPDs are used here: one (SR-CPD2), expressed as an expectation value over the many-body wavefunction, describes the space and spin correlation of a particle pair, for a given location and spin of one of them (referred to as an observation point) and the other being any other particle of the $N-1$ remaining ones (with a specified spin). The other probability distribution (SR-CPDN), given by the modulus square of the many-body function, expresses the spatial probability of finding the $N$-th particle with a specified spin when one fixes the positions and spins of the other $N-1$ particles (for the mathematical definition of SR-CPDs, see Appendix B).

Figs. 2(c,d) display the SR-CPD2s of a 4-fermion plaquette for the CI states (GS and 2EX) associated with the two stars in Fig. 2(a). The GS SR-CPD2 in Fig. 2(c) reveals the presence (89.4%) of a $d$-wave RVB component, and the 2EX SR-CPD2 in Fig. 2(d) reveals the presence (84.2%) of an $s$-wave RVB component (see the analysis in Appendix C); these percentages are derived from the double-occupancy fraction at the fixed-point site [i.e., the hump volume above the arrows in Figs. 2(c,d)].

The above SR-CPD2 analysis points to a deficiency of the $t$-$J$ model because it excludes triv04 double occupancies. We further illustrate the limitations of the $t$-$J$ model for the full range of the repulsive interaction by plotting for all three approaches (CI, Hubbard, $t$-$J$) the 4-fermion and 2-fermion total energies ($E$) for a single plaquette, as well as the double occupancy ($doc$), as a function of $U/t$; see Figs. 2(e) and 2(f), respectively. Evidently, the CI and Hubbard energies, apart from a constant shift, are in very good overall agreement, whereas the $t$-$J$ values deviate greatly. The exact CI results provide a here-to-date lacking benchmark for assessments of the validity of the $t$-$J$ model and its variants ande08 ; rice88 ; fuku08 ; spal17 , as well as the Hubbard-calculated double occupancies.

The hole-pairing gap may be defined in terms of the energies of a 4-site plaquette, as follows tsai06 :

[TABLE]

where $N$ is the total number of fermions on the plaquette.

CI-calculated pairing gaps are shown in Fig. 3 for two interwell distances: $d=6$ $\mu$m [in (a,b)], and $d=2.5$ $\mu$m [in (c,d)]. Fig. 3(a) displays $\Delta_{p}$ as a function of the repulsion sterngth, $\lambda$, and the interwell parameter $\epsilon_{b}$ in a single plaquette. The gap maximum ($\sim$ 0.72 Hz) occurs at ($\lambda\sim 0.30$ $l_{0}^{2}\hbar\omega$, $\epsilon_{b}\sim 0.43$).

To compare with the results of the Hubbard model, Fig. 3(b) displays the same CI results for $\Delta_{p}$, but with all energies expressed in units of the intersite tunneling parameter $t$ and the interaction strength $\lambda$ expressed as the ratio of $U/t$, with $U$ being the Hubbard on-site repulsion extracted from the CI calculation. The maximum of $\Delta_{p}$ is now $0.045t$, occurring at ($U/t=2.6$, $\epsilon_{b}\geq 0.7$).

When $\epsilon_{b}>0.4$ yann15 ; yann09 ; yann99 , the maximum value (i.e., $0.045t$) of $\Delta_{p}$ and the range $0\leq U/t<4.8$, corresponding to $\Delta_{p}>0$ (i.e., hole pairing), agree with those found for a single-plaquette pure Hubbard model tsai06 ; tsai08 ; schum02 . The additional dependence of $\Delta_{p}$ on $\epsilon_{b}$ [see Fig. 3(b)] cannot be described by the standard Hubbard model; it reflects the effect of Hamiltonian terms that are present in the CI calculation, but are absent in the standard Hubbard model, pointing to possible applications of the extended Hubbard Hamiltonian fali69 ; lueh15 in the optical traps assemblage. Particularly relevant here is the off-site repulsion $V$ which effectively reduces fali69 the on-site $U$. However, $V$ decreases strongly for increasing intersite barrier heights $\epsilon_{b}$ (see Fig. 4), whereas $U$ is highly insensitive. When $V$ becomes sufficiently small relative to $U$, the standard Hubbard single-plaquette results are recovered [see inset in Fig. 3(b)].

$U$ in Fig. 4 is the on-site repulsion as determined from the ground state energy evolution as a function of the interaction strength $\lambda$ in a single well. $V$ is the inter-site repulsive interaction, determined through the matrix element of $V({\bf r}_{i},{\bf r}_{j})$ (see Eq. 1) taken for the ground state wavefunctions, described by Gaussians determined through fitting to the CI SPDs. $t$ is the inter-well tunneling parameter, as determined through the energy difference between the two lowest energy states in the $y$ direction (the tunneling split). Note the insensitivity of $U$ vs $\epsilon_{b}$ compared to the strong variation in $V$. The tunneling parameter $t$ exhibits a reduced sensitivity compared to $V$. The $V$ term constitutes a perturbation to the pure Hubbard model. When the $V$ term becomes of the order of magnitude of $\Delta_{p}$, i.e., $0.05t$, the pairing gap vanishes; this happens for values of $\epsilon_{b}\lesssim 0.3$ for the plaquette considered in this figure.

The role played by the off-site interaction $V$ is further illustrated in Figs. 3(c,d), where $\Delta_{p}$ is plotted for the case when the intersite distance in the single plaquette is $d_{w}=2.5$ nm [compared with $d_{w}=6$ nm in Fig. 3(a,b)]. Clearly much higher values of the intersite potential barriers (i.e., $\epsilon_{b}>6$) are required to reach the $\Delta_{p}>0$ region.

III Weakly-coupled two-plaquette molecule

The pairing-gap behavior estimated from a single 4-well optical plaquette should be reflected in the properties of a two-plaquette molecule (TPM) when the two plaquettes are weakly coupled (wc). For the single plaquette [Eq. (2) and Fig. 3, $\Delta_{p}$ indicates that the GS of the wc-TPM is associated with a (left,right) (4,2) or a (2,4) particle distribution, with the equal-particle arrangement between the (left, right) plaquettes, i.e., (3,3), corresponding to an excited state; $\Delta_{p}<0$ indicates reverse energy ordering. The GS and lowest excited state wave functions of the wc-TPM show complex behavior corresponding to entangled two mirror-reflected charge-density-wave (CDW) or spin-density-wave (SDW) symmetry-broken configurations. Experimental probing and quantitative analysis of such entangled states has recently been demonstrated joch15.2 , based on inducing particle escape from the optical wells by lowering one side of the trapping potential.

Fig. 5 displays CI results for the spectra of $N=6$ 6Li atoms in two weakly-coupled plaquettes separated by $D=18$ $\mu$m. The well-separation within the identical plaquettes is $d_{w}=6$ $\mu$m [Fig. 5(I)], so that the intraplaquette barriers are sufficiently high to yield a vanishing off-site $V$, see above. The corresponding low-energy spectrum is plotted versus $-1/\lambda$; the entire spectrum is plotted in Fig. 5(a) and an enlarged view of the spectrum in the interval $-9.5<1/\lambda<-7$ is displayed in Fig. 5(b). We focus on the two (lowest) pairs of degenerate states denoted as (${\it A}1$, ${\it A}2$) and (${\it B}1$, ${\it B}2$) in Fig. 5(b).

The CI method preserves all the quantum numbers, which are explicitly denoted in Fig. 5; i.e., ($\pm 1$) associated with the parities $P_{x}$ and $P_{y}$ along the $x$ and $y$ directions, and the total spin values $S=0$ (for ${\it A}1$, ${\it A}2$, ${\it B}2$) or $S=1$ (for ${\it B}1$).

The pairs ${\it A}$ and ${\it B}$ constitute the four lowest-in-energy states in the whole range $-10<-1/\lambda<-1$ that we investigated [Fig. 5(b)]. The crossing at $-1/\lambda_{c}=-2.8$ is reflected in the modified spectrum generated by subtracting the energy $E({\it B}1)$ at each point $-1/\lambda$; see the pocket in Fig. 6(a), where the horizontal axis is expressed now in units of $U/t$ and the crossing occurs at $U/t=4.8$. For $0<U/t<4.8$, the ${\it A}$ pair is lower in energy, with the $P_{x}=1$, $P_{y}=1$, and $S=0$ state (denoted as ${\it A}1$) being the ground state; for $4.8<U/t$, the ${\it B}$ pair is lower in energy, with the $P_{x}=1$, $P_{y}=1$, and $S=1$ state (denoted as ${\it B}1$) being the ground state.

To investigate the microscopic structure of the states in the ${\it A}$ and ${\it B}$ pairs, we first display the CI SPDs [$\rho({\bf r},\uparrow)+\rho({\bf r},\downarrow)$] and spin-polarization densities [SPOLs, $\rho({\bf r},\uparrow)-\rho({\bf r},\downarrow)$] for the states ${\it A}1$ and ${\it B}1$ in Fig. 6(b); see definitions in Appendix B. These CI SPDs and SPOLs are identical for both states, since both parities $P_{x}$ and $P_{y}$, as well as the total spin $S$, are preserved. In particular, the SPDs exhibit eight humps of equal height, integrating to a total number of $N=6$ particles. Furthermore, the SPOLs are structureless plane surfaces. The SPDs and SPOLs for the CI states ${\it A}2$ and ${\it B}2$ exhibit a similar behavior.

Direct insight into the nature of the GS structure is afforded by the SR-CPDN shown in Fig. 6(c), where the appearance of a peak in the left plaquette (see red peak), integrating to very-close to unit probability for finding an up-spin particle there when fixing the positions and spins of the other five atoms (see black arrows), indicates a dominant (4,2) wc-TPM configuration. Little double occupancy is portrayed by small red peaks on diagonal sites (left plaquette).

Although the CI-calculated densities appear uniform across the two plaquettes [see Fig. 6(b), the microscopic structures of the ${\it A}$ and ${\it B}$ states are different: they involve different superpositions of many primitive basis functions, $\Omega_{i}$, as exemplified by the $|\Phi\rangle$ superpositions in Figs. 6(I) and 6(II). These superpositions contain also primitives with doubly-occupied sites (denoted as $d$ earlier, not shown explicitly). For illustrative purposes, the pictorial representations in these panels describe visually one of these primitives associated with the (4,2) [Fig. 6(I)] and the (3,3) [Fig. 6(II)] configurations.

As aforementioned, detailed information can be extracted from the CI wave functions with the help of the SR-CPDs yann15 ; yann16 . However, for $N=6$ fermions in $M=8$ sites, calculating the large number of needed CPDs is time consuming. Nevertheless, key features of the hidden anatomy of the entangled states in the ${\it A}$ and ${\it B}$ pairs can be revealed through the implementation of a forced breaking of their $P_{x}$ parity symmetry, i.e., by subtracting and adding their member states and constructing the four auxiliary states ${\it A}^{SB}_{\mp}={\it A}1\mp{\it A}2$ and ${\it B}^{SB}_{\mp}={\it B}1\mp{\it B}2$.

The expectation value of $P_{x}$ in the ${\it A}^{SB}_{\mp}$ and ${\it B}^{SB}_{\mp}$ states vanishes; this symmetry-breaking results in different characteristics of the above two symmetry-broken states. The SPD for ${\it A}^{SB}_{-}$ exhibits a (4,2) particle configuration [Fig. 6(d), left], while the associated SPOL is featureless [Fig. 6(d), right]. The SPD for ${\it A}^{SB}_{+}$ corresponds to the mirror-reflected (2,4) configuration [Fig. 6(e), left], with the associated SPOL remaining structureless [Fig. 6(e), right]. The two fermions (or two holes) in the two-particle plaquette are not localized on specific sites. Indeed after the symmetry breaking, the wave functions are still an entangled superposition of many $\Omega_{i}$’s, thus distributing the particles (or holes) with equal probability over all sites in each plaquette. In a plaquette with two fermions, the volumes under the four density humps integrate to $N=2$, whereas the volumes under the four density humps in a 4-fermion plaquette integrate to $N=4$. Naturally, the two broken-symmetry states ${\it A}^{SB}_{\mp}$ can be characterized as a charge density wave (CDW), because they exhibit a modulation in the SPDs, but none in the spin polarizations.

The SPDs and spin polarizations of the ${\it B}^{SB}_{\mp}$ states are displayed in Figs. 6(f) and 6(g). In both cases, the SPDs are symmetric with respect to the left and right plaquettes. Furthermore, the volumes under the humps in each plaquette integrate to $N=3$, indicating that these states have a (3,3) configuration. However, the left-right asymmetry (due to the broken $P_{x}$ symmetry) emerges now as an asymmetry in the spin polarization. The three fermions (with total spin $S=1/2$) in one plaquette have a spin peojection $S_{z}=\pm 1/2$, while the remaining three fermions in the other plaquette have the opposite spin projection $S_{z}=\mp 1/2$. The ${\it B}^{SB}_{\mp}$ broken-symmetry states exhibit spin-density wave (SDW) characteristics, see right panels in Figs. 6(f,g).

In addition to the above analysis, the (4,2) or (2,4) components of the ${\it A}$ states, but not the (3$\uparrow$,3$\downarrow$) or (3$\downarrow$,3$\uparrow$) components of the ${\it B}$ states can be separated in actual experiments by lifting the left/right degeneracy with the help of two processes: (1) Tilting one plaquette with respect to the other [Fig. 7(a-d)], and (2) distorting one of the two plaquettes Fig. 7(e-h)]. From the CI spectrum depicted in Fig. 7(b), the degeneracy of the states in the A pair is lifted for a tilt $\Delta=0.1t$, with one of the states becoming the GS at $U/t=6$ (see vertical bar). Moreover the corresponding CI SPD plotted in Fig. 7(d) shows that the tilt induced the symmetry breaking discussed earlier, i.e., a (4.2) configuration; for smaller tilt results, see Figs. 7(a,c). Similar symmetry-breaking results (i.e., spectral energy gap corresponding to the emergence of a (3,3), CDW configuration) are found also for distortions of the plaquette-landscape [see Figs. 7(e-h)].

While we focused here on conditions allowing exploration of pure interplaquette entanglement effects, the sensitivity of the energy spectrum and entanglement characteristics to the interplaquette distance (i.e., increasing tunneling between the plaquettes) is illustrated in Fig. 8(a-c), where the spectrum for a double-plaqutee with $d_{w}=6$ $\mu$m and $D=8$ $\mu$m is displayed. The spectrum in Fig. 8(a) shows formation of an energy gap between the two lowest states [compare with Fig. 6(a)]. The SR-CPDN in Fig. 8(b), corresponding to the GS at $U/t=3.5$ (see red star), reflects contributions from both the (4,2) and (3,3) configurations (see the large red peak in the left plaquette and a smaller one in the right plaquette, corresponding to the probabilities of finding spin-up atoms in these locations). On the other hand the result [Fig. 8(c)] for the 1EX state [gray star in Fig. 8(a)] portrays formation of a pure (4,2)$-$(2,4) configuration [compare to the (4,2)$+$(2,4) GS of the wc-TPM in Fig. 6(c); the small difference between the $+$ and $-$ cases is not visible in the plots].

IV Summary

In this work we have developed and implemented a configuration-interaction-based computational methodology for obtaining exact solutions to the microscopic many-body quantum Hamiltonian describing ulracold fermionic atoms (here 6Li atoms) moving under the influence of an optically-induced multi-well confining potential surface and with short-range interatomic repulsive interactions; we remark here that the same CI methodology can be extended in a straightforward manner to treat interatomic dipolar interactions, or confined electrons interacting via long-range Coulomb interactions. A similar type of potential energy surface, made of the assembly of plaquettes – each comprising four sites (potential wells) arranged in a square geometry – has been introduced in early investigations scal96 ; tsai06 of electronic high-temperature superconductivity in cuprate materials, in the context of the two-dimensional Hubbard model often mentioned as a starting point for formulating a theoretical understanding of unconventional superconductivity lee06 ; lee08 .

The Hubbard model has been originally derived for the description of the behavior of strongly correlated electrons in solids hubb78 . However, variants of this model have also been implemented, for over a decade, targeting investigations of strongly-interacting ultracold atoms (bosonic or fermionic) in optical lattices, with the interest in such systems originating from remarkable advances in atom cooling and optical trapping techniques, that opened heretofore untapped prospects of preparing, emulating, and measuring the behavior of strongly-interacting quantum many-body systems under pristine, defect-free, environments. Analysis of such ultracold-atom emulations of interacting quantum many-body systems revealed on various occasions lueh15 that many Hubbard models that were simulated with ultracold atoms were found not to have a standard form – meaning that the corresponding Hamiltonians, required for analysis of these results, frequently had to incorporate terms that are not included in the standard Hubbard model and its customary variants. These extra terms may include correlated and occupation-dependent tunneling within the lowest band, as well as correlated tunneling and occupation of higher bands lueh15 .

The above findings connote that progress in developing future ultracold-atom emulations of interacting quantum many-body systems (including quantum magnetism and the origins of unconventional pairing mechanisms), and theoretical methods for the analysis of such emulations (including the development and implementation of effective models, such as extended and non-standard Hubbard Hamiltonians), require benchmark exact calculations (such as the ones introduced in this work), which owing to their ab-initio nature incorporate (with no restriction or approximation) all contributions to the many-body microscopic Hamiltonian. With such computational methodology, one can then proceed to construct, assess, and improve effective models for a reliable analysis of the new forthcoming data. Indeed, in this work we focused on providing benchmark exact quantum-mechanical solutions, aiding and enabling a bottom-up approach aimed at ultracold atoms experiments, and their analysis, starting from a single ultracold atom plaquette as a building block, and progressing in a systematic manner to double plaquettes (with variable inter-plaquette couplings) and larger, more complex, systems.

Our exact Schrödinger-equation ultracold-fermionic-atoms-plaquette simulations demonstrated entangled $d$- and $s$- wave RVB states coexisting with partial double occupancies, uncovered hole-pairing phase diagrams, and explored the robustness of the energy spectrum and entanglement of a double-plaquette system (assembled from entangled-multipartite-plaquettes) through symmetry-breaking via interplaquette tilting and distortion and as a function of the strength of interplaquette coupling. These calculations may serve as a launchpad for further exact calculations, the development of approximate treatments baer09 ; yoko06 ; mezz16 , and experiments (including the use of site resolved microscopy enabling direct observation of charge and spin correlations bloc16 ; grei16 ; zwie16 ) on systems with hierarchically increasing complexity rega16 ; barr16 ; endr16 .

Acknowledgements.

Acknowledgment: Work supported by the Air Force Office of Scientic Research under Award No. FA9550-15-1-0519. Calculations were carried out at the GATECH Center for Computational Materials Science. We thank R. Schumann for correspondence concerning the analytic Hubbard-model solutions of the 4-site plaquette.

Appendix A The configuration-interaction method

As aforementioned, we use the method of configuration interaction for determining the solution of the many-body problem specified by the $N$-fermion general Hamiltonian

[TABLE]

where ${\bf r}_{i}$, ${\bf r}_{j}$ denote the vector positions of the $i$ and $j$ fermions (e.g., 6Li atoms). This Hamiltonian is the sum of a single-particle part $H(i)$, which includes a kinetic energy term and a single-particle external confinement potential (see Fig. 1) that expresses formation of an interwell barrier between the individual wells, and the two-particle interaction $V({\bf r}_{i},{\bf r}_{j})$.

For the case of 2D ultracold atoms, the two-body repulsion is taken as a Gaussian given by Eq. (1).

In the CI method, one writes the many-body wave function $\Phi^{\text{CI}}_{N}({\bf r}_{1},{\bf r}_{2},\ldots,{\bf r}_{N})$ as a linear superposition of Slater determinants $\Psi^{N}({\bf r}_{1},{\bf r}_{2},\ldots,{\bf r}_{N})$ that span the many-body Hilbert space and are constructed out of the single-particle spin-orbitals

[TABLE]

and

[TABLE]

where $\alpha(\beta)$ denote up (down) spins. Namely

[TABLE]

where

[TABLE]

and the master index $I$ counts the number of arrangements $\{j_{1},j_{2},\ldots,j_{N}\}$ under the restriction that $1\leq j_{1}<j_{2}<\ldots<j_{N}\leq 2K$. Of course, $q=1,2,\ldots$ counts the excitation spectrum, with $q=1$ corresponding to the ground state.

The many-body Schrödinger equation

[TABLE]

transforms into a matrix diagonalizatiom problem, which yields the coefficients $C_{I}^{q}$ and the eigenenergies $E^{\text{CI}}_{N,q}$. Because the resulting matrix is sparse, we implement its numerical diagonalization employing the well known ARPACK solver arpack .

The matrix elements $\langle\Psi_{N}^{I}|{\cal H}|\Psi_{N}^{J}\rangle$ between the basis determinants [see Eq. (7)] are calculated using the Slater rules szabobook . Naturally, an important ingredient in this respect are the matrix elements of the two-body interaction,

[TABLE]

in the basis formed out of the single-particle spatial orbitals $\varphi_{i}({\bf r})$, $i=1,2,\ldots,K$. In our approach, these matrix elements are determined numerically and stored separately.

The Slater determinants $\Psi^{N}_{I}$ [see Eq. (7)] conserve the third projection $S_{z}$, but not the square $\hat{\bf S}^{2}$ of the total spin. However, because $\hat{\bf S}^{2}$ commutes with the many-body Hamiltonian, the CI solutions are automatically eigenstates of $\hat{\bf S}^{2}$ with eigenvalues $S(S+1)$. After the diagonalization, these eigenvalues are determined by applying $\hat{\bf S}^{2}$ onto $\Phi^{\text{CI}}_{N,q}$ and using the relation

[TABLE]

where the operator $\varpi_{ij}$ interchanges the spins of fermions $i$ and $j$ provided that their spins are different; $N_{\alpha}$ and $N_{\beta}$ denote the number of spin-up and spin-down fermions, respectively.

When there is no tilt [i.e., $\Delta=0$, see Fig. 1(b)], the $xy$-parity operator associated with reflections about the origin of the axes is defined as

[TABLE]

and has eigenvalues $\pm 1$. With the separable in $x$ and $y$ basis that we use, it is easy to calculate the parity eigenvalues for the Slater determinants, Eq. (7), that span the many-body Hilbert space. The many-body Hamiltonian used in this paper (without an applied magnetic field) conserves also the partial $\hat{\cal P}_{x}$ and $\hat{\cal P}_{y}$ parities.

Appendix B Single particle densities and conditional probability distributions – SPD, SPOL, SR-CPD2, and SR-CPDN

Denoting the CI wave function as $|\Phi^{\text{CI}}\rangle$, the single-particle density (SPD) is defined as

[TABLE]

Furthermore the single-particle spin polarization (SPOL) is defined as

[TABLE]

where $\uparrow$ and $\downarrow$ denote the up and down values, respectively, of the $\sigma$ spin variable. The SPOL is the difference between the up and down spin densities.

We probe the intrinsic structure of the CI eigenstates using the spin resolved conditional probability distribution $P({\bf r}\sigma,{\bf r}_{0}\sigma_{0})$ defined by the expression yann09 ; yann07.2

[TABLE]

where $\delta$ is the Dirac delta function, $N$ the total number of particles, $i,j$ are particle indices, $\sigma_{i}$ and $\sigma_{j}$ represent the spins of particles $i$ and $j$. The position of the fixed point is ${\bf r}_{0}$ and the spin of the fixed particle is $\sigma_{0}$; $r_{0}$ (together with the associated spin at that location, $\sigma_{0}$) is referred to as the observation point. $\sigma$ is the spin of the particle whose spatial distribution we want to know. Since the spin-resolved conditional probability distribution [in Eq. (14)] is represented by a two-body operator we denote it as SR-CPD2.

However, we found that sometimes SR-CPD2s alone are not sufficient to fully decipher the intrinsic configuration of the emerging quantum states. In such a case, one needs to calculate higher correlation functions. In this paper, we use the SR-CPDN ($N$-point correlation function) defined as the modulus square of the full many-body CI wave function, i.e.,

[TABLE]

where one fixes the positions and spins of $N-1$ particles and inquiries about the (conditional) probability of finding the $N$th particle with spin $\sigma$ at any position ${\bf r}$.

Appendix C Spin eigenfunction determination

C.1 RVB state in a single plaquette

In this section we provide a detailed description of the procedure used to identify the d-wave RVB state in our spin resolved SR-CPD2s.

As first described in Ref. yann09 one can map SR-CPD2s to spin eigenfunctions by analyzing the volumes underneath the SR-CPD2s. The versatility of this method has been demonstrated several times yann09 ; yann16 and detailed explanations of the methodology can be found in lidis09 ; yann16 ; yann09 . This procedure was carried out by us to identify the RVB states within our plaquettes. The starting point is the general spin eigenfunction for four spin $1/2$ fermions trapped in a single plaquette with quantum numbers $S=0,S_{z}=0$, which is given as yann09 :

[TABLE]

This function is parameterized by the angle $\theta$ and yields the quantum numbers $S=0,S_{z}=0$ for all values of $\theta$. The angle $\theta$ is a free parameter that for certain values gives spin functions with specified characteristics (e.g., the orthogonal functions mentioned below), or can be treated as a “fitting” parameter (see the discussion below Eq. (20) in connection with $\theta_{CI}$). Eq. (16) expresses in a compact form the fact that the dimension of the total-spin space for $N=4$ and $S=0$, $S_{z}=0$ is 2; see the branching diagram in Ref. paunczbook and Ref. yann16 . For example, two orthogonal basis functions can be obtained by setting $\theta=0$ and $\theta=\pi/2$. From the general spin function, one can read off the spin components contributing to a specific SR-CPD2. For instance, for a spin down SR-CPD2 with a fixed spin-up fermion on position 1 (counting from left to right in the corresponding kets), the spin components contributing to the conditional probability densities are the following three:

[TABLE]

The volume under the hump in position 2 of such a SR-CPD2 is therefore proportional to

[TABLE]

This is referred to as the partial conditional probability $\Pi_{\uparrow\downarrow}(1,2)$ (partial because it corresponds to a part of the full SR-CPD2, i.e. one peak for a specific spin configuration). Using a normalized SR-CPD2. one can directly equate the volume underneath the hump in position 2 to $\Pi_{\uparrow\downarrow}(1,2)$ and determine the angle $\theta_{CI}$. Due to the involved squares, this procedure is not necessarily unique. The unique solution can be found by comparing another hump (i.e., at position 3) to the corresponding partial conditional probability $\Pi_{\uparrow\downarrow}(1,3)$.

However it is important to note that this procedure is only exact as long as (i) The overlap between sites is sufficiently small and (ii) The amount of double occupancy is small.

From Fig. 2(c,d), it is clear that (i) is fulfilled. To assure that the possible error due to (ii) is as small as possible, we minimize

[TABLE]

where $\mathit{Vol}(i)$ represents the volume under the SR-CPD2 hump at position $i$.

C.2 $d$-wave RVB determination

For the ground state SR-CPD2 shown in Fig. 2, this procedure yielded an angle $\theta_{CI}=2.618614\pi$. It is known that the spin function for the $d$-wave RVB is given as nasc12

[TABLE]

which corresponds to an angle of $\theta_{d-RV\!B}=\frac{5\pi}{6}$. We can therefore conclude that our CI results show the presence of a $d$-wave RVB state with

[TABLE]

fidelity. However while the CI result shows the presence of an RVB with high fidelity, the SR-CPD2 also undoubtedly shows the presence of double occupancy. The volume of the hump at position 1 amounts to $10.39\%$. This might seem at first contradictory but is in fact the correct solution for an interaction strength of $U/t=2.461$. To show this we calculated the ground state solution of the Hubbard Model

[TABLE]

which is given as

[TABLE]

The first two lines of Eq. (25) are the $d$-wave RVB component and appear with dominant coefficients, however the wavefunction clearly contains contributions from doubly occupied states. Summing the squared coefficients of the spin primitives that contain doubly occupied sites yields a double occupancy of $10.17\%$ in excellent agreement with the CI result.

C.3 s-wave RVB determination

The same analysis can be performed for the $s$-wave RVB state. The spin function for an $s$-wave RVB state is given in literature asnasc12 ; rey09 :

[TABLE]

This corresponds to an angle of $\theta_{s-RV\!B}=-\frac{2\pi}{3}$. The angle determined from our CI is $\theta_{CI}=-0.66632\pi$ corresponding to a fidelity of

[TABLE]

However just like the ground state SR-CPD2, the SR-CPD2 for the second excited state shows a non-zero double occupancy. Therefore we computed the Hubbard model solution for the 2nd excited state at $U/t=2.461$:

[TABLE]

This eigenfunction predicts a double occupancy on site 1 of $14.52\%$. The volume of the hump on position 1 indicating double occupancy in our CI amounts to $15.78\%$ in good agreement with the Hubbard result.

Appendix D External Confining Potentials for Plaquettes

D.1 The external potential for a single plaquette

For our single plaquette calculations we combine the double well external potentials from li_artificial_2009 ; brandt_double-well_2015 ; yannouleas_ultracold_2016 in both the $x$ and $y$ directions.

[TABLE]

where $v_{dw}$ is the double well external potential (see Fig. 1 in the main paper). This results in a two-dimensional confining potential that has four minima and smooth connecting necks between them. We have full control over the distances between the minima in $x$ and $y$ and over the barriers between them. The set of parameters used to control the double-well potential are $V_{b}$, $\Delta$, and $d_{w}$, where $V_{b}$ is the height of the smooth inter-well barrier, $\Delta$ is the tilt between the left and right well, and $d_{w}$ is the distance between the wells li_artificial_2009 ; brandt_double-well_2015 ; yannouleas_ultracold_2016 . For a graphical illustration of the double-well potential, including the definitions of $V_{b}$, $\epsilon_{b}=V_{b}/V_{0}$, and $V_{0}$, see Fig. 1 of the main paper.

D.2 The external potential for the double (multi) plaquette

For our double-plaquette calculations we use a repetition of the single plaquette external potential from Sec. D.1. We keep the potential along $y$ as is. For the potential along $x$ we start out with three double-well potentials and use Heaviside theta functions $\Theta$ to cut the inward facing potential walls to the right and/or left of the minima. We then combine the resulting potentials. This procedure is illustrated in Fig. 9. The above procedure applies directly to large inter-plaquette distances. For smaller inter-plaquette distances end-effects of the two-plaquette assembly are alleviated through ”symmetrization” of the trapping potential at the left and right edges of the double-plaquette assembly so that it takes the shape of the intraplaquette interwell potential from the bottom of the well to the midpoint between sites.

Just like in Sec. D.1 this gives us full control over the positions of the minima in $x$ and $y$ and allows us to independently vary the barriers in the $x$ and $y$ directions. Unless otherwise stated the parameters for our double-plaquette are: $\hbar\omega=1$ kHz, $l_{0}=1.296$ $\mu$m, $d_{w}=6$ $\mu$m; (intraplaquette), $V_{b}=1.34$ kHz (intraplaquette) $(\epsilon_{b}=0.5)$, $\sigma=0.1833$ $\mu$m for all potential wells, and total-spin projection $S_{z}=0$. The $V_{b}$ between the plaquettes (interplaquette) varies depending on the interplaquette distance ($\epsilon_{b}=0.5$ is kept constant).

D.3 Solution of the single particle Hamiltonian

For rapid convergence of the full CI calculations, it is often beneficial to use a basis that already diagonalizes the single particle part of the Hamiltonian. Since the single-particle Hamiltonians used in this paper are separable in $x$ and $y$ we have to solve two one-dimensional Schrödinger equations. A software package that does so, with high performance and high accuracy, is Chebfun driscoll_chebfun_2014 . It uses Chebychev polynomials as an efficient and highly accurate approximation of the true wave functions. From the documentation:

”The implementation of Chebfun is based on the mathematical fact that smooth functions can be represented very efficiently by polynomial interpolation in Chebyshev points, or equivalently, thanks to the Fast Fourier Transform, by expansions in Chebyshev polynomials”; and ”Chebfun makes use of adaptive procedures that aim to find the right number of points automatically so as to represent each function to roughly machine precision, that is, about 15 digits of relative accuracy.”

We have confirmed the accuracy of these solutions with both symbolic and numeric calculations in MATHEMATICA and the agreement was excellent at significantly higher performance.

Appendix E Analytical and numerical solutions of the Hubbard Hamiltonian

To obtain numerical solutions to the Hubbard model we use the software package SNEG zitko_sneg_2011 , which is a ”MATHEMATICA package for symbolic calculations with second-quantization-operator expressions”. From the documentation:

”SNEG library is a MATHEMATICA package that provides a framework for performing calculations using the operators of the second quantization with an emphasis on the anti-commuting fermionic operators in the context of solid-state and atomic physics. It consists of a collection of transformation rules that define the algebra of operators and a comprehensive library of utility functions.”

For the case of a single square plaquette, the analytical results in schumann_2012 were also consulted.