Hubbard models and state preparation in an optical Lieb lattice

Inspired by the growing interest in probing many-body phases in novel two-dimensional lattice geometries we investigate the properties of cold atoms as they could be observed in an optical Lieb lattice. We begin by computing Wannier functions localised at individual sites for a realistic experimental setup, and determining coefficients for a Hubbard-like model. Based on this, we show how experiments could probe the robustness of edge states in a Lieb lattice with diagonal boundary conditions to the effects of interactions and realise strongly correlated many-body phases in this geometry. We then generalise this to interacting particles in a half-filled 1D Lieb ladder, where excitations are dominated by flat band states. We show that for strong attractive interactions, pair correlations are enhanced even when there is strong mixing with the Dirac cone. These findings in 1D raise interesting questions about the phases in the full 2D Lieb lattice which we show can be explored in current experiments.


I. INTRODUCTION
In recent years, there has been much progress in realising and exploring novel lattice geometries using ultracold atoms in optical lattices [1,2], and especially lattices with topological properties [3][4][5][6][7][8][9][10][11][12][13]. While a lot of work has focused on few-particle properties, there is increasing interest in the interplay between topological band structures and strong interactions. For example, recent studies have explored how lattice topology can enhance superfluid phases with cold gases [13][14][15] or potentially superconducting properties in solid-state systems [16][17][18], as a result of the interplay between interactions and flat energy bands.
In this work, we consider the potential experimental realisation of interacting particles in an optical Lieb lattice, as shown in Fig. 1(a). This can be realised in experiments either by using a quantum gas microscope, or following schemes such as that proposed by Di Liberto et al. in Ref. [13] and sketched in Fig. 1(b). We begin by numerically determining localised Wannier functions on individual lattice sites and calculating coefficients for a Hubbard model based on experimentally accessible parameters. We then consider interacting particles in these systems in two example cases of weak and strong interactions. We first probe the robustness of edge and corner states in a Lieb lattice with diagonal boundary conditions [19] when weak interactions are introduced, and then we investigate the many-body phases manifested in a half-filled 1D Lieb ladder for strong attractive interactions. These calculations provide a potential roadmap for experimental investigation of many-body phases for cold atoms in a Lieb lattice.
As depicted in Fig. 1(c), the Lieb lattice has a flat energy band that touches a dispersive band at a single k-point, which becomes a Dirac cone [20] for uniform onsite energies. Dispersionless (or flat) single particle energy bands arise through a destructive interference effect of the single particle wavefunction, [21][22][23][24][25][26][27], leading to an infinite energy degeneracy and means that conventional  [13]. The incoming solid red laser (for VA) is polarised along the y-direction, while the dashed patterned laser (for VB) is polarised in the direction orthogonal to the x-y plane. (c-e) Single particle energy band structure for a tight-binding model on the Lieb lattice for different ratios of the onsite energies, {EA, EB, EC }, for each site of the unit cell (EC = EB for all), given in units of the nearest-neighbour tunnelling rate, J, and the unit cell spacing, a.
single-particle tunnelling is suppressed for particles in superpositions of states belonging to this flat band. It has been shown however, that interactions allow the atoms to form pairs that are mobile [28]. For bosons this can push the ground state of the system towards forming a pair condensate if the lowest energy band is flat [15,[29][30][31][32] and for fermions can result in an enhanced critical temperature for superconductivity [33,34], with phases that have large pair correlations [16,35]. But there are additional open questions here regarding the corrections to previous theoretical predictions, particularly when the interaction strengths are strong enough to mix states in different bands.
In order to analyse the effects of strongly interacting systems in novel band structures it is necessary to derive effective models which can describe the dynamics induced in these systems. While it is now routine to realise Hubbard models in simple periodic potentials, up to well controlled approximations [2,[36][37][38], it can still be a challenging task to derive the coefficients in these models for more complicated lattices. In order to allow future experiments to more easily explore these questions in the Lieb lattice, we first derive experimentally realistic values for the coefficients in a Hubbard model as a function of the laser intensity through calculating the single particle Wannier functions localised on each site [39]. We then show that one could utilise this knowledge in probing the effects of interactions on the edge states that are manifested in a rotated Lieb lattice, where we find that for sufficiently dilute and weakly interacting atoms that they are robust to these effects. We also consider the many-body fermion case at half-filling where we show that pair correlations appear quasi-long range. These strongly dominate over exponentially suppressed single particle correlations and are sufficiently large to be resolved experimentally. Furthermore, we show that the pair correlations become further enhanced when there is greater mixing between states in the flat band and in the higher dispersive bands due to increasing the interaction strength. Finally, in order to be able to realise these phases experimentally we describe and demonstrate an experimentally feasible adiabatic ramp process for preparing these ground states.
The remainder of this article is organised as follows. In section II we compute Wannier functions and determine Hubbard model parameters from potential experimental parameters with cold atoms. In section III, we then discuss geometrically protected edge states, and the effects of interactions on mixing of those states, before discussing strongly interacting systems of fermions at halffilling in a Lieb ladder in section IV. We discuss how these states might be prepared in experiments in section V before giving a summary and outlook in section VI.

II. LOCALISED WANNIER FUNCTIONS AND DERIVATION OF HUBBARD COEFFICIENTS
In the Lieb lattice for the case where every site is degenerate the system manifests a three band touching point, with a Dirac cone that crosses the flat band, but by tuning the ratio of the onsite energies it is possible to realise different energy gaps, see Fig. 1. Here we consider producing the Lieb lattice experimentally by superimposing two optical potentials [13], where the wavevector k i = 2π/λ i , and λ 1 ≈ 2λ 2 . In Fig. 1(a) we illustrate the direction of the beams and the values required for the wave-vectors. Usually, V 1 , V 2 are expressed in units of the recoil energy, E Ri =h 2 k 2 i /2m, where m is the mass of the atoms. The above potential can be produced experimentally using the configuration shown in Fig. 1(b) which follows from the proposal in Ref. [13]. For V A (red) we use laser light that is polarised in the plane of the lattice and for V B (blue) we use a polarisation that is orthogonal to the plane. We note that this geometry has recently been realised by instead superimposing three optical potentials [40,41]. The main differences between this approach and the one that we have considered, is that it requires three optical potentials, instead of only two in this scheme, but has the advantage of not requiring orthogonal polarisations between the beams. Also note that these lattices could alternatively be produced directly through spatial light modulation in a quantum gas microscope [42,43].
By tuning the ratio R = V 2 /V 1 (achieved experimentally by varying the intensities of the two superimposed lasers) around 1, the relative depth of the trap at the A sites and the B, C sites can be varied, which then modifies the single particle energy band structure, see Fig. 1(c-e). Under well controlled approximations it has been shown that atoms confined in an optical lattice can be described with a Hubbard model [1,2,36] for sufficiently large potential depths. We will see below that for the Lieb lattice we require laser intensities of around V 1 , V 2 ∼ 40E R2 so that the next nearest-neighbour tunnelling rates are below 1% of the nearest-neighbour components. In this case the Hubbard model on this geometry for two species fermions takes the form (h ≡ 1) σ,A,n a σ,A,n + E B a † σ,B,n a σ,B,n + E C a † σ,C,n a σ,C,n − J a † σ,A,n a σ,B,n + a † σ,A,n a σ,C,n + a † σ,B,n a σ,A,n− x +a † σ,C,n a σ,A,n+ y + h.c + U Anσ,A,nnσ,A,n + U Bnσ,B,nnσ,B,n + U Cnσ,C,nnσ,C,n , where a † σ,κ,n creates an atom at the unit cell n on site κ andn σ,κ,n = a † σ,κ,n a σ,κ,n . Additionally, x ( y) denotes the translation of one unit cell in the x(y)-direction. For fermions, σ denotes the opposite spin species.
The coefficients in this model can be calculated through the single particle Wannier functions, which can be found by first diagonalising the single particle Hamiltonian, to obtain the Bloch functions, φ m k ( r), associated with each band m. We can then calculate the Wannier functions, where U k n,m is a k-dependent 3 × 3 unitary matrix. For the case where each U k is diagonal and we find the ap- propriate choice for the phase factors then we can obtain unique and real Wannier functions that are exponentially localised to a single unit cell and each represent particles in an individual energy band [44]. This phase is chosen so that the Bloch functions are as smoothly varying as possible around the Brillouin zone. In this way, the Wannier functions for the Lieb lattice can be found as shown in Fig. 2. However, calculating the terms in the Hubbard model using these ordinary Wannier functions would result in a model with many long range tunnelling and interaction processes. It is then advantageous to calculate a new set of basis states which themselves are localised on individual sites, and will allow for the realisation of a Hubbard model of the form of Eq. 2 with only nearest-neighbour processes. In general, this can be a non-trivial task particularly when the lattice geometry has a complicated spatial structure and connectivity. To this end, we use the method presented in Ref. [39] which finds the optimal superposition of the ordinary Wannier functions, by iteratively optimising each U k n,m , such that we obtain a new set of orthogonal functions, but which are much more localised. The effect of the off-diagonal elements of U k n,m is to mix the Bloch functions of each energy band, thus the new basis set corresponds to particles localised on individual lattice sites, but which exist in a superposition of the three energy bands.
While we encountered no technical problems in applying this iterative method in this case, it is possible for the algorithm to find local minima giving rise to a suboptimal Wanner basis. Note that an alternative method for calculating the maximally localised onsite Wannier function has been proposed in Ref. [45] which allows for a non-iterative way of calculating the band mixing matrix U k n,m in Eq. 3, which may lead to advantages for certain systems.
Our calculated maximally localised onsite basis is most convenient for deriving the coefficients of a Hubbard model because we can use them to obtain the tunnelling amplitudes between nearest-neighbouring sites and the onsite energy and interaction coefficients as a function of the optical lattice laser potential values which is important for experimentally realising theoretical predictions that require knowledge of the interaction strength and are paramount for being able to reliably explore the physics of many-body systems in complex geometries. These coefficients are shown in Fig. 3 for the regimes where the tight-binding spectrum is well approximated, where we can see that the onsite interaction strengths can be tuned to values close to the tunnelling rates for moderate scattering lengths. Comparing these maximally localised Wannier functions to those calculated in Ref. [40] for the alternative experimental realisation, we find very similar results. In particular, we also find that the Wannier states centred on the B/C sites to be more localised than those centred on the A sites, which we believe will result in the same asymmetry in the onsite interaction values found in this description.
Note that compared to a conventional square lattice, which requires intensities of V 0 > 5E R in order to accurately approximate a tight-binding model, we find here that we must use a larger laser intensity, V 0 > 25E R in order to reproduce a spectrum similar to the tight-binding case in Fig. 1(c-e). However, as can be seen from Fig. 2(a) the potential is greatly dominated by the large peaks which creates the distinct pattern of the Lieb lattice and the potential barriers between nearest-neighbouring sites are much lower, giving values for the nearest-neighbour tunnelling rates on the order of J ∼ 100 Hz for all parameters considered in Fig. 3. This means that while we require a higher laser power to produce the tight-binding spectrum compared to conventional lattices, the dynamics manifested in this system will be at a similar rate.
We are also able to calculate corrections to a standard Hubbard model, such as longer range off site processes. In Fig. 3(b) we plot the next-nearest-neighbour single particle tunnelling rate and the largest nearest-neighbour density assisted tunnelling coefficient, which is explicitly given by where g = 4πh 2 a s /2m, with a s the characteristic scattering length and w A ( r), w B ( r) are the Wannier functions centred on an A site and B site respectively. We can see that these values for all of these correction terms are on the order of 1% for the potential values illustrated here. Additionally, the increasing behaviour of the density assisted interaction term as the potential strength is increased arises from the non-trivial interplay between FIG. 3. Local Hubbard coefficients produced from Eq. 1 as a function of the laser intensity, calculated with a maximally localised Wannier basis for each site. In the red zone the flat band is degenerate with the higher dispersive band (see Fig. 1(d)), in the blue zone it is degenerate with the lower energy band (see Fig. 1(e)) and at the boundary the band structure manifests the three band touching point and Dirac cone (see Fig. 1(c)). the two optical potentials. As such, increasing the potential barrier height has a great effect on the large potential peaks producing the Lieb lattice pattern, which increases the confinement of the B/C sites, in the directions orthogonal to the tunnelling between these sites and the nearest A sites (see the elliptical shape of these sites in Fig 2(a)). This means that in Eq. 4 the weight of the onsite terms |w B ( r)| 2 w B ( r) increases and will lead to an increase in U tun . In a conventional lattice, this effect is usually counteracted by a much smaller overlap between nearest-neighbour Wannier functions which then results in these tunnelling rates decreasing as the potential is increased, however for the case considered here, the barrier height between the A and B sites, is only weakly affected by an increase in the potential and so the overlap between w A ( r) and w B ( r) does not decrease rapidly enough to counteract the effects of the increased confinement of the B sites. Thus resulting in a density assisted tunnelling rate that slowly grows with increasing potential depth for the regime considered here.
The most dominant effect of the next nearestneighbour tunnelling term is to introduce a small curvature and bandwidth for the flat energy band, however, this is also much smaller than the dominant nearestneighbour tunnelling rates (< 1%) indicating that these corrections will have a negligible effect on the dynamics as they are significantly dominated by the conventional Hubbard terms.

III. GEOMETRICALLY PROTECTED EDGE STATES
A first application of studying interactions in this setup is to understand how single-particle properties are affected by interactions, and to what extent these are modified in experiments with unavoidable weak interactions. A good example of interesting single-particle phenomena for the Lieb lattice occurs when we consider a rotated Lieb lattice with a boundary as shown in Fig. 4(a), which would require additional optical potentials to act as a hard boundary. This model has edge states for all values of dimerisation J 1 /J 2 and corner states for the C 4 symmetric regime at J 1 = J 2 . In this system, the corner states are similar to edge states in the sense they are localised at the systems boundary and have an exponentially decaying amplitude in the bulk towards the centre of the lattice. The main differences are that the corner states also have exponential suppression for amplitudes towards the middle of the edges and so are localised to the corners. Additionally, they do not exist in an energy gap, and are instead degenerate with the bulk bands. The subset of edge and corner states that we will consider here, are those that are produced through a geometrical frustration effect similar to the mechanism that gives rise to the bulk flat band states, where as the particle attempts to move away from the edge (or corner) the various tunnelling components from the eigenstate onto the nearest bulk state interfere destructively, suppressing transport and creating localised states. We illustrated in Ref. [19] by considering periodic boundary conditions in the y-direction, that these geometrically enhanced edge states exist at the ak y = π point in the Brillouin zone.
Also, in Ref. [19] we showed that for non-interacting particles, these edge states can be prepared by an adiabatic ramping protocol, where we simulated a timedependent variation of the value of the tunnelling dimerisation J 1 /J 2 . In Fig. 4(b) we illustrate the energy spectrum throughout this time-dependent process, where we begin in an edge state at J 1 /J 2 = 2, with similar spatial structure to the one shown in blue in Fig. 4(a) and linearly ramp this ratio to J 1 = J 2 in a set ramp time. In Ref. [19] we found that even though the time-dependent state crosses through many bulk energy states (blue lines) that there were no avoided crossings and we were able to prepare the final corner state, indicating that the geometrical frustration mechanism suppresses transport away from the edge and protects mixing with the bulk states even in the presence of time-dependent perturbations.
Here, we can consider how interactions affect this protocol, and how small the interactions should be so as to ensure that the final state can be reached with highfidelity in this ramping protocol. We employ the same ramp process, but now with two (bosonic and same species) atoms in the initial edge state that can interact with an onsite interaction strength, U . In Fig. 4(c) we plot the final projection of our produced state onto the corner eigenstates of the final Hamiltonian where U = 0 for a range of ramp times and onsite interaction strengths. As expected it is apparent that as we increase the interaction strength that the final fidelity decreases, indicating that the interactions can cause mixing between the edge states and the bulk states somewhat breaking the geometrical frustration. It is also clear that beyond weak interactions U > 0.1J 2 , increasing the ramp time beyond some value dependent on U leads to lower fidelities. This is because we calculate the fidelities with respect to the U = 0 eigenstates, so this feature indicates that the eigenstates for U = 0 have a decreasing overlap with the non-interacting states as U is increased. In some cases, this can then result in a better preparation of the U = 0 corner states for non-adiabatic processes at shorter ramp times. We further probe this in Fig. 4(d) where we calculate the time-dependent projection throughout a single ramp (J 2 T = 100) of the current state onto the product state where both particles are in the single particle eigenstate highlighted on the energy spectrum in red in Fig. 4(b). Comparing different strengths of interactions we illustrate that for weak interactions (U < 0.2J), the majority of the mixing occurs when the edge state crosses the bulk states for J 1 /J 2 < √ 2, whereas for stronger interaction strengths (U ∼ 0.4J) there is also large mixing with the band of higher energy states at short times when J 1 /J 2 > √ 2. However, even for moderate interaction strengths U = 0.2J 2 it is still possible to achieve fidelities F ∼ 0.9 for experimentally feasible ramp times, J 2 T ∼ 100. This indicates that even with experimental errors in the tuning of the interaction strength, it is possible to prepare these states for sufficiently dilute and weakly interacting samples.

IV. MANY-BODY CORRELATIONS AT HALF-FILLING
There have been several previous studies on twodimensional systems with flat bands, for example [25][26][27]46], which utilise either mean-field theory, quantum Monte-Carlo, or dynamical mean-field theory techniques. Nevertheless, there are still open questions here on the validity of a conventional BCS treatment in a flat band, as the interaction strength dominates the dynamics resulting in strongly correlated phases even for weak interactions.
A full quantum treatment for fermions in onedimensional flat band systems has also been explored [47,48] with matrix product states (MPS) [49], but usually only in the isolated flat band approximation where it has been shown that the ground state can coincide with the BCS ground state [48]. There are also additional questions here as to the effects on these properties in systems that have a flat band as well as dispersive bands, such as the Lieb lattice, and in particular when the interaction strengths are strong enough to mix states in different bands. However, current state of the art classical methods for quantitatively analysing these phases, such as MPS are only practical for one-dimensional systems, where it is known that Bose-Einstein condensation or superconductivity cannot occur in the ground state and instead at best the system can develop quasi-longe range order in the correlations [50]. While there have been some recent numerical analysis on a minimal quasi-2D model for superconductivity [51], this is an interesting opportunity to potentially use quantum simulators to improve the understanding of strongly correlated phases in 2D.
In order to perform a quantitative analysis of the many-body phases produced in the Lieb lattice we focus on a one-dimensional ladder cut, shown in Fig. 5, and calculate the ground state using MPS techniques. This allows us to quantitatively calculate the properties of the ground state in the strongly interacting regimes while also preserving the main qualitative features of the full Lieb lattice and providing a way to benchmark a quantum simulation of the full 2D geometry. For this ladder system we have a similar energy band structure with a flat band and a Dirac cone (see Fig. 5). In the fol- lowing we consider the case of fermions at half-filling as in this regime the lower two energy bands are filled and the flat band is half filled, meaning that all excitations and dynamics around the Fermi level will be dominated by states in the flat band. Here we consider the case of uniform interaction strength on all sites of the unit cell, U = U A = U B,C , for simplicity, but note that we only expect small quantitative corrections to the correlation functions below as the low energy properties are dominated by the flat band Wannier basis, which have no amplitude on the A sites. This means that including an interaction imbalance on the A site should not strongly affect the behaviour we find below. Of course, for increasing interaction strength, there will be strong mixing with Wannier states in higher energy bands which do have some amplitude on the A sites, meaning that these corrections may grow as the interaction strength is increased. We calculate the ground state for attractive interactions as the interaction strength is increased for the two band structures illustrated in Fig. 5. We calculate the single particle correlations, where we first apply the creation operator on the A site in the centre of the system and then apply the destruction operator at every other site to the right in the system.
We compare these to the Cooper pair correlations where ∆ r = a ↑,r a ↓,r , and we first apply the pair creation operator to the B site in the centre of the system. We plot the results in Fig. 5, where we can see that for non-zero interactions the pair correlations dominate over the single particle correlations for both regimes. Due to the spatial distribution of the flat band eigenstates (shown in Fig. 5(a)) we find that the dominant pair correlations occur between B and C sites, where the noticeable dips in the plots for lower interaction strengths occur when the subsequent destruction operator is applied to an A site. As the interaction strength is increased these features are washed out, indicating that there is strong mixing with Wannier states associated with the higher energy bands, but intriguingly, the overall decay of the correlations seems to actually be enhanced by this mixing, whereas the single particle correlations become further suppressed.
We can also see that in the regime E A = E B = E C with the Dirac cone (a-b), that the pair correlations are suppressed for weak interactions, but for strong interactions the pair correlations seem to become even larger than the case with a gap in the band structure (E A = J) and for U = −4J appear to remain almost constant over the entire length of the system that we have considered. Of course, with these long range correlations, finite size effects become important, making it difficult to really unequivocally determine if these correlations decay algebraically, or simply have a correlation length that is larger than the system. Nevertheless, these results are very suggestive of quasi-long range superconducting pairing, and also demonstrate that strongly mixing with dispersive energy bands, rather than destroying the novel pair dominated phases in flat bands, can actually enhance these properties.
A further interesting question is how these features will change upon including a magnetic flux through the loops in each plaquette of the ladder, which we include in the Hamiltonian with a Peierls substitution [52]. Including these features in an experimental implementation requires modifying the above proposal, but can be achieved using Floquet engineering or laser induced tunnelling for example [53]. In the Lieb lattice, this has been shown to introduce energy gaps between the flat band and the higher and lower dispersive bands [54] which we explicitly show in Fig. 6(a). In Fig. 6(b) we calculate the single particle Wannier functions (using Eq. 3 on an infinite tightbinding model) where we illustrate that including a magnetic flux qualitatively changes the long-range behaviour of this basis where compared to an exponential decay in the φ = 0 limit they develop algebraically decaying tails with degree ∼ 2 which in higher dimensions would suggest a topological Chern insulator [55,56]. Additionally, we can see that at short distances these basis states become strongly extended to several nearest-neighbouring unit cells [52,57]. We investigate how these features affect the correlations upon including strong onsite interactions by performing a similar analysis as before and we compare the single particle and pair correlations in Fig. 6(c-f) for a range of values for the magnetic flux, φ. We compare the case for U = −0.5J, which is an interaction strength that is smaller than the maximum energy gap at φ = π, to U = −4J which is sufficiently strong to allow mixing between single particle states even at this maximum gap opening point. From this analysis, we can see that the overall trend in the single particle correlations (c,e) remains unaffected and these retain their exponentially decaying behaviour, indicating that the many-body effects dominate over the algebraic decay of the Wannier states even for weak interactions. However, the pair correlations also become exponentially suppressed as φ is increased, reaching a maximum suppression at φ ∼ π/2. Further increasing the flux, the correlations begin to increase again until the maximum gap opening point at φ = π, where there is a similar algebraic like decay compared to the φ = 0 case, but now with a distinctly different spatial behaviour, reflecting the change in the short distance behaviour of the single particle basis states.

V. STATE PREPARATION
In this section, we now consider how these many-body phases might be prepared in an experimental realisation of the Lieb lattice with an optical lattice. This preparation scheme is based on adiabatic state preparation [58,59] and utilises the ability of the experimental systems to time-dependently vary the onsite energies, interaction strengths and tunnelling amplitudes. In order to achieve this experimentally, it requires precise knowledge of the relationship between the parameters of the optical potential and the coefficients in the effective Hubbard model which can be understood with the help of the analysis presented in the previous sections.
We require the preparation of an initial product state where a single fermion is projected onto each of the blue sites shown in Fig. 7(a). Our initial state in this process is the case where all tunnelling rates are zero which can be achieved with strong laser intensity creating large potential barriers and where the onsite interaction strengths are also zero, which can be tuned with a magnetic field around a Feshbach resonance [60]. Additionally, we require a large detuning between onsite energies with much lower energies for the sites that are initially populated, which can be achieved by varying the relative intensity of the two optical potentials creating our lattice in Eq. 1. We then simultaneously time dependently ramp the tunnelling, the onsite energy and the onsite interaction strength to the desired value.
First in Fig. 7(b) we consider the case for zero interactions, with E A = J; E B = E C = 0, where we linearly ramp the parameters in the Hubbard model. For the non-interacting case the ground state is highly degenerate so we plot the projection of the final state onto the ground state manifold for small system sizes using exact diagonalisation. We can see that for ramp times that are achievable before decoherence effects begin to lead to errors (T J ∼ 200) we can achieve fidelities > 1 − 1 × 10 −3 for the system sizes considered here.
Next we investigate the more interesting case of interacting ground states. We begin with the same initial state, but we now apply an exponential ramp for the parameters, from t = 0 to t = T . As the energy gap between the ground state and first excited state decreases as we ramp these parameters, with the final state being gapless in the thermodynamic limit we find it optimal to apply this exponential ramp which ensures that there is a slower change in the parameters as the gap decreases. This ensures that the adiabatic principle is not as strongly violated compared to the case for a linear ramp. In Fig. 7(c) we plot the projection of the final state onto the ground state using exact diagonalisation on a small system for both band structures with E A = J and E A = 0. We consider a range of attractive interactions where we can see that as we increase the interaction strength, the projection on the ground state improves. It is also apparent that we can achieve near perfect fidelities for experimentally achievable ramp times (T J ∼ 200) for strong interaction strengths, which are the phases with the most dominant pair correlations, see Fig. 5.
Of course, as this analysis was carried out on small system sizes, we are only able to prepare these states so accurately due to finite size effects introducing effective energy gaps in the spectra. So, we also simulate the same ramp processes on larger system sizes using time-evolution techniques with MPS [61,62], where in Fig. 8 we compare the correlations produced through the time-dependent preparation to those for the ground state, again for both band structures, E A = J (a) and We can see that as we increase the total ramp time, the correlations of the prepared state quickly approach those of the ground state, where for experimentally feasible ramp times (T J ∼ 100) the deviations become too small to be resolved experimentally. Note that in all cases, the single particle correlations are exponentially suppressed. This analysis indicates that while preparing these ground states with near perfect fidelities may not be possible using adiabatic preparation schemes, due to the fact that the states are gapless in the thermodynamic limit, it is clear that we can feasibly prepare states with strong pair correlations that have a similar dependence to those of the true ground state.
This analysis touches upon an important technical consideration in these experimental architectures, namely reducing the effective temperature of the ultra-cold fermionic gas. As these atomic systems are incredibly well isolated, such that on typical experimental timescales (T J < 500) they can be considered closed systems, an equivalent consideration is in reducing the entropy of the many-body system [42,63,64]. The simulations presented in Fig. 8 indicate that, if the initial product state can be accurately prepared in the experiment (assumed to have zero entropy in the simulation) then the proposed experimental sequence can in principle be prepared with low enough entropy (and therefore effective temperature) in timescales (T J ∼ 200) where additional heating effects can be neglected, such as spontaneous emission [65].
Of course there are still open questions here regarding the critical temperatures required for these pair correlations to dominate. As for typical cases with attractive pairs in optical lattices, in order to realise the proposed many-body ground states we require the temperature to be small compared to the effective dispersion relation for these pairs. In Ref. [15] we showed that a flat band system can host stable bound pairs with a significantly enhanced pair kinetic energy which grows with interaction strength, suggesting that the critical temperature required to realise pair dominated phases may be larger compared to conventional systems and may even also grow with increasing interactions. These are important and relevant considerations which will be interesting to explore in the future.

VI. CONCLUSION
We have considered an experimental realisation of the Lieb lattice, determining realistic parameters in a Hubbard model, and showing how to prepare states with weak interaction, and also phases where strong interactions give rise to enhanced pairing correlations. By performing numerical calculations on a one-dimensional ladder system we quantitatively demonstrated that we can prepare the half filled fermion ground state in experimentally realistic timescales. Additionally we provided predictions for the ground state correlations for attractive interactions where we found that due to the flat energy band the pair correlations dominate and can be further enhanced with a timestep Jτ = 0.1 and D = 512. Note that the lattice sites in these plots are ordered according to the site index numbers in Fig. 5(a) by increasing the interaction strength such that there is strong mixing with the higher dispersive energy bands. We also found that in the strongly interacting regime, mixing between the flat band states and states close to the Dirac cone resulted in correlations that were almost constant throughout the full system that we considered. These calculations provide a starting point for realisation of these models using parameters available in current experiments. On the theoretical side this work has further implications for understanding strongly interacting phases in flat band systems and how mixing with additional dispersive bands affects the properties of these phases. Additionally, our experimental proposal demonstrates that it is, in principle, possible to realise these phases with strong pairing correlations in current cold atom experiments.