Topological repulsively bound atom pairs in optical lattice ladders

There is a growing interest in using cold-atom systems to explore the effects of strong interactions in topological band structures. Here we investigate interacting bosons in a Cruetz ladder, which is characterised by topological flat energy bands where it has been proposed that interactions can lead to the formation of bound atomic pairs giving rise to pair superfluidity. By investigating realistic experimental implementations, we understand how the relatively large effective pair-tunnelling in these systems can lead to robust pair superfluidity, and we find lattice supersolid phases involving only pairs. We identify schemes for preparation of these phases via time-dependent parameter variation and look at ways to detect and characterise these systems in a lattice. This work provides a starting point for investigating the interplay between the effects of topology, interactions and pairing in more general systems, with potential future connections to quantum simulation of topological materials.

There is a growing interest in using cold-atom systems to explore the effects of strong interactions in topological band structures. Here we investigate interacting bosons in a Cruetz ladder, which is characterised by topological flat energy bands where it has been proposed that interactions can lead to the formation of bound atomic pairs giving rise to pair superfluidity. By investigating realistic experimental implementations, we understand how the relatively large effective pair-tunnelling in these systems can lead to robust pair superfluidity, and we find lattice supersolid phases involving only pairs. We identify schemes for preparation of these phases via time-dependent parameter variation and look at ways to detect and characterise these systems in a lattice. This work provides a starting point for investigating the interplay between the effects of topology, interactions and pairing in more general systems, with potential future connections to quantum simulation of topological materials.
Introduction. Recent experiments have demonstrated the utility of ultra-cold atoms in optical lattices to explore the physics of topological quantum systems [1][2][3][4][5][6][7][8][9][10][11][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][16][17][18][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][21][22][23][24][25][26][27][28][29][30]. In this work, we investigate interacting bosons in a topological band structure where the single-particle kinetic energy is completely frustrated [31][32][33][34][35][36][37][38][39][40][41][42][43][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 Fig. 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][36][37][38][39][40][46][47][48][49][50][51][52][53]. There is growing interest in repulsively interacting bound pairs in general cold-atom systems [54,55], 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.
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]. Furthermore, for high densities we identify new lattice supersolid phases [56][57][58][59][60][61][62][63][64][65][66], 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 [67,68], 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.

Model.
In Fig. 1 (a) we include the Creutz ladder which is characterised by complex tunnelling amplitudes that, through a destructive interference effect, gives rise to a dispersion relation with only flat energy bands (see Fig. 1 (b)). There are various ways to produce this experimentally [45,69,70], where most proposals utilise a synthetic dimension [71][72][73][74][75][76][77][78][79] where each leg of the ladder corresponds to two atomic internal states. However including interactions will then require care- ful 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, 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 Supplementary Material). Explicitly the Hamiltonian including an onsite twobody contact interaction, U A for the A-sites and U B for the B-sites, is given by ( = 1), whereâ † n (b † n ) 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 Fig. 1 (a) and the transformation is, where theŴ ± n 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, whereŴ n =Ŵ + n +Ŵ − n and W n =Ŵ + n −Ŵ − 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, W † nŴ † n W n+1 W n+1 , as well as nearest-neighbour interactions,Ŵ † nŴn W † n+1 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.
Topological bound pairs. 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 Supplementary Material). In Fig. 1 (c) we plot the resulting dispersion of the lowest bound state band (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 Fig. 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. The asymmetry in Fig. 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 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 Fig. 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 2 /U [55]. 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 topological pairs in the Creutz ladder have a kinetic energy that is larger than conventional pairs by nearly an order of magnitude. 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 topological 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. Phase diagram.
We now consider the many-body bosonic case and characterise the phases that are manifested by the topological 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,80]. 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 [81], note that we increase the local dimension such that all bosonic fluctuations are captured. We include a chemical potential term, n , in Eq. 3 and calculate the pair correlations in the site basis, â † rb † râ0b0 , 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 râ0b0 | (solid lines) and the single-particle correlations, Φ = | â † râ0 | (crosses) for an imbalance between the onsite interactions, UA = UB, in the ρ < 1 pair Luttinger liquid (PLL) phase, µ/UA = 0.1. lack of single-particle dispersion in this system. 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 CDW phase respectively, examples of which we have included in Fig. 3(a-b). The phase diagram is shown in Fig. 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 [82,83], which is given by, We agree with the predictions of Ref. [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 Fig. 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 singleparticle 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.
Experimental preparation and detection.
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 [67,68] (see Supplementary Material 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 Fig. 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 Fig. 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, T J, is varied. It is clear that we can produce a many-body state with significant pair correlations, â † rb † râ 0b0 and vanishing single-particle correlations, b † r b 0 and a † r a 0 , in experimentally feasibly timescales (T J = 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 (T J = 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.
We also consider the effects of a finite temperature on the pair correlations in the system and we use an imagi- nary time MPS algorithm which utilises a purification of the density matrix to calculate the state at a given temperature [84,85]. In Fig. 5 (b) we plot the pair correlations at varying temperature and find that the correlations are exponentially suppressed at high temperatures (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. Fig. 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 <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 temperaturesT < 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 [86,87].
Explicitly we calculate the unequal time two-point correlator, A(r, t) = Ψ 0 |δn r (t)δn 0 (0)|Ψ 0 , where |Ψ 0 is the initial state, δn r (t) =Ŵ − † r (t)Ŵ − r (t) − ρ where ρ is the density. We then take the Fourier transform which we plot in Fig. 6 for the different phases. Note that in all cases, U A = U B .
We can see for the PLL phase (a) there is a dominant linear excitation originating from the k = π mode matching the predictions from Luttinger liquid theory (see Supplementary Material) 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 (Fig. 6(c-d)) we can see gapless excitations that are linear for low energies, 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 Fig. 6(c-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 Fig. 6(d) there is a local minimum in the low energy dispersion which resembles the helium roton spectrum [88][89][90], 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.
Conclusion. 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].
Acknowledgements. This work was supported in part by the EPSRC Programme Grant DesOEQ (EP/P009565/1), and by the European Unions Horizon 2020 research and innovation program under grant agreement No. 817482 PASQuanS. S.F. 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 Knowledge-Base at https://doi.org/10.15129/81ecb9ec-7701-4940-b496-9f2785079198.

Appendix A: Experimental Realisation
Experimental realisations of the Creutz ladder have been proposed in Ref. [45,69], 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 singleparticle 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 Fig. 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 Fig. 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 Fig. 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 Fig. 7(a) is able to quantitatively reproduce all features presented in the main text, while also offering a more viable experimental implementation.
For this scheme we require a particular separation of the onsite energy levels, shown in Fig. 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 (Fig. 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, [91] where φ * ( r) are the onsite Wannier functions, δ k = k 1 − k 2 is the difference between the wave vectors of the two lasers in a single Raman process, | k i | = 2πa/λ i , δ is the detuning between ω 1 , ω 2 and the excited internal state (see Fig. 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 Fig. 7(c) is, where all terms correspond to tunnelling events from left to right in the lattice, along the directions of the arrows in Fig. 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, |ω 1 − ω 2 | ≪ ω 1 , ω 2 , where ω 1 > ω 2 then we can assume that the wavelengths for each component has the same magnitude, λ ≡ λ 1 ≈ λ 2 , 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) where on the right we have shown the values that would produce the Creutz ladder shown in Fig. 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  8. Possible angles and wavelengths for the two fields that are responsible for producing the tunnelling processes in the dimerised Creutz ladder shown in Fig. 7 (c). Wavelengths given in units of the lattice spacing a.
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, This requirement, means that motion in the x-direction does not give rise to any phase change, as cos θ 2 − cos θ 1 = 0. And so all phase dependence comes from moving in the y-direction with sin θ 2 − sin θ 1 = −2 sin θ 1 , 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 θ 1 and θ 2 must be tuned with the wavelength, λ in order to obtain the correct phases. As an example, we plot one possible choice in Fig. 7(c) where θ 1 = −π/2 and θ 2 = π/2, which would then require a wavelength, λ = 4a. And in Fig. 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 Eq. 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.
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.

Appendix B: Bound State Model
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 Eq. 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, and the second set corresponds to the two atoms on neighbouring unit cells, 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 Eqs. B1 & B2 form an orthonormal basis set and we can solve the system exactly. Explicitly, the momentum space Hamiltonian is, where and In Fig.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.
Appendix C: Universal behaviour in the pair superfluid phases As we have a 1D superfluid, we expect to be able to describe the superfluid phases by mapping to a homogeneous Luttinger liquid model [92], ( = 1) where the bosonic field, φ(r), and its conjugate momentum density, Π(r), satisfy the commutation relation, [Π(r), φ(r ′ )] = iδ(r − r ′ ). 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, 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, where we have evaluated the dρ/dµ numerically from our data presented in Fig. 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 Ref. [93] 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 Ref. [94] 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 [95,96], through, Below, we perform this analysis for both pair superfluid phases in density ranges ρ < 1 and then ρ > 1.
We begin with the ρ < 1 superfluid phase and apply the Luttinger liquid formalism. Ref. [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 [92,97]. 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 Fig. 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.
We also plot the inverse effective mass, 1/m * in Fig. 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, 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 Fig. 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. A similar analysis was carried out in Ref. [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. the case for hard-core bosons [92], indicating that we can realise a regime where the physical onsite interactions are weak but the effective onsite interactions are infinite. The values for the inverse effective mass approach the same magnitude as the u values for smaller chemical potential values, indicating that there may be additional features present that are beyond Luttinger liquid theory. In particular we can see from Fig. 6(a) in the main text that the excitation spectrum in this phase is dominated by a linear branch, which is predicted from Luttinger liquid theory, but there is also clearly a sinusoidal branch present further indicating that there are additional features present here.
We also calculate the density-density structure factors and find that it is peaked at k = π for all µ/U . In Fig. 9 (d) we plot the fraction of the population that is in this peak, which indicates that the system is strongly condensed for larger µ/U and when K ≈ 1 the condensate peak is somewhat suppressed by the strong effective interactions.

Appendix D: Phase Separation
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 = E(µ 1 )/2 + E(µ 2 )/2, with E(µ 0 ) such that the density, ρ(µ 1 )/2 + ρ(µ 2 )/2 = ρ(µ 0 ). Where we choose µ 0 to be in the ρ < 2 PSS phase, µ 1 , to be in the ρ > 1 PLL phase and µ 2 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.

Appendix E: Adiabatic Preparation Scheme
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 [67,68], 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 0 /J = −µ 0 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, 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 T J = 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 µ 0 = 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.