Enhanced repulsively bound atom pairs in topological optical lattice ladders

Recent experiments have demonstrated the utility of ultra-cold atoms in optical lattices to explore the physics of topological quantum systems [1–12]. These systems have band structures characterised by a non-local order parameter resulting in novel global features that are in a separate classification from conventional phases [13, 14]. While single-particle properties are generally well understood and have recently been measured experimentally [15–19], there are still many open questions relating to interacting quantum systems in these band structures, questions that cold-atom systems are perfect for exploring [20–30]. In this work, we investigate interacting bosons in a topological band structure where the single-particle kinetic energy is completely frustrated [31–44], and find that the topology enhances the formation of bound pairs allowing them to remain stable for higher temperatures. Going beyond the regime of perturbative interactions we find that pair superfluid phases can be engineered, prepared and detected in current optical-lattice experiments. This opens up ways of exploring the complex interplay between topological band structures and strongly interacting systems allowing for investigations into the effects on the many-body phases and on the resulting dynamical properties.

Specifically, we analyse the properties of bosons in a Creutz ladder (shown in figure 1), which is characterised by complex tunnelling amplitudes along the legs of the ladder while also having diagonal tunnelling components between the legs [45]. In this system geometrical frustration results from the combination of these tunnelling terms where there is a destructive interference effect that completely suppresses the single-particle kinetic energy and gives rise to flat energy bands. However, it has been previously shown that including an onsite interaction can lead to the formation of bound pairs that are stable even for infinitesimal interaction strength which now have dispersion completely dictated by the interactions [35–40, 46–56]. There is growing interest in repulsively interacting bound pairs in general cold-atom systems [57, 58] as well as in more novel band structures [59–64], but pairs are usually only stable for large interaction strengths compared to the tunnelling. By analysing the dispersion relation for single bound pairs in the Creutz ladder beyond the limit of weak interactions we find that in contrast to those formed in conventional lattices, the pair kinetic energy grows with increasing interaction strength. This shows that the topology of the lattice plays an important role and enhances the formation and properties of these pairs.

Zoom In Zoom Out Reset image size

Previous ground state analysis of these systems has identified many-body phases where the correlations are dominated by superfluidity of these pairs [35, 36]. Here we study the excitation spectrum and investigate the robustness of pair correlations to temperature and to excitations in time-dependent preparation. This connects directly to questions of the temperature-dependence of bound fermions in flat-band geometries, which has been recently discussed in the context of topological superconductions [40, 65]. Furthermore, for high densities we identify new lattice supersolid phases [66–76], corresponding to the coexistence of a charge-density wave (CDW) and superfluidity, but where there is no single-particle superfluidity, but rather only pair superfluid correlations. Additionally, we offer new perspectives in the ability to prepare and detect these phases, by first proposing an experimental preparation scheme for a pair condensate using adiabatic manipulations of the optical-lattice potential [77, 78], which can be achieved in timescales that are reachable in current experiments. Finally we consider experimental detection through measurements of the dynamics induced after a local quench through calculations of the dynamical structure factors, where we find substantial qualitative differences between each phase indicating that they can be resolved experimentally.

In figure 1(a) we include the Creutz ladder which has a topological band structure characterised by bands with a non-trivial Zaks phase and (for these tunnelling phases) only flat energy bands, see figure 1(b). This suppression of the dynamics arise through a destructive interference effect as atoms in the Wannier states attempt to travel through the lattice. There are various ways to produce this experimentally [45, 79, 80], where most proposals utilise a synthetic dimension [81–89] where each leg of the ladder corresponds to two atomic internal states. However including interactions will then require careful tuning of the inter-component and intra-component strengths. If there are non-zero inter-component or unequal intra-component interactions then this will result in additional terms appearing [90], complicating the simple bound state picture presented below, however if the imbalance is small, will lead to qualitatively the same features which we demonstrate explicitly below. We also propose an alternative realisation, requiring only atoms in a single internal state confined in a dimerised ladder optical potential, so is not affected by these concerns. All the tunnelling amplitudes, with the correct phase relations, are produced from only two applied fields to facilitate multiple two-photon Raman-assisted tunnelling processes (see appendix A).

Explicitly the Hamiltonian including an onsite two-body contact interaction, U_A for the A-sites and U_B for the B-sites, is given by (ℏ = 1),

$$\begin{equation}\begin{aligned}\hfill H=& \sum _{n}J\left[{\hat{b}}_{n}^{{\dagger}}{\hat{a}}_{n+1}+{\hat{a}}_{n}^{{\dagger}}{\hat{b}}_{n+1}+i{\hat{b}}_{n}^{{\dagger}}{\hat{b}}_{n+1}-i{\hat{a}}_{n}^{{\dagger}}{\hat{a}}_{n+1}+\text{h.c}.\right]\hfill \\ \hfill & +\sum _{n}\left[\frac{{U}_{\text{A}}}{2}{\hat{a}}_{n}^{{\dagger}}{\hat{a}}_{n}^{{\dagger}}{\hat{a}}_{n}{\hat{a}}_{n}+\frac{{U}_{\text{B}}}{2}{\hat{b}}_{n}^{{\dagger}}{\hat{b}}_{n}^{{\dagger}}{\hat{b}}_{n}{\hat{b}}_{n}\right],\hfill \end{aligned}\end{equation} \tag{ 1 }$$

where ${\hat{a}}_{n}^{{\dagger}}$ (${\hat{b}}_{n}^{{\dagger}}$) creates a particle on the A (B) site in the nth unit cell. In order to analyse the bound states in this system, it is advantageous to apply a basis transformation to the local Wannier basis that diagonalises the single-particle Hamiltonian. This basis is shown in blue in figure 1(a) and the transformation is,

$$\begin{equation}\begin{aligned}\hfill {\hat{a}}_{n}& =\frac{1}{2}{\hat{W}}_{n}^{+}+\frac{1}{2}{\hat{W}}_{n}^{-}-\frac{i}{2}{\hat{W}}_{n+1}^{+}+\frac{i}{2}{\hat{W}}_{n+1}^{-}\hfill \\ \hfill {\hat{b}}_{n}& =-\frac{i}{2}{\hat{W}}_{n}^{+}-\frac{i}{2}{\hat{W}}_{n}^{-}+\frac{1}{2}{\hat{W}}_{n+1}^{+}-\frac{1}{2}{\hat{W}}_{n+1}^{-},\hfill \end{aligned}\end{equation} \tag{ 2 }$$

where the ${\hat{W}}_{n}^{{\pm}}$ annihilate a boson at unit cell n in the higher/lower band. This transformation allows us to explicitly see that the single-particle dynamics are suppressed,

$$\begin{align}\hfill H=& -2J\sum _{n}{\hat{W}}_{n}^{-{\dagger}}{\hat{W}}_{n}^{-}+2J\sum _{n}{\hat{W}}_{n}^{+{\dagger}}{\hat{W}}_{n}^{+}+\frac{1}{8}\sum _{n}\left[\frac{{U}_{\text{A}}+{U}_{\text{B}}}{4}{\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}{\hat{W}}_{n}\right.+\frac{{U}_{\text{A}}+{U}_{\text{B}}}{4}{\tilde {W}}_{n}^{{\dagger}}{\tilde {W}}_{n}^{{\dagger}}{\tilde {W}}_{n}{\tilde {W}}_{n}\hfill \\ \hfill & +\left({U}_{\text{A}}+{U}_{\text{B}}\right){\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}{\tilde {W}}_{n+1}^{{\dagger}}{\tilde {W}}_{n+1}-\frac{{U}_{\text{A}}+{U}_{\text{B}}}{4}{\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}^{{\dagger}}{\tilde {W}}_{n+1}{\tilde {W}}_{n+1}+\text{h.c}.-i\frac{{U}_{\text{A}}-{U}_{\text{B}}}{2}{\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}{\tilde {W}}_{n+1}+\text{h.c}\hfill \\ \hfill & \left.-i\frac{{U}_{\text{A}}-{U}_{\text{B}}}{2}{\hat{W}}_{n}^{{\dagger}}{\tilde {W}}_{n+1}^{{\dagger}}{\tilde {W}}_{n+1}{\tilde {W}}_{n+1}+\text{h.c}.\right],\hfill \end{align} \tag{ 3 }$$

where ${\hat{W}}_{n}={\hat{W}}_{n}^{+}+{\hat{W}}_{n}^{-}$ and ${\tilde {W}}_{n}={\hat{W}}_{n}^{+}-{\hat{W}}_{n}^{-}$.

This non-local Wannier function basis makes two novel features of this system apparent: firstly it illustrates the vanishing single-particle kinetic energy upon diagonalising the single-particle Hamiltonian, and secondly it has illuminated the existence of strong pair-tunnelling terms, ${\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}^{{\dagger}}{\tilde {W}}_{n+1}{\tilde {W}}_{n+1}$, as well as nearest-neighbour interactions, ${\hat{W}}_{n}^{{\dagger}}{\hat{W}}_{n}{\tilde {W}}_{n+1}^{{\dagger}}{\tilde {W}}_{n+1}$, which are both proportional to the onsite interaction strength. The effect of an imbalance in the interactions, U_A ≠ U_B, is to introduce complex single-particle density-assisted tunnelling terms.

If we consider the case of, U ≡ U_A = U_B and only two particles, then we can apply an additional basis transformation into a two atom bound state picture. In this basis we find that the Hamiltonian can be diagonalised exactly (see appendix B). In figure 1(c) we plot the resulting dispersion of the lowest bound state band, J_eff (blue) as a function of the onsite interaction strength U where there is an asymmetry between repulsive and attractive interactions, which arises through different couplings between the two single-particle bands depending on the sign of the interactions. For increasing attractive interactions the kinetic energy increases to a maximum and then begins to decrease again, agreeing with previous predictions for excitations on top of a fermion background in the Creutz ladder [40]. This is in contrast to the case of repulsive interactions where the kinetic energy asymptotically approaches a value very close to the single-particle tunnelling amplitude. In figure 2 we plot the energy spectrum for repulsive and attractive onsite interaction strength, U, where changing the sign of the interactions inverts the dispersion relation meaning that the lowest energy band in the repulsive case corresponds to that of the highest band in the attractive interacting case. There are also overlaps between two qualitatively different sets of energy bands, one set dispersive, corresponding to a bound state where each single-particle is in a Wannier state centred on the same unit cell, and a set of dispersionless bands, corresponding to each single-particle centred on a nearest neighbouring unit cell. Note that for a single bound state, there are no terms that can mix these two types of state, allowing us to consider them in isolation. We will see below that in the many-body case there is mixing between these two states, the strength of which varies with density, resulting in a complex phase diagram.

Zoom In Zoom Out Reset image size

The asymmetry in figure 1(c) arises from a qualitative difference in the dispersion relations for the lowest dispersive bands (blue) depending on the sign of the interactions, where for attractive interactions a Dirac cone forms in the lower bands, see figure 2(a), at the point of maximum bound state kinetic energy and then further increasing the strength of the interaction the width of this band begins to decrease due to the strong interactions mixing states within the two single-particle bands, explaining the kinetic energy observations for attractive interactions. For repulsive interactions the lowest band is not affected by this mechanism, and although there are still strong mixing between the single-particle bands, this only imposes an upper bound to the dispersion of the bound states making it possible to realise a large kinetic energy for a wide range of interaction strengths.

In figure 1(c) we also include the effective tunnelling for pairs in a simple lattice (purple), which can be calculated in the limit that U ≫ J through second order perturbation theory and is J²/U [58]. We only expect these conventional pairs to be stable in the region |U| > 10J, and we can see for repulsive interactions in this range that the topologically bound pairs in the Creutz ladder have a kinetic energy that is larger than conventional pairs by nearly an order of magnitude. This difference in behaviour arises because the pairing in conventional lattices occurs with both atoms on a single site, whereas in the case of the Creutz ladder, both atoms exist in the same Wannier basis function, shown in figure 1(a), allowing components to exist where each atom is bound but on different sites. These offsite binding components can then tunnel through the system due to a strong interaction induced nearest neighbour coupling between Wannier states centred at different unit cells, see equation (3). This has consequences for the critical temperature for superfluidity, which because it is proportional to the tunnelling amplitudes, means that the temperatures required to produce a superfluid with the pairs in this system are similar to those needed for a single-particle condensate and is an order of magnitude larger than those needed to prepare a superfluid consisting of pairs in a conventional lattice.

Even though these bound pairs exist in a novel effective topological band structure, see figure 2, the many-body phases will be dominated by pairs condensing in the lowest bands, hence our focus on the kinetic energy of this branch. While in principle there are interesting aspects to the higher energy bands, particularly the Dirac cone at U = ±4J, we focus on the low-energy properties of the repulsively bound pairs and leave the analysis of these higher-energy features for a future study. Nevertheless, we will see below that the enhanced properties of the lowest effective two-atom energy band leads to shorter timescales for the dynamical features, for example allowing for the preparation of states with (quasi) long range pair correlations in experimentally feasible timescales. Note that previous studies on the properties of interacting pairs in other flat band systems do not observe a strong enhancement of kinetic energy [46], even in those that are also characterised by a topological invariant [51]. This indicates that the enhancement is not a general feature for flat band systems but is unique to the Creutz ladder geometry due to the unique form of the underlying Wannier basis, which are completely localised on only two neighbouring unit cells.

We now consider the many-body bosonic case and characterise the phases that are manifested by the pairs as we vary the density. For the moment we restrict the study to interaction strengths that are symmetric U ≡ U_A = U_B and that are weak compared to the separation of the single-particle energy bands, U ≪ 4J, by employing a perturbative Schrieffer–Wolff transformation [35, 36, 38, 91]. This approximation allows us to focus on the lowest flat band in isolation but qualitatively preserves the main features of the full model.

We variationally calculate the ground state directly in the thermodynamic limit using a matrix product state (MPS) algorithm that assumes an infinite and uniform ansatz [92], note that we increase the local dimension such that all bosonic fluctuations are captured. We include a chemical potential term, $\left(2-\mu \right){\sum }_{n}{\hat{W}}_{n}^{-{\dagger}}{\hat{W}}_{n}^{-}$, in equation (3) and calculate the pair correlations in the site basis, $\langle {\hat{a}}_{r}^{{\dagger}}{\hat{b}}_{r}^{{\dagger}}{\hat{a}}_{0}{\hat{b}}_{0}\rangle$, for a range of densities controlled by the ratio of the chemical potential to the onsite interaction strength, μ/U. In all cases the single-particle correlation function is exponentially suppressed, reflecting the lack of single-particle dispersion in this system. Additionally, we find that these off-site pair correlation functions dominate over terms of the form, $\langle {\hat{a}}_{r}^{{\dagger}}{\hat{a}}_{r}^{{\dagger}}{\hat{a}}_{0}{\hat{a}}_{0}\rangle$, reflecting that the pairs exist in a superposition of separate sites, in contrast to conventional repulsively bound doublons [58]. The different phases are then characterised by algebraically or exponentially decaying pair correlations for the pair superfluid phase, which we refer to as a pair Luttinger liquid (PLL), and the pair CDW phase respectively, examples of which we have included in figures 3(a) and (b). The phase diagram is shown in figure 3(c) where we have highlighted the PLL phases and the CDW phases and we have also plotted the value of the CDW order parameter [93, 94], which is given by, ${O}_{\text{DW}}={\mathrm{lim}}_{r\to \infty }{\left(-1\right)}^{r}\langle \delta {\hat{n}}_{r}\delta {\hat{n}}_{0}\rangle$, where $\delta {\hat{n}}_{i}={\hat{W}}_{i}^{-{\dagger}}{\hat{W}}_{i}^{-}-\rho$. Note that analysing the density–density correlation functions (with the exception of the ρ < 1 PLL phase, see appendix C) the component that oscillates with wave-vector ka = π is strongly dominant. This is due to the strong nearest neighbour interaction term in equation (3) which, as there are no longer range interactions present, ensures that this commensurate order is stabilised for all densities considered.

Zoom In Zoom Out Reset image size

We agree with the predictions of reference [35, 36] which analyse the same system but restrict their analysis to low densities, where we find a PLL phase for densities between 1/2 < ρ < 1 (per unit cell) and a phase transition to a CDW for ρ = 1. We then investigate larger densities where we find a second PLL phase, indicating that the ρ = 1 CDW is unstable to the addition of more pairs. And upon further increasing the density we find large regions at incommensurate density where distinct phases exist that share features of both the PLL and CDW. We denote these phases lattice pair supersolid (PSS) and are characterised by algebraically decaying pair correlations but with a non-zero density wave order parameter. Note that if we go beyond the weakly interacting regime such that, U ∼ J, the coexistence of these phases is suppressed and we either have a PLL or a CDW phase with a clear phase transition point.

Now we include an imbalance between the onsite interaction strengths, U_A ≠ U_B, and perform the same analysis on the ρ < 1 PLL phase. We plot the resulting single-particle and pair correlations in figure 4 where we see that there are still dominant algebraically decaying pair correlations, but now there are exponentially decaying single-particle correlations. For values of U_A/U_B close to 1, the pair correlations greatly dominate over the single-particle correlations, indicating that the novel properties of this phase survive well into the imbalanced interaction regime. We can see that as U_A/U_B decreases the single-particle correlations decay with a smaller correlation length indicating the onset of a conventional superfluid phase for U_A/U_B → 0.

Zoom In Zoom Out Reset image size

Here we present a scheme to prepare a many-body eigenstate with strong pair superfluid correlations which can be achieved in a cold-atom experiment by varying the relative intensity of the lasers that create the optical potential. We begin the experimental sequence by applying a large dimerisation to the optical potential such that sites that are populated with atoms are at a much lower energy than neighbouring sites leading to atoms that are strongly localised. We then adiabatically vary the optical potential in order to slowly remove this dimerisation, which amounts to a ramp of onsite energy and allows the atoms to gradually delocalise throughout the time-dependent ramp which then prepares the eigenstate of the final Hamiltonian if the ramp time is long enough [77, 78] (see appendix E for details). We consider the case of large interaction strength, U ≡ U_A = U_B = 6J, so that we have a large pair kinetic energy (see figure 1(c)), allowing the correlations to spread to the entire system in timescales that are sufficiently fast so that we can ignore heating and dissipation effects. In principle, once a condensate has been prepared, we can ramp the interaction strength to weak values in order to prepare the phases predicted in the previous section.

In figure 5 we consider a ladder of M = 192 sites and for N = 72 bosons and plot the results of this process as the total ramp time, TJ, is varied. It is clear that we can produce a many-body state with significant pair correlations, $\langle {\hat{a}}_{r}^{{\dagger}}{\hat{b}}_{r}^{{\dagger}}{\hat{a}}_{0}{\hat{b}}_{0}\rangle$ and vanishing single-particle correlations, $\langle {b}_{r}^{{\dagger}}{b}_{0}\rangle$ and $\langle {a}_{r}^{{\dagger}}{a}_{0}\rangle$, in experimentally feasibly timescales (TJ = 210). However, as we are attempting to prepare a phase that is gapless in the thermodynamic limit we expect that as we increase the system size that the total ramp time to achieve the same level of correlation decay will continue to increase. The analysis for short ramp times (TJ = 40) indicates that for larger systems we could still produce a pair superfluid in experimentally achievable timescales with the cost of introducing effective finite size effects in the correlations.

Zoom In Zoom Out Reset image size

We also consider the effects of a finite temperature on the pair correlations in the system and we use an imaginary time MPS algorithm which utilises a purification of the density matrix to calculate the state at a given temperature [95, 96]. In figure 5(b) we plot the pair correlations at varying temperature and find that the correlations are exponentially suppressed at high temperatures ($\tilde {T}{ >}20J$) where they become numerically indistinguishable from the exponentially small single-particle correlations. Note that at short distances the single-particle correlations remain qualitatively the same as those of the ground state (see figure 5(a)) where they dominate over a distance of several unit cells, simply because the particles are spread into the Wannier basis states. For intermediate temperatures ($2J{< }\tilde {T}{< }15J$), the correlations begin to approximate those of the ground state at short distances but still decay exponentially, with a correlation length that grows as the temperature is decreased. And for temperatures $\tilde {T}{< }2J$ the pair correlations very closely match the zero temperature case, indicating that these pair superfluid phases are robust to finite temperatures. There are small discrepancies for the longer range tails but these will be experimentally indistinguishable.

We move on and consider experimental detection of these phases. One possible way to do this is through measurements of the time-dependent onsite particle number after a local quench, in this case after the application of the number operator on a single unit cell. To this end we have calculated the dynamical structure factors for each phase (within the weakly interacting and isolated flat band limit) using an MPS algorithm for time evolving infinite systems after a local perturbation [97, 98]. Explicitly we calculate the unequal time two-point correlator, $A\left(r,t\right)=\langle {{\Psi}}_{0}\vert \delta {\hat{n}}_{r}\left(t\right)\delta {\hat{n}}_{0}\left(0\right)\vert {{\Psi}}_{0}\rangle$, where |Ψ₀⟩ is the initial state, $\delta {\hat{n}}_{r}\left(t\right)={\hat{W}}_{r}^{-{\dagger}}\left(t\right){\hat{W}}_{r}^{-}\left(t\right)-\rho$ where ρ is the density. We then take the Fourier transform which we plot in figure 6 for the different phases. Note that in all cases, U_A = U_B.

Zoom In Zoom Out Reset image size

We can see for the PLL phase (a) there is a dominant linear excitation originating from the ka = π mode matching the predictions from Luttinger liquid theory (see appendix C) and for the CDW phase (b) there are well defined gapped energy branches, which can be interpreted as collective quasi-particle excitations. For the lattice supersolid phases (figures 6(c) and (d)) we can see gapless excitations that are linear for low energies and originate from ka = 0 and ka = π which is set by the underlying CDW order of the ground state, but there are also sinusoidal excitation branches present at higher energies similar to those in the CDW phase indicating the coexistence of the two phases also exists in the excitation spectrum. In figures 6(c) and (d) we also observe in the lower energy branches a breaking of translation symmetry as there is a doubled periodicity in k space indicating that the low energy superfluid features exist on top of a dimerised ground state. Additionally in figure 6(d) there is a local minimum in the low energy dispersion which resembles the helium roton spectrum [99–101], with a roton mode local minimum close to zero gap. These calculations clearly illustrate that the distinct features of each phase are manifested in the excitation spectrum offering a way to experimentally resolve the phases through measurements of the dynamics produced after a local quench.

We have considered the experimental opportunities of using the Creutz ladder to investigate the interplay between topological band structures and strong interactions. By analysing the properties of single repulsively bound pairs we found that the topology greatly enhances the stability and kinetic energy of formed pairs making it possible to realise and investigate pair superfluid phases in experiments with cold-atoms. We considered the ability to prepare and detect these phases where we illustrated an experimentally feasible preparation scheme allowing us to prepare a pair superfluid in realistic timescales and demonstrated that these phases can be resolved through measurements of the dynamical properties. This opens up opportunities for understanding and exploring the unique many-body phases that can be produced with strong interactions in more general topological band structures [40].

While previous analysis on other systems with flat energy bands have found similar superfluid phases dominated by bound pairs that share qualitative features as the pair superfluid phases presented here, such as in the sawtooth and Kagome lattice [38] or the diamond chain [56] they have not been considered in the context of practicalities for experimental detection. In particular, the latter system must have a very precise value for the tunnelling phases in order for it to manifest a pair superfluid, in contrast to the robustness of these phases in the Creutz ladder to similar phase deviations [36] and even to an interaction imbalance or finite temperatures as analysed in this paper. The stability of these features in the Creutz ladder is due to the enhanced properties for single repulsively bound pairs, which in contrast to the sawtooth lattice [46] or the diamond chain [51], we find have a kinetic energy that grows with increasing interaction strength, due to the particular topology of the Creutz ladder. This distinct behaviour for the Creutz ladder opens up many questions relating to the features of other flat band systems, in particular those with additional dispersive energy bands [65].

This work was supported in part by the EPSRC Programme Grant DesOEQ (EP/P009565/1), and by the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 817482 PASQuanS. SF acknowledges the financial support of the Carnegie trust. Results were obtained using the ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk) based at the University of Strathclyde. All data underpinning this publication are openly available from the University of Strathclyde KnowledgeBase at https://doi.org/10.15129/81ecb9ec-7701-4940-b496-9f2785079198.

Experimental realisations of the Creutz ladder have been proposed in references [45, 79], but require the manipulation of two internal states of an atom where each represent each leg of the ladder. In order to produce the complex tunnelling terms in these methods we need multiple Raman-assisted tunnelling processes and Floquet driving elements. The main barrier to these approaches for realising the physics analysed in this work is that the use of a synthetic dimension, with two atomic internal states, would in most cases lead to strong inter-component interactions, U_AB and/or an asymmetry between the intra-component interactions for each species, U_AA and U_BB. The former would lead to a strong nearest-neighbour interaction between sites in the unit cell and the latter would result in density-assisted single-particle tunnelling terms appearing between nearest-neighbour unit cells in our model. These additional processes would complicate the simple bound state picture that we discussed in the main text. However if an atomic species could be found that has a zero crossing for the inter-component interaction, such that U_AB ≈ 0 in the vicinity of a Feshbach resonance for the intra-component interaction, which results in a small interaction strength asymmetry relative to the overall magnitude, U_AA − U_BB ≪ U_AA, then we explicitly illustrated that the main features of the phases survive. If these considerations can be satisfied then these schemes present viable options for realising the many-body phases proposed here.

Here we propose a different experimental realisation with only a single internal atomic state, allowing us to easily satisfy the considerations above. Our scheme requires two fields to facilitate all Raman-assisted tunnelling processes and a two-site superlattice see figure 7. In order to produce this, we find it convenient to redefine the tunnelling amplitudes through a gauge transformation, and we include the resulting tunnelling terms in figure 7(a). This ladder with new tunnelling components results in the same topological physics as the ladder considered in the main text because the phases accumulated when moving around loops in the lattice are unmodified. Now the tunnelling elements along each leg of the ladder are real and have the same phase while it is the diagonal tunnelling terms that carry the complex phase factors. Notice also that there is now a dimerisation of the diagonal phase resulting in a doubled periodicity and a larger unit cell and therefore a different local Wannier function basis. The new basis is shown in figure 7(b) and is very similar to the one used in the main text where it is also perfectly localised to two unit cells as before. This new basis results in the same single-particle spectrum and transforming the many-body Hamiltonian into this basis results in the same model that is analysed in the main text. We confirm that the Creutz ladder shown in figure 7(a) is able to quantitatively reproduce all features presented in the main text, while also offering a more viable experimental implementation.

Zoom In Zoom Out Reset image size

For this scheme we require a particular separation of the onsite energy levels, shown in figure 7(d) (note that there is a flexibility in the energy differences). This energy level separation can be produced with a simple 2D superlattice potential, however, additional potential barriers must be applied to ensure that there is no tunnelling between sites A–B and C–D.

We then apply a Raman-assisted tunnelling process, requiring two fields (figure 7(e)), to create the two complex diagonal tunnelling processes and the real processes along the legs. This can be achieved with a single laser with a sideband allowing the necessary phase relations to be easily enforced. The effects of these applied fields is to induce tunnelling processes between off-resonant sites [102],

$$\begin{equation}{J}_{\alpha }=\frac{{{\Omega}}_{1}{{\Omega}}_{2}^{{\ast}}}{\delta }\int \mathrm{d} \overrightarrow {r}\enspace {\phi }^{{\ast}}\left( \overrightarrow {r}\right){\mathrm{e}}^{\text{i}\delta { \overrightarrow {k}}_{\alpha } \overrightarrow {r}}\phi \left( \overrightarrow {r}-a \overrightarrow {{R}_{\alpha }}\right),\end{equation} \tag{ A1 }$$

where ${\phi }^{{\ast}}\left( \overrightarrow {r}\right)$ are the onsite Wannier functions, $\delta \overrightarrow {k}={ \overrightarrow {k}}_{1}-{ \overrightarrow {k}}_{2}$ is the difference between the wave vectors of the two lasers in a single Raman process, $\vert { \overrightarrow {k}}_{i}\vert =2\pi a/{\lambda }_{i}$, δ is the detuning between ω₁, ω₂ and the excited internal state (see figure 7(e)) and Ω_i is the Rabi frequency of the applied laser with frequency ω_i .

Assuming that the distance between sites, a, is the same in both directions, then the phase factor in the tunnelling amplitudes between each site labelled in figure 7(c) is,

$$\begin{equation}\begin{aligned}\hfill {J}_{\text{CB}}& \propto \mathrm{exp}\left(i\frac{a}{2}\left(\delta {k}_{x,\text{CB}}+\delta {k}_{y,\text{CB}}\right)\right)\hfill \\ \hfill {J}_{\text{BC}}& \propto \mathrm{exp}\left(i\frac{a}{2}\left(\delta {k}_{x,\text{BC}}-\delta {k}_{y,\text{BC}}\right)\right)\hfill \\ \hfill {J}_{\text{DA}}& \propto \mathrm{exp}\left(i\frac{a}{2}\left(\delta {k}_{x,\text{DA}}-\delta {k}_{y,\text{DA}}\right)\right)\hfill \\ \hfill {J}_{\text{AD}}& \propto \mathrm{exp}\left(i\frac{a}{2}\left(\delta {k}_{x,\text{AD}}+\delta {k}_{y,\text{AD}}\right)\right)\hfill \\ \hfill {J}_{\text{DB}}& \propto \mathrm{exp}\left(i\frac{a}{2}\delta {k}_{x,\text{DB}}\right)\hfill \\ \hfill {J}_{\text{CA}}& \propto \mathrm{exp}\left(i\frac{a}{2}\delta {k}_{x,\text{CA}}\right).\hfill \end{aligned}\end{equation} \tag{ A2 }$$

where all terms correspond to tunnelling events from left to right in the lattice, along the directions of the arrows in figure 7(a). The right to left processes are then the complex conjugates.

If we assume that the difference in the frequencies of the two fields in a Raman pulse is much smaller than the magnitudes, |ω₁ − ω₂|≪ ω₁, ω₂, where ω₁ > ω₂ then we can assume that the wavelengths for each component has the same magnitude, λ ≡ λ₁ ≈ λ₂, when calculating the phases appearing in the tunnelling amplitudes. Then assuming each field is applied in a general direction in the x–y plane, results in the phases, (where the angles are given relative to the x-axis)

$$\begin{equation}\begin{aligned}\hfill {\phi }_{\text{CB}}& =\frac{a}{2}\left(\delta {k}_{x,\text{CB}}+\delta {k}_{y,\text{CB}}\right)\approx \frac{a\pi }{\lambda }\left(\mathrm{cos}\enspace {\theta }_{2}-\mathrm{cos}\enspace {\theta }_{1}+\mathrm{sin}\enspace {\theta }_{2}-\mathrm{sin}\enspace {\theta }_{1}\right)\equiv \pi /2\hfill \\ \hfill {\phi }_{\text{BC}}& =\frac{a}{2}\left(\delta {k}_{x,\text{BC}}-\delta {k}_{y,\text{BC}}\right)\approx \frac{a\pi }{\lambda }\left(\mathrm{cos}\enspace {\theta }_{1}-\mathrm{cos}\enspace {\theta }_{2}-\mathrm{sin}\enspace {\theta }_{1}+\mathrm{sin}\enspace {\theta }_{2}\right)\equiv \pi /2\hfill \\ \hfill {\phi }_{\text{DA}}& =\frac{a}{2}\left(\delta {k}_{x,\text{DA}}-\delta {k}_{y,\text{DA}}\right)\approx \frac{a\pi }{\lambda }\left(\mathrm{cos}\enspace {\theta }_{2}-\mathrm{cos}\enspace {\theta }_{1}-\mathrm{sin}\enspace {\theta }_{2}+\mathrm{sin}\enspace {\theta }_{1}\right)\equiv -\pi /2\hfill \\ \hfill {\phi }_{\text{AD}}& =\frac{a}{2}\left(\delta {k}_{x,\text{AD}}+\delta {k}_{y,\text{AD}}\right)\approx \frac{a\pi }{\lambda }\left(\mathrm{cos}\enspace {\theta }_{1}-\mathrm{cos}\enspace {\theta }_{2}+\mathrm{sin}\enspace {\theta }_{1}-\mathrm{sin}\enspace {\theta }_{2}\right)\equiv -\pi /2\hfill \\ \hfill {\phi }_{\text{DB}}& =\frac{a}{2}\delta {k}_{x,\text{DB}}\approx \frac{a\pi }{\lambda }\left(\mathrm{cos}\enspace {\theta }_{2}-\mathrm{cos}\enspace {\theta }_{1}\right)\equiv 0\hfill \\ \hfill {\phi }_{\text{CA}}& =\frac{a}{2}\delta {k}_{x,\text{CA}}\approx \frac{a\pi }{\lambda }\left(\mathrm{cos}\enspace {\theta }_{1}-\mathrm{cos}\enspace {\theta }_{2}\right)\equiv 0,\hfill \end{aligned}\end{equation} \tag{ A3 }$$

where on the right we have shown the values that would produce the Creutz ladder shown in figure 7(a). It is then simply an exercise in finding the optimal direction for the applied fields in order to produce the desired phase differences between each of the tunnelling terms. There is a huge flexibility over the relative angles and the value of the applied λ. The relationship between the angles must satisfy,

$$\begin{equation}{\theta }_{1}=-{\theta }_{2}.\end{equation} \tag{ A4 }$$

This requirement, means that motion in the x-direction does not give rise to any phase change, as cos θ₂ − cos θ₁ = 0. And so all phase dependence comes from moving in the y-direction with sin θ₂ − sin θ₁ = −2 sin θ₁,

$$\begin{equation}\begin{aligned}\hfill {\phi }_{\text{CB}}& =-\frac{2a\pi }{\lambda }\enspace \mathrm{sin}\enspace {\theta }_{1}\equiv \pi /2\hfill \\ \hfill {\phi }_{\text{BC}}& =-\frac{2a\pi }{\lambda }\enspace \mathrm{sin}\enspace {\theta }_{1}\equiv \pi /2\hfill \\ \hfill {\phi }_{\text{DA}}& =\frac{2a\pi }{\lambda }\enspace \mathrm{sin}\enspace {\theta }_{1}\equiv -\pi /2\hfill \\ \hfill {\phi }_{\text{AD}}& =\frac{2a\pi }{\lambda }\enspace \mathrm{sin}\enspace {\theta }_{1}\equiv -\pi /2\hfill \\ \hfill {\phi }_{\text{DB}}& =0\hfill \\ \hfill {\phi }_{\text{CA}}& =0.\hfill \end{aligned}\end{equation} \tag{ A5 }$$

To clarify, the phases picked up from tunnelling (left to right) from site B to C, ϕ_CB, and from site D to A, ϕ_AD, correspond to tunnelling along the same direction (just translated in the x-direction). But, because B to C requires a decrease in energy while D to A requires an increase, this ensures a difference in phase between these components. Similarly, the tunnelling (left to right) from site C to B, ϕ_BC, and tunnelling from site D to A, ϕ_AD, although they both require an increase in energy they involve motion in different directions, so also have a phase difference.

Now the angle θ₁ and θ₂ must be tuned with the wavelength, λ in order to obtain the correct phases. As an example, we plot one possible choice in figure 7(c) where θ₁ = −π/2 and θ₂ = π/2, which would then require a wavelength, λ = 4a. And in figure 8 we include all possible values for these angles and the corresponding values that λ must be set to in order to achieve the phase values in equation (A3). The values for these wavelengths are between 0 and 4a, so are at the same order of magnitude as the lasers responsible for producing the lattice.

Zoom In Zoom Out Reset image size

As a final comment, if it is not possible to realise these exact phase relations (or terms with equal magnitude) for the various tunnelling terms in an experimental setting then this will most likely result in a single-particle energy band structure with bands that are not perfectly flat. If the curvature of the resulting bands are small, compared to the energy band gap and the onsite interaction strength, then some of the many-body phases illustrated in the main text have already been shown to survive into this regime [36]. And we believe that the new phases predicted here should also survive and that our considerations on experimental preparation and detection to still be relevant.

In this section we present the derivation of the effective single-particle Hamiltonian for interacting bound sates, which allows us to analytically calculate the single bound state dispersion relation presented in the main text. If we consider the case where the onsite interaction strengths are uniform, U_A = U_B, then the Hamiltonian in equation (3) does not contain any terms that correspond to the motion of a single-particle. This means that we have two types of bound state that cannot mix with one another. One set of states corresponds to two atoms on the same unit cell,

$$\begin{equation}\begin{matrix}\hfill {\hat{\alpha }}_{i}=\frac{1}{\sqrt{2}}{\hat{W}}_{i}^{-}{\hat{W}}_{i}^{-},\quad {\hat{\beta }}_{i}=\frac{1}{\sqrt{2}}{\hat{W}}_{i}^{+}{\hat{W}}_{i}^{+},\quad {\hat{\gamma }}_{i}={\hat{W}}_{i}^{+}{\hat{W}}_{i}^{-},\hfill \end{matrix}\end{equation} \tag{ B1 }$$

and the second set corresponds to the two atoms on neighbouring unit cells,

$$\begin{equation}\begin{matrix}\hfill {\tilde {\alpha }}_{i}={\hat{W}}_{i}^{-}{\hat{W}}_{i+1}^{-},\quad {\tilde {\gamma }}_{i}={\hat{W}}_{i}^{+}{\hat{W}}_{i+1}^{-}\hfill \\ \hfill {\tilde {\beta }}_{i}={\hat{W}}_{i}^{+}{\hat{W}}_{i+1}^{+},\quad {\tilde {\kappa }}_{i}={\hat{W}}_{i}^{-}{\hat{W}}_{i+1}^{+}.\hfill \end{matrix}\end{equation} \tag{ B2 }$$

Transforming to this basis gives rise to a Hamiltonian containing only quadratic terms.

If we consider the case where we only have two atoms in the system, then the basis defined by equations (B1) and (B2) form an orthonormal basis set and we can solve the system exactly. Explicitly, the momentum space Hamiltonian is,

$$\begin{equation}H=\sum _{k}{{\Psi}}_{k}^{{\dagger}}{H}_{k}{{\Psi}}_{k},\end{equation} \tag{ B3 }$$

where

$$\begin{equation}{{\Psi}}_{k}=\left(\begin{matrix}\hfill {\hat{\alpha }}_{k}\hfill \\ \hfill {\hat{\beta }}_{k}\hfill \\ \hfill {\hat{\gamma }}_{k}\hfill \\ \hfill {\tilde {\alpha }}_{k}\hfill \\ \hfill {\tilde {\beta }}_{k}\hfill \\ \hfill {\tilde {\gamma }}_{k}\hfill \\ \hfill {\tilde {\kappa }}_{k}\hfill \end{matrix}\right).\end{equation} \tag{ B4 }$$

and

$$\begin{align}\hfill {H}_{k}=\left(\begin{matrix}\hfill -4+U/4\left(1-\mathrm{cos}\enspace k\right)\hfill & \hfill U/4\left(1-\mathrm{cos}\enspace k\right)\hfill & \hfill iU\sqrt{2}/4\enspace \mathrm{sin}\enspace k\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill U/4\left(1-\mathrm{cos}\enspace k\right)\hfill & \hfill 4+U/4\left(1-\mathrm{cos}\enspace k\right)\hfill & \hfill iU\sqrt{2}/4\enspace \mathrm{sin}\enspace k\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -iU\sqrt{2}/4\enspace \mathrm{sin}\enspace k\hfill & \hfill -iU\sqrt{2}/4\enspace \mathrm{sin}\enspace k\hfill & \hfill U/2\left(1+\mathrm{cos}\enspace k\right)\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -4+U/4\hfill & \hfill -U/4\hfill & \hfill U/4\hfill & \hfill -U/4\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -U/4\hfill & \hfill 4+U/4\hfill & \hfill -U/4\hfill & \hfill U/4\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill U/4\hfill & \hfill -U/4\hfill & \hfill U/4\hfill & \hfill -U/4\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -U/4\hfill & \hfill U/4\hfill & \hfill -U/4\hfill & \hfill U/4\hfill \end{matrix}\right).\end{align} \tag{ B5 }$$

In figure 2 we plot the energy spectrum as a function of the onsite interaction strength, U, where we can see that at U = 4J the formation and then separation of a Dirac cone, which is a signature of a topological transition. Analysing the topology of these bound states will form part of our future analysis.

As we have a 1D superfluid, we expect to be able to describe the superfluid phases by mapping to a homogeneous Luttinger liquid model [103], (ℏ = 1)

$$\begin{equation}H=\frac{u}{2\pi }\int \mathrm{d}r\left[K{\left(\pi {\Pi}\left(r\right)\right)}^{2}+\frac{1}{K}{\left(\nabla \phi \left(r\right)\right)}^{2}\right],\end{equation} \tag{ C1 }$$

where the bosonic field, ϕ(r), and its conjugate momentum density, Π(r), satisfy the commutation relation, $\left[{\Pi}\left(r\right),\phi \left({r}^{\prime }\right)\right]=i\delta \left(r-{r}^{\prime }\right)$. All low energy properties of Luttinger liquids are completely known once the two parameters, u and K, are obtained, hence the benefit of mapping our phases into this model. In the following we extract these quantities by first fitting an algebraically decaying function to the off-diagonal pair correlation functions to obtain K,

$$\begin{equation}\langle {W}_{r}^{{\dagger}}{W}_{r}^{{\dagger}}{W}_{0}{W}_{0}\rangle \sim {r}^{-1/2K}.\end{equation} \tag{ C2 }$$

We use the pair correlation function because we know that the fundamental particles in the superfluid are pairs, and also that the single-particle off-diagonal correlation function is zero for all distances. Note that the algebraic decay persists for close to 1000 unit cells before numerical errors destroy this behaviour, allowing us to very accurately extract the decay exponents. We then plug the value for K into the expression for the compressibility to obtain u,

$$\begin{equation}\kappa =\frac{\mathrm{d}\rho }{\mathrm{d}\mu }=\frac{K}{u\pi },\end{equation} \tag{ C3 }$$

where we have evaluated the dρ/dμ numerically from our data presented in figure 3 from the main text.

The parameter u is the effective speed of sound in the condensate which is the gradient of the linear dispersion of the excitations. The parameter K controls the thermodynamic properties of the system, for example, it has been shown in reference [104] that if K > 1 then the density transport in the system is completely robust against a single impurity, but for K < 1 the effect of an impurity is to completely suppress transport. And in reference [105] it was shown that K controls the thermal conductivity.

There has also been recent interest in going beyond the assumptions in the Luttinger liquid theory, namely assuming that the excitations follow a linear dispersion relation. We can compute the leading order correction to this, the effective mass, m*, for a non-linear Luttinger liquid theory [106, 107], through,

$$\begin{equation}\frac{1}{{m}^{{\ast}}}=\frac{u}{K}\frac{\mathrm{d}}{\mathrm{d}\mu }\left(u\sqrt{K}\right).\end{equation} \tag{ C4 }$$

Below, we perform this analysis for both pair superfluid phases in density ranges ρ < 1 and then ρ > 1.

C.1.  ρ < 1

We begin with the ρ < 1 superfluid phase and apply the Luttinger liquid formalism. Reference [35] maps this phase onto a spin 1/2 system where |2⟩ = |↑⟩ and |0⟩ = |↓⟩, and while they test the validity of this mapping for small system sizes, we find that the mapped spin 1/2 system does not yield the same physics as the full boson model for the infinite system considered here. We find that we need to account for the possibility of up to four bosons existing on a given site in order to properly account for all the bosonic fluctuations. Our predictions are not qualitatively different, but the critical value for the chemical potential at the phase transition is shifted. However, if this mapping was valid then this commensurate-incommensurate quantum phase transition at ρ = 1 would be mathematically equivalent to that for gapped spin 1/2 chains in a magnetic field [103, 108]. By comparing the correlation functions of the bosonic ground state to these predictions we can quantify the deviations away from the spin 1/2 regime. If the mapping was valid, then we would be able to exactly derive the critical exponents of the phase transition on the incommensurate side simply by knowing the value of the Luttinger liquid parameters at the phase transition point. Because the deviations away from the spin 1/2 regime are small, we can still estimate these quantities.

We extract the Luttinger liquid parameters numerically through the procedure described above and plot these in figure 9(a). The Luttinger liquid parameter K is always less than one, indicating that the superfluid is dominated by charge fluctuations induced by the effective nearest-neighbour interactions, and the value at the phase transition point can be extrapolated, K* ≈ 0.3.

Zoom In Zoom Out Reset image size

We also plot the inverse effective mass, 1/m* in figure 9(a) where we see that for smaller values of the chemical potential it is much smaller than the other parameters, indicating that there are only small corrections to the conventional Luttinger liquid theory. However, for larger chemical potentials—and so larger densities—the values become negative with magnitudes that are larger than the u, signifying that there may be features beyond that of a Luttinger liquid.

The long distance behaviour of the density-density correlation functions for a spin 1/2 system in the presence of an applied magnetic field is given by,

$$\begin{equation}\langle {\hat{n}}_{r}{\hat{n}}_{0}\rangle \sim \mathrm{cos}\left(\pi \left(1+2\delta \rho \right)r\right){r}^{-2K},\end{equation} \tag{ C5 }$$

where δρ is the deviation in the spin magnetisation (or the density for the bosonic system considered here) from the commensurate regime. This correlation function predicts peaks in the Fourier transform (static structure factor) at values of q_± = π(1 ± 2δρ). We calculate these structure factors for our system and plot these in figure 9(c) for a range of μ/U and we find that the peaks do indeed correspond to the values predicted by the Luttinger liquid theory. Even though there are terms present in the model which give rise to an effective nearest neighbour interaction (see equation (3)) these are clearly negligible in this phase as the features can be predicted from the Luttinger liquid theory. This can be understood by realising that the relative magnitudes of the terms in the Hamiltonian will be dependent on density, for example $\langle {W}_{n}^{{\dagger}}{W}_{n}{W}_{n+1}^{{\dagger}}{W}_{n+1}\rangle \approx {\rho }^{2}$, and so the nearest neighbour interactions are dominated by the other processes in this phase. A similar analysis was carried out in reference [37] on a Sawtooth lattice but the commensurate phase is a ρ = 1/2 CDW made stable by the dispersionless energy band and the incommensurate phase occurs for increasing densities. Here we find qualitatively the same features but our gapless commensurate phase is a CDW stabilised by the effective nearest-neighbour interactions.

C.2.  ρ > 1

There is still the question of whether the regimes where the two phases exist simultaneously are really lattice supersolids. In the main text we calculate the energy dispersion relations above the ground state in an attempt to detect supersolid signatures in the excitations, but there is also the question of phase separation. Will the CDW and PLL phases exist uniformly throughout the whole system or will there be distinct regions of one and separate regions of the other? In an attempt to answer this, we consider the ρ < 2 PSS phase and compare the free energy, = E(μ₁)/2 + E(μ₂)/2, with E(μ₀) such that the density, ρ(μ₁)/2 + ρ(μ₂)/2 = ρ(μ₀). Where we choose μ₀ to be in the ρ < 2 PSS phase, μ₁, to be in the ρ > 1 PLL phase and μ₂ to be in the ρ = 2 CDW phase. This quantifies if it is more energy favourable, for a system with a given total particle number, to have regions of lower and regions of higher density, or for it to have uniform density. The results of this analysis are that it is more favourable to have uniform density, indicating that phase separation in the ρ < 2 PSS phase does not occur. However, finite size effects in a real experimental setting may alter this behaviour so future care must be taken.

In cold-atom experiments, if the temperature is much smaller than the critical temperature for condensation then effectively the system is at zero temperature and we can model the dynamics as a pure state. The main consideration in the preparation of low energy eigenstates is in reducing the overall entropy of the many-body state. In adiabatic state preparation [77, 78], this is achieved by first producing a low entropy initial state, in this case by projecting single atoms onto single sites, and ensuring that the state is the ground state of the initial Hamiltonian, in this case by having the onsite energies of populated sites at a much lower energy than the others. This also ensures that the trapped atoms have no dynamics. We then ramp the parameters of the lattice so as to create the final Hamiltonian that we are interested in. And if this ramp process is slow enough, so as not to induce unwanted heating effects, but fast enough so that decoherence effects can still be ignored then we can produce the desired low energy eigenstate with also a low entropy.

We begin with atoms populating only particular unit-cells, where the number of populated unit cells is chosen to give the required density and are equally separated. On a populated unit cell, we have a single atom on each of the two sites, and we have the onsite energies of these populated sites at a significantly lower energy, E₀/J = −μ₀ compared to the rest, E = 0. We then ramp the energy of the populated sites to the value of the other sites, using the following exponential ramp,

$$\begin{equation}E\left(t\right)=-{\mu }_{0}\frac{{\mathrm{e}}^{5\left(1-t/T\right)}-1}{{\mathrm{e}}^{5}-1},\end{equation} \tag{ E1 }$$

where T is the total time for the ramp. We ensure that the initial state is an eigenstate by beginning with a product state where we have the atoms localised on these initial sites and with all nearest-neighbour tunnelling amplitudes set to zero. We then linearly ramp all tunnelling terms from zero to one in a time TJ = 10. The resulting state only has a small amplitude on sites around the ones with a lower onsite energy and has an overlap with the initial product state >0.9, but there are important phases in these new components which ensure that all tunnelling processes to sites at higher energy (although very highly suppressed in the product state) exactly cancel on the Creutz ladder geometry for the eigenstate. In the main text, we set the initial energy offset to μ₀ = 20J.

Note that the correlations of the time-dependently produced state in the main text decay faster than the correlations of the ground state. This discrepancy arrises from the breakdown of the adiabatic principle due to the energy gap closing at the end of the ramp process. However the state that is produced is an eigenstate of the many-body Hamiltonian with a low energy variance, and from exact diagonalisation analysis of smaller systems, we find that the prepared state is the first excited state.