Realizing the Haldane Phase with Bosons in Optical Lattices

TL;DR

This paper proposes a method to realize and detect the Haldane phase in bosonic atoms within optical lattices, highlighting the potential for observing topological order and edge states through quantum gas microscopy.

Contribution

It introduces an experimentally feasible model for the Haldane phase with bosons and details how to detect its non-local order and topological features.

Findings

Haldane phase can be realized with bosons in optical lattices

Quantum gas microscopy can detect non-local string order and edge states

Supersolid correlations are possible but challenging to observe experimentally

Abstract

We analyze an experimentally realizable model of bosons in a zig-zag optical lattice, showing that by rapidly modulating the magnetic field one can tune interaction parameters and realize an analog of the Haldane phase. We explain how quantum gas microscopy can be used to detect this phase's non-local string order and its topological edge states. We model the detection process. We also find that this model can display supersolid correlations, but argue that they only occur at parameter values which would be challenging to realize in an experiment.

Figures

Equations

$$H$$

$$-t^{\prime }\sum _{j,s}(c_{s,j+1}^{†}{c}_{s,j}+c_{s,j}^{†}{c}_{s,j+1})+\frac{U}{2}\sum _{j,s}{n}_{s,j}(n_{s,j}-1)$$

$$+U_2\sum _j[n_{1,j}+n_{1,j-1}]n_{-1,j}.$$

$$C_{ij}^{\mathrm{str}}$$

$$C_{ij}^{\mathrm{MI}}$$

$$C_{ij}^{\mathrm{DW}}$$

$$|i{⟩}_L^j=\sum _{k\alpha }{\Gamma }_{ik\alpha }^{Lj}|k{⟩}_L^{j-1}|\alpha {⟩}_a^{j-1},|i{⟩}_R^j=\sum _{k\beta }{\Gamma }_{ik\beta }^{Rj}|k{⟩}_R^{j+1}|\beta {⟩}_b^{j+1}.$$

$$\sum _{k\alpha }{\Gamma }_{ik\alpha }^{Lj}({\Gamma }_{i^{\prime }k\alpha }^{Lj}{)}^{\ast }={\delta }_{i{i}^{\prime }},\sum _i{\Gamma }_{ik\beta }^{Lj}({\Gamma }_{i{k}^{\prime }{\beta }^{\prime }}^{Lj}{)}^{\ast }={\delta }_{k{k}^{\prime }}{\delta }_{\beta {\beta }^{\prime }},$$

$$|\psi ⟩=\sum _{\ell r\alpha \beta }{\psi }_{\ell r\alpha \beta }^{(j)}|\ell {⟩}_L^j|r{⟩}_R^j|\alpha {⟩}_a^j|\beta {⟩}_b^j.$$

$${\psi }_{\mu \nu \beta t}^{(j+1)}=\sum _{\ell r\alpha }{\psi }_{\ell r\alpha \beta }^{(j)}{\Gamma }_{r\nu t}^{Rj}{\left({\Gamma }_{\mu \ell \alpha }^{L(j+1)}\right)}^{\ast }.$$

$${\rho }_{\ell {\ell }^{\prime }}^{Lj}=\sum _{r\alpha \beta }{\psi }_{\ell r\alpha \beta }^{(j)}{\left({\psi }_{{\ell }^{\prime }r\alpha \beta }^{(j)}\right)}^{\ast }.$$

$${\rho }_{\ell \alpha ,{\ell }^{\prime }{\alpha }^{\prime }}^{L(j+1)}=\sum _{r\beta }{\psi }_{\ell r\alpha \beta }^{(j)}{\left({\psi }_{{\ell }^{\prime }r{\alpha }^{\prime }\beta }^{(j)}\right)}^{\ast },$$

$${\Gamma }_{s\ell \alpha }^{L(j+1)}={\varphi }_{\ell \alpha }^s,$$

$$P_{n_1=n}=\sum _{r\alpha \beta }\sum _{\ell {\ell }^{\prime }}({\psi }_{\ell r\alpha \beta }^{(1)}{)}^{\ast }{\psi }_{{\ell }^{\prime }r\alpha \beta }^{(1)}(P_n^1{)}_{\ell {\ell }^{\prime }}.$$

$${\psi }_{\ell r\alpha \beta }^{(1)}\to \sum _{{\ell }^{\prime }}(P_{n_1}^1{)}_{\ell {\ell }^{\prime }}{\psi }_{{\ell }^{\prime }r\alpha \beta }^{(1)}.$$

$$P_{n_2=n}=\sum _{r\ell \beta }\sum _{\alpha {\alpha }^{\prime }}({\psi }_{\ell r\alpha \beta }^{(1)}{)}^{\ast }{\psi }_{\ell r{\alpha }^{\prime }\beta }^{(1)}(P_n^2{)}_{\alpha {\alpha }^{\prime }}.$$

$${\psi }_{\ell r\alpha \beta }^{(1)}\to \sum _{{\alpha }^{\prime }}(P_{n_2}^2{)}_{\alpha {\alpha }^{\prime }}{\psi }_{\ell r{\alpha }^{\prime }\beta }^{(1)}.$$

$${\psi }_{\ell r\alpha \beta }^{(2)}=\sum _{{\ell }^{\prime }{r}^{\prime }{\alpha }^{\prime }}{\psi }_{{\ell }^{\prime }{r}^{\prime }{\alpha }^{\prime }\alpha }^{(1)}{\Gamma }_{r^{\prime }r\beta }^{R1}{\left({\Gamma }_{\ell {\ell }^{\prime }{\alpha }^{\prime }}^{L2}\right)}^{\ast }.$$

$$H$$

$$+\bar{t}\sum _i({\bar{b}}_{i+1}^{†}{\bar{b}}_i+{\bar{b}}_i^{†}{\bar{b}}_{i+1})-{\bar{t}}^{\prime }\sum _i({\bar{b}}_{i+2}^{†}{\bar{b}}_i+{\bar{b}}_i^{†}{\bar{b}}_{i+2})+\frac{\bar{U}}{2}\sum _i{\bar{n}}_i({\bar{n}}_i-1)+{\bar{U}}_2\sum _i{\bar{n}}_{i+1}{\bar{n}}_i$$

$$+2\tilde{U}\sum _i{n}_i{\bar{n}}_i+{\tilde{U}}_2\sum _i(n_i{\bar{n}}_{i+1}+{\bar{n}}_i{n}_{i+1})+\frac{\tilde{U}}{2}\sum _i({\bar{b}}_i^{†}{\bar{b}}_i^{†}{b}_i{b}_i+b_i^{†}{b}_i^{†}{\bar{b}}_i{\bar{b}}_i)$$

$$+\lambda \sum _i(b_{i+1}^{†}{\bar{b}}_i+{\bar{b}}_i^{†}{b}_{i+1}-{\bar{b}}_{i+1}^{†}{b}_i-b_i^{†}{\bar{b}}_{i+1})+\delta \sum _i{\bar{n}}_i,$$

$$t=\Omega \int w_1^{\ast }(x-a/4)w_1(x+a/4)dx,\bar{t}=\Omega \int w_2^{\ast }(x-a/4)w_2(x+a/4)dx,$$

$$t^{\prime }=-\int w_1^{\ast }(x)H_0{w}_1(x+a)dx,{\bar{t}}^{\prime }=-\int w_2^{\ast }(x)H_0{w}_2(x+a)dx,$$

$$U=U_0\int |w_1(x){|}^{4}dx,\bar{U}=U_0\int |w_2(x){|}^{4}dx,$$

$$U_2=U_0\int |w_1(x){|}^2|w_1(x+a/2){|}^{2}dx,{\bar{U}}_2=U_0\int |w_2(x){|}^2|w_2(x+a/2){|}^{2}dx,$$

$$\tilde{U}=U_0\int |w_1(x){|}^2|w_2(x){|}^{2}dx,{\tilde{U}}_2=U_0\int |w_1(x){|}^2|w_2(x+a/2){|}^{2}dx,$$

$$\lambda =\Omega \int w_1^{\ast }(x-a/4)w_2(x+a/4)dx,$$

$$\delta =\int w_2^{\ast }(x)H_0{w}_2(x)dx-\int w_1^{\ast }(x)H_0{w}_1(x)dx,$$

$$C_{ij}^{\mathrm{str}}$$

$$C_{ij}^{\mathrm{MI}}$$

$$C_{ij}^{\mathrm{DW}}$$

Full text

In the past 30 years, one of the dominant themes in condensed matter theory has been the search for models where the collective excitations behave unlike any known fundamental particle. While many such fractionalized and topologically ordered models have been found topology , very few of them have been experimentally realized. Here we show how to build on a setup proposed by the NIST cold atom experimental group Spielman to explore one of the iconic fractionalized phases, the Haldane phase of a spin-1 chain Haldane .

In 1983, Haldane showed that the properties of integer and half-integer spin chains can be profoundly different Haldane ; review . Over the following decade, several researchers explored the rich properties of the integer spin chain, finding half integer spin edge modes AKLT ; Hagiwara ; Glarum , and non-local string order Nijs ; Girvin ; Kennedy . More recently, Dalla Torre, Berg, and Altman noted that similar physics should occur for spinless bosons hopping on a one-dimensional lattice: the occupation numbers on each site plays the role of the different spin states Altman ; Altman2 . Subsequently, analogs of the Haldane phase have been predicted for a number of one-dimensional Bose-Hubbard models with off-site interactions Rossini ; Batrouni ; Ejima ; Lange . One enlarges the parameter range over which the Haldane phase is stable if there is a constraint on the maximum number of particles per site. By combining a number of experimental techniques, we show how to realize a model which would be expected to support the Haldane phase. We use Density Matrix Renormalization Group (DMRG) techniques to calculate the properties of this model White ; Schollwock , and explain how to detect the exotic signatures of the Haldane phase.

In a system of one-dimensional lattice bosons, the Haldane (HI) phase lies at the intersection of the density wave (DW) phase, where double occupied sites (doublons) alternate with empty sites (holons), the Mott insulator (MI) phase, where each site is occupied by a single atom, and the superfluid (SF) phase, where the quasiparticles (doublons and holons) are free to move around. In the Haldane phase the quasiparticles are fluid but ordered: their spacing varies, but as one moves from left to right the next quasiparticle after a doublon is a holon, and vice-versa. This ground state is four-fold degenerate in a large but finite system with hard-wall boundary conditions – corresponding to the flavors of the leftmost and rightmost quasiparticles – which are bound to the edges of the system. This four-fold degeneracy was also found in the original spin context, corresponding to two spin-1/2 degrees of freedom, one sitting at each boundary.

One-dimensional bosonic system have been realized by trapping cold atoms in elongated optical traps exp1 ; exp2 ; exp3 . Anisimovas et al. showed that by using a one-dimensional (1D) spin dependent optical lattice and Raman induced hopping, one could produce the zig-zag lattice illustrated in Fig. 1(a), described by the tight-binding model Spielman

[TABLE]

Here $s=\pm 1$ labels the spin state of the atoms and $j$ is the position of the atom along the lattices. The hopping between and within these two spin states are characterized by $t$ and $t^{\prime}$. The on-site and nearest interspecies interaction are described by $U$ and $U_{2}$. Following Anisimovas et al. and shown in Fig. 1(a), one can think of Eq. (Realizing the Haldane Phase with Bosons in Optical Lattices) in one of two ways: either as a two-leg “zig-zag” lattice, or a 1D chain with next-nearest neighbor hopping. The latter is closer to the actual physical system. Following that interpretation, we introduce operators $b_{2j}=c_{-1,j}$ and $b_{2j+1}=c_{+1,j}$, in which case $t,t^{\prime}$ are nearest and next-nearest neighbor hopping parameters, while $U,U_{2}$ are on-site and nearest neighbor interaction parameters, i.e. $H=t\sum_{i}(b^{\dagger}_{i+1}b_{i}+b^{\dagger}_{i}b_{i+1})-t^{\prime}\sum_{i}(b^{\dagger}_{i+2}b_{i}+b^{\dagger}_{i}b_{i+2})+\frac{U}{2}\sum_{i}n_{i}(n_{i}-1)+U_{2}\sum_{i}n_{i+1}n_{i}$. Aside from the longer range hopping, Eq. (Realizing the Haldane Phase with Bosons in Optical Lattices) maps onto the model introduced by Dalla Torre et al. in considering polar molecules in optical lattices. Unfortunately, based on their analysis, one expects that the Haldane phase is not stable when $U/U_{2}$ is large – which is the physically relevant regime considered in Spielman . (Our numerics confirm this expectation.) Here we argue that by using a Feshbach resonance feshbach , one can reduce $U$, driving the system into the Haldane phase. The lossy nature of bosonic Feshbach resonances aids us, as the quantum Zeno effect converts the resulting large 3-body recombination rate into a suppression of the probability of having more than two particles on any given site – further stabilizing the Haldane phase.

More concretely, we consider the $F=1$, $m_{F}=1,0$ states of 87Rb. The coefficient $U$ is proportional to the scattering length associated with two atoms in the same magnetic sublevel. As illustrated in Fig. 1(b), this scattering length can be manipulated by applying a magnetic field. Near $B_{1}\sim 661.43$G there is a zero-crossing where the interactions between two $m_{F}=0$ atoms vanish, while near $B_{2}\sim 685.43$G there is a similar zero crossing for $m_{F}=1$ Marte . We envision rapidly switching the magnetic field between these two fields, as illustrated in Fig. 1(c). As long as the switching time is short compared to the other scales in the problem ($h/U,h/t\sim h/E_{R}\approx 0.27$ms for laser wavelength $\lambda=789$nm) the effective interaction in each spin channel will be given by time-averaging the instantaneous Hamiltonian $U_{\rm eff}=\int_{0}^{t}U(\tau)d\tau$ Blanes ; Bukov ; Eckardt . Even though at any given time the interactions in the two channels will be different, this time averaged interaction is the same for each spin species, and $U$ will be the same on all sites. This technique effectively halves the strength of the on-site interaction as in Fig. 1(b). The coefficient $U_{2}$ is largely unaffected. We find that one can achieve a ratio of on-site to nearest neighbor interaction of $U/U_{2}\approx 1.6$, for a lattice depth of $V_{0}=2E_{R}$ (for the effect of a higher band, see the Supplemental Material Sup ). By appropriately tuning the transverse confinement, one can take $U_{2}/t=2.5$, yielding $U/t=4$. Figure 2 shows the phase diagram for this model, revealing that these parameters place the system within the Haldane phase regime.

The zero crossings are very close to Feshbach resonances, and hence induce a large 3-body loss rate $K_{3}$. In the present circumstance this is advantageous. Following the logic in Zoller , when $K_{3}$ is large, there is a strong suppression of the process in which a third particle hops onto a site containing two other particles. This suppression can be modeled by a complex on-site three-body repulsion of strength $U_{3b}\sim-ihK_{3}n^{2}/12$. We estimate that one can get $|U_{3b}|\sim 10E_{R}$ for a typical on-site particle density $n\sim 2.1\times 10^{15}$cm*-3* and typical on-resonance three-body rate $K_{3}\sim 10^{-25}$cm6/s. Since $|U_{3b}|$ is larger than the other scales in the problem, it can be replaced by a constraint that no more than two particles can occupy any site.

With this constraint we use the DMRG to calculate the properties of the model in Eq. (Realizing the Haldane Phase with Bosons in Optical Lattices). We start with a infinite DMRG algorithm to grow the system to desired size, and then do finite DMRG sweeps until we reach convergent. This technique is typically understood as systematically optimizing a variational wavefunction in the form of a matrix product state Schollwock . The degree of approximation is controlled by the bond dimension $d$. We have considered systems as large as $L=512$ sites, and bond dimensions as large as $d=500$. The algorithm is more efficient if we alter the boundary conditions to break the potential four-fold degeneracy of the Haldane-phase groundstate, and the potential two-fold degeneracy of the density wave groundstate. In particular we used boundary conditions which pin a vacancy at the left-most site, and doublon at the right-most site. We analyze convergence with bond dimension and system size in the supplementary information. From these studies we expect that experiments on systems of size $L\sim 60$ will see significant finite-size effects near the phase boundaries, but the bulk physics is unchanged, and such experiments will be able to unambiguously observe all of the relevant physical phenomena.

The order in the Haldane, Mott insulator, density wave phases are encoded in string (str), parity (MI) and density wave (DW) correlation functions Nijs ; Altman2

[TABLE]

where $\delta n_{k}=n_{k}-1$. The phase factor $e^{i\pi\sum_{i<k<j}\delta n_{k}}=\pm 1$ depends on if the number of quasiparticles between sites $i$ and $j$ is even or odd. In the superfluid phase all of the three correlation functions fall to zero as $i$ and $j$ are separated. In the Haldane/Mott insulator phase, only $C^{\rm str}/C^{\rm MI}$ has long-range order, while in the density wave phase all the three correlation functions are nonzero.

We additionally study the single particle density matrix $C^{\rm SF}_{ij}=\langle b_{i}^{\dagger}b_{j}\rangle$ and the bipartite number fluctuations, $D_{j}=\langle N_{i<j}^{2}\rangle-\langle N_{i<j}\rangle^{2}$, where $N_{i<j}=\sum_{i=1}^{j-1}n_{i}$ is the number of particles to the left of site $j$. The superfluid phase is characterized by power-law behavior of the density matrix, and enhance bipartite number fluctuations when compared to the incompressible insulating phases. We found that these number fluctuations were the most reliable way to extract the phase boundaries between the superfluid and insulating phases Sup . In particular, due to the different scaling with system size, the number fluctuations in half the chain, $D_{L/2}$, form plateaus in each of the phases, and the phase boundaries correspond to peaks in the slope $dD_{L/2}/dt^{\prime}$. We use the full width half max of these peaks as an estimate of the accuracy of these boundaries. This approach is adapted from fluckt , and is similar to finding phase boundaries from peaks in specific heat.

Additionally, the DW to HI transition can be accurately determined from the properties of $C^{\rm DW}$. Applying finite size scaling fsc to the asymptotic behavior of this correlation function yields a DW-HI boundary which agrees with our calculation using the number fluctuations. The various superfluid-insulator transitions are not amenable to this standard finite size scaling analysis: they have behavior related to Kosterlitz-Thouless transitions, and are harder to determine. In addition to our technique of looking at number fluctuations, these transitions can be identified by looking at the excitations spectrum Rossini or superfluid stiffness Stiff , by taking moments of $C^{\rm SF}$ pai ; roomany , or by comparing the power law decay of $C^{\rm SF}$ to a Luttinger liquid model pow ; minguzzi . More discussion of the critical behavior can be found in kurdestany . Due to the significant finite size effects, it is unlikely that an experiment would be able to accurately determine these phase boundaries.

The correlation functions in Eq. (S2) are directly measurable via a quantum gas microscope Ott ; Bakr ; Bakr2 ; Sherson ; Endres . One projects the quantum state into one in which there is a definite number of particles on each site – giving a single realization of $\{n_{i}\}$. Repeating the measurement many times allows one to extract the expectation values in Eq. (S2). This technique has already been used to measure the parity order Endres .

In addition to showing the phase boundaries, Fig. 2 shows the size of correlations on a length-scale of 80 sites. In a significant part of the phase diagram, the string correlations are large but all other correlation functions vanish. This corresponds to the desired Haldane phase. The $HI$, $DW$, and $MI$ correlations at shorter lengths scales are of similar strength, but display less sharp boundaries. As would be expected, the $SF$ correlations are strongly length-dependent. They are also extremely hard to measure in an experiment. The simplest experimental knob for moving through this phase diagram is the strength of the Raman beams, which changes $t$ while leaving all other scales unchanged. Additionally, both $t$ and $t^{\prime}$ are exponentially sensitive to the lattice depth, while the ratio between $U$ and $U_{2}$ is controlled by modifying the time dependence of the magnetic field.

In addition to exploring the expectation values of the various correlation functions, we use a novel Monte-Carlo sampling algorithm to stochastically generate ‘typical’ cold-gas microscope images Sup . Given the DMRG wavefunction $|\psi\rangle$, we first calculate the probability that site-1 had $0,1$ or $2$ particles on it. We use these probabilities to choose one of these sectors, and project the wavefunction into that sector. This calculation is then repeated on site-2, using the new wavefunction… Figure 3 shows configurations generated by this algorithm, which should be representative of what is seen in an experiment. We emphasize that these are not cartoons, but rather are unbiased samples. As expected, in the HI phase the doublons and holons alternate, with a variable number of singly-occupied sites between them. This can be contrasted with the SF phase, where there is no ordering of the doublons and holons. In the DW phase, doublons and holons alternate. In these images one sees a small number of defects in the order – as should be expected. In the MI phase the images show very few holons and doublons – and those which exist are tightly bound together. In this figure we also show images with supersolid (SS) correlations that can appear when we relax the constraint forbidding double occupancy. The physics of this regime will be discussed below.

To illustrate the role of the three-body constraint, we repeated our calculations, allowing the on-site particle number to be as large as 3. Figure 4 shows the analog of Fig. 2. All correlations, except those corresponding to SF order, are much weaker. Short and medium range HI, MI, and DW correlations are detectable, but our scaling analysis suggest that for these parameters there is no long-range DW or HI order.

For small on-site interaction $U$ and finite next nearest hopping $t^{\prime}$, we find a superfluid region with short-range density wave order, which is suggestive of proximity to a supersolid phase (SS). Such a phase would be more familiar in the language of the “zig-zag” ladder picture: The atoms form a superfluid which preferentially sits on one leg of the ladder. An alternative cartoon can be constructed from the DW state “2020202020.” Because of the next-nearest neighbor hopping, one can produce a triplon-singlon pair “2020103020,” and these defects may be mobile. Forbidding triple occupancy eliminates these excitations, and prevents the occurrence of this phase in the constrained model. In addition to such triplons and singlons, the configurations in Fig. 3 display defects where atoms have hopped from even to odd sublattices. These defects are responsible for the short-range nature of the correlations.

To summarize, we have proposed a way to realize the Haldane phase in a gas of 87Rb atoms trapped in a zigzag optical lattice, where a different atomic spin state is trapped on each leg of the ladder. One reduces the on-site interactions (relative to the nearest neighbor interactions) by rapidly sweeping the magnetic field between two zero-crossings associated with Feshbach resonances in each of the spin states. The proximity to the Feshbach resonances introduces large three-body loss, which via the quantum Zeno effect prevents triple-occupation. We calculate the phase diagram of this model, and find that the Haldane phase is experimentally realizable. We modeled a quantum-gas microscope experiment, and found that one can readily identify the string order of the Haldane phase in individual images. More quantitative tests require averaging over several images. Such averaging has been used to identify other non-local order parameters Endres . We further show that without the constraint on particle number, this model shows hints of a supersolid phase (cf. Esslinger ; Ketterle ). One would need to use other techniques, however, to experimentally reach this supersolid regime. Seeing the string order in the Haldane phase would be a remarkable triumph in engineering quantum matter.

J.X. acknowledges helps from Matthew Reichl on the numerical calculation. In addition to our custom-built DMRG software, we performed many of the calculation using the ITensor library (http://itensor.org). E.J.M. acknowledges help from Joshua Squires on aspects related to the phase boundaries. This research is supported by the NBRPC (No. 2013CB922002), NNSFC (No. 11504021, 11574028), NSF (No. PHY-1508300), and FRFCU (No. FRF-TP-17-023A2).

S1 DMRG Algorithm

To study our effectively one-dimensional system we use the Density Matrix Renormalization Group (DMRG) method. There are a number of excellent reviews of the technique Schollwock and pedagogical sample code which makes it easier to understand the mechanics codeexamples . Here we give a quick introduction, and elaborate on our explanation of how we extract cold-gas microscope configurations, as in Fig. 3 of the main text.

S1.1 Configurations and blocks

The DMRG algorithm involves sequentially repartitioning the system into configurations of different shape, refining a variational wavefunction as one sweeps through the configurations. For a chain of length $L$, we consider $L-3$ different configurations. In the $j$’th configuration, the left block contains $j$ sites, the right block contains $L-j-2$ sites, and two extra sites – labeled $a$ and $b$ – lie between them. The basis for the left block is $|1\rangle_{L}^{j},|2\rangle_{L}^{j}\cdots|d_{L}\rangle_{L}^{j}$, for the right block is $|1\rangle_{R}^{j},|2\rangle_{R}^{j}\cdots|d_{R}\rangle_{R}^{j}$, and the central two sites are $|1\rangle_{a}^{j}\cdots|m\rangle_{a}^{j}$ and $|1\rangle_{b}^{j}\cdots|m\rangle_{b}^{j}$. Here $m$ encodes the size of the Hilbert space for a single site: $m=3$ or 4 when we truncate to 2 or 3 particles per site. In the absence of any approximations $d_{L}=m^{j}$ and $d_{R}=m^{L-j-2}$.

The bases for the various partitions are related by tensors $\Gamma^{L}$ and $\Gamma^{R}$,

[TABLE]

Here $i$ indexes the basis states of the Left/Right block in partition $j$, while $k$ runs over the basis in partition $j-1$ or $j+1$, while $\alpha$ and $\beta$ index the basis states for the $a$ and $b$ site in partition $j$. In the absence of approximations, the tensors are unitary,

[TABLE]

and similar with $L\to R$.

A wavefunction in the $j$’th partition, can be written as a rank 4 tensor of size $d_{L}\times d_{R}\times m\times m$,

[TABLE]

Using the relationship between the basis functions, one has

[TABLE]

S1.2 Approximation

In the DMRG one truncates the basis used for each block, keeping only the $d$ “most important” basis states, where $d$ is referred to as the bond dimension. The “important” states correspond to the $d$ largest eigenvalues of the reduced density matrix for that block. For example, given a wavefunction $\psi^{(j)}_{\ell r\alpha\beta}$, one can construct a reduced density matrix for the left block by tracing out the other degrees of freedom,

[TABLE]

The eigenvalues of this matrix are denoted $\lambda_{i}$, and $\sum_{i}\lambda_{i}=1$.

Rather than storing the basis vectors, in the DMRG one stores matrix elements of the various operators which are needed to construct the Hamiltonian. One only stores operators which act on a single block – which in the truncated basis requires at most $d^{2}$ numbers.

The DMRG is a systematic way of generating a sequence of truncations. In the next few paragraphs we explain how one generates a new, improved truncation from a prior non-optimal truncation.

Consider the $j$’th configuration of the non-optimal truncation. One uses the stored information to generate the Hamiltonian as a $(d^{2}m^{2})\times(d^{2}m^{2})$ matrix. One finds the eigenvector of this matrix corresponding to the smallest eigenvalue, which yields an approximate $\psi^{(j)}_{\ell r\alpha\beta}$. One then finds the reduced density matrix associated with the leftmost $j+1$ sites,

[TABLE]

where the $m\times d$ pairs $(\ell,\alpha)$ provide indices for the basis state of these $j+1$ sites. One then finds the $d$ largest eigenvalues, and the corresponding wavefuctions, $\phi^{s}_{\ell\alpha}$, with $s=1,2,\cdots d$. These become the basis states for left block of the $j+1$’st configuration. These states also act as the transformation matrix

[TABLE]

and are used to transform the matrix elements of the operators which either act on the left block, or the $a$ site. These new matrix elements are stored. Note, the constructed $\Gamma$ tensors are no longer unitary, but rather act as projection operators. This new truncation is better than the old.

The process can then be repeated, using the matrices from the $j+1$’st configuration to update the left-block matrices for the $j+2$’nd configuration. One “sweeps” from left to right until all configurations are updated. Sweeping right to left then updates the right-block matrices. These sweeps are repeated until a fixed point is reached.

There are several tricks to speed up the calculation. First, symmetries and conservation laws (such as number conservation in our problem) imply that the $(d^{2}m^{2})\times(d^{2}m^{2})$ Hamiltonian matrix is block diagonal. One need only diagonalize the block corresponding to the physical quantum numbers. Second, one can use iterative algorithms for the diagonalization. The wavefunctions from prior steps in the sweep can be used to find good starting points for the iterations. As emphasized in the main text, appropriate boundary conditions can also eliminate degeneracies or near degeneracies, and accelerate convergence.

The approximate wavefunction generated by this technique is a matrix-product-state, and the DMRG algorithm can be interpreted as a means of finding the optimal variational matrix-product-state wavefunction mps .

S1.3 Cold gas microscopy

In cold atoms one has a unique probe where one counts how many particles are on each lattice site Ott ; Bakr ; Bakr2 ; Sherson ; Endres . This is a projective measurement which not only tells one about the average occupation of sites, but also the fluctuations and their correlations.

All measurement processes are stochastic, and in each run of the experiment one will see a different image in the microscope. We can use the information stored in our DMRG calculation to stochastically generate configurations which are drawn from the same distribution as the experimental images.

In particular, we need the following: The set of transformation matrices, $\Gamma^{Lj}_{\ell k\alpha}$, $\Gamma^{Rj}_{rk\alpha}$, the wavefunction in the first configuration $\psi^{(1)}_{\ell r\alpha\beta}$, and matrix elements for each of the operators which project into the space where a given number of particles sit on a given site. We let $\hat{P}^{k}_{n}$ be the operator that projects into the space where the $k$’th site has exactly $n$ particles. If $2\leq k\leq L-2$, then this operator acts on the $a$ site in the $j=k-1$’th configuration. If $k=1,L-1$, or $L$, it respectively acts on the left block of the $j=1$’th configuration, the $b$ site of the $j=L-3$ configuration, and the right block of the $j=L-3$ configuration.

We first calculate the probabilities that the first site has $n$ particles on it – for $n=0,1,...n_{\rm max}$, with $n_{\rm max}=2$ or $3$, depending on the calculation. These probabilities are extracted from the wavefunction via

[TABLE]

We then use a random number generator to choose one of these options with the appropriate probability. Denoting our choice $n_{1}$, we project our wavefunction into that space via

[TABLE]

We then use this new wavefunction to calculate the probability for the various values of $n_{2}$, conditioned on the measured $n_{1}$. For notational simplicity, we suppress reference to the prior measurements, writing,

[TABLE]

Again we use a random number generator to choose one of these options with the appropriate probability, and project via

[TABLE]

Next we use the transformation matrices to convert this wavefunction to the next sector,

[TABLE]

We then use this wavefunction to calculate the probabilities for various values of $n_{3}$, conditioned on the prior measurements. One continues in this manner until all occupation numbers are determined.

The reason that we shift between the configurations is so that our projectors always act on blocks which are not truncated, and therefore introduce no approximations. One may be concerned that with our truncated basis, transformations such as Eq. (S12) are not unitary. We verify that the norm of our wavefunction does not decrease as long as we have a sufficiently large bond dimension. Physically, the ultimate state that we project into is “close enough” to the grounds state that it can be well-approximated using the truncated basis.

S2 Higher Band Effects

In our approach, the Haldane phase appears at relatively low lattice depths, $V_{0}\sim 2E_{R}$. Here we model higher band effects, and show that they are negligible.

Using similar notation as in the main text, we consider a two-band tight binding model,

[TABLE]

where the bars label the second band, and the tilde labels the interaction between these two bands. These coefficients are connected with the Wannier functions for the two bands $w_{1}(x)$ and $w_{2}(x)$ via Spielman

[TABLE]

where $H_{0}=p^{2}/(2m)+V_{0}\cos(2k_{L}x)/2$ is the Hamiltonian of one spin component in the absence of the Raman laser with Rabi frequency $\Omega$. The laser recoil energy is $E_{R}=\hbar^{2}k_{L}^{2}/(2m)$, and $U_{0}$ denotes the reduced interaction energy, related to the scattering length and the transverse confinement. We write the Wannier functions in terms of Mathieu functions, and numerically perform the integrals.

We fix $t^{\prime}/t=0.3$, $U_{2}/t=2.5$, and $V_{0}=2E_{R}$. We then obtain $\bar{t}/t=-1.00$, $\bar{t}^{\prime}/t=1.07$, $U/t=3.95$, $\bar{U}/t=2.53$, $\tilde{U}/t=1.91$, $\bar{U}_{2}/t=2.48$, $\tilde{U}_{2}/t=3.45$, $\lambda/t=0.70$, $\delta/t=4.84$, $\Omega/t=4.23$, $U_{0}/t=10.00$.

Of particular note is the fact that for these parameters, the interaction between neighboring atoms in the same band ($U_{2}$ and $\bar{U}_{2}$) are smaller than the interactions between atoms in different bands, $\tilde{U}_{2}$. This feature has to do with the shape of the Wannier states for each bands, and will tend to discourage states where there is occupation of both bands.

We carry out DMRG simulations of this two-band model, using the same boundary conditions as in our single band calculations. We find results which are nearly indistinguishable from our one-band model. To further illustrate this equivalence, we extract order parameters

[TABLE]

Here the particle number is defined as the total number $\tilde{n}_{i}=n_{i}+\bar{n}_{i}$. Our two-band results are shown on the right in Fig. S1 with the single band results on the left. Despite the relatively shallow potential, these results are nearly identical.

Previous studies of higher band effects have seen similar results. For example, Xu et al. used a Quantum Monte-Carlo approach to solve the problem of one-dimensional (1D) interacting Bosons in a 1D optical lattice xu . Although they found that the relationship between the scattering length and the Hubbard $U$ was renormalized at small lattice depths, the single band Bose-Hubbard model remained predictive in this regime. A recent review of extended Hubbard models contains some discussion of multi-band effects dutta .

S3 Convergence Study with bond dimension

As discussed in section S1, the DMRG can be interpreted as a variational technique, which rapidly optimizes the parameters in a matrix-product-state wavefunction. The size of the matrices used is called the bond-dimension $d$. Physically $\log(d)$ is the maximum amount of entanglement entropy which is captured by this ansatz. This technique is exact in the limit $d\to\infty$. We take chains of length $L=60$, and vary the bond dimension from $d=10$ to $500$. Correlation functions are measured over the central 30 sites. For the calculations with large bond dimensions we use the ITensor package, as it is more efficient than the custom-built code that we used for many of our studies.

Figure S2 uses parameters corresponding to the Haldane Insulator phase. This is a gapped phase, where the entanglement entropy is expected to be modest, and independent of system size. Consequently, convergence is extremely rapid. Figure S3 uses parameters corresponding to a region of the superfluid phase which displays strong medium-range density wave correlations. This phase is gapless, and one expects more modest convergence with bond dimension. Moreover the entanglement entropy should grow with system size, Nonetheless the results appear to be quite accurate even for relatively small $d$.

S4 Convergence studies with system size, and extracting phase boundaries

To understand finite size effects in our simulations, we repeat our calculations for different length chains $L$. We calculate the correlation functions $O_{L}=C_{i=L/4,j=3L/4}$. When $t^{\prime}=0$ and $U_{2}/t=2.5$ we find that we can readily extrapolate to the thermodynamic limit, modeling $O_{L}=O_{\infty}+Ae^{-\kappa L}/L^{\alpha}$, and fitting the various coefficients. Figure S4 shows the resulting extrapolation for the density wave, parity, and string correlation functions. For our extrapolation we use data with $L=32,64,128,256,512$, but for clarity do not plot the $L=512$ data. Notice the rapid convergence with system size, away from the phase boundaries. We find equivalent results if we fix $L$, and fit $C_{ij}=O_{\infty}+Ae^{-\kappa|i-j|}/|i-j|^{\alpha}$.

At non-zero $t^{\prime}/t$, the same procedure works for identifying the MI-HI transition. For the other transitions, however, the results somewhat ambiguous, as the correlation length $\xi=1/\kappa$ becomes larger than the system size well before where $O_{\infty}$ vanishes, signaling a large critical region. This behavior is consistent with the expected Kosterlitz-Thouless type physics, and we do not believe that we can accurately extrapolate to the thermodynamic limit.

We can, however, estimate the locations of the phase boundaries by looking for peaks in the curvature $\chi=t^{2}(\partial^{2}O_{L}/\partial t^{\prime 2}+\partial^{2}O_{L}/\partial U^{2})$. Figure S5 shows these susceptibilities for the string and parity correlation functions.

As discussed in the main text, we also estimate the locations of phase boundaries from the bipartite number fluctuations $D_{j}=\langle N_{i<j}^{2}\rangle-\langle N_{i<j}\rangle^{2}$, where $N_{i<j}=\sum_{i=1}^{j-1}n_{i}$ is the number of particles to the left of site $j$. In figure S6 we show $D_{L/2}$ for $L=256$. Plateaus are indicative of bulk phases, while the regions of sharp variation represent our best estimates of the locations of the phase boundaries. As we change system size, we find these transition regions become sharper, but we have not established if the width vanishes in the thermodynamic limit. Nonetheless the width represents a conservative upper bound to our uncertainty about the location of the phase boundary.