# Staggered Ground States in an Optical Lattice

###### Abstract

Non-standard Bose-Hubbard models can exhibit rich ground state phase diagrams, even when considering the one-dimensional limit. By utilising a self-consistent Gutzwiller diagonalisation approach, we study the mean-field ground state properties of a long-range interacting atomic gas in a one-dimensional optical lattice. We first confirm that the inclusion of long-range two-body interactions to the standard Bose-Hubbard model introduces density wave and supersolid phases. However, the introduction of pair and density-dependent tunnelling can result in new phases with two-site periodic density, single-particle transport and two-body transport order parameters. These staggered phases are potentially a mean-field signature of the known novel twisted superfluids found via a DMRG approach [PRA 94, 011603(R) (2016)]. We also observe other unconventional phases, which are characterised by sign staggered order parameters between adjacent lattice sites.

## I Introduction

The Bose-Hubbard model and its extensions has long been a theoretical workhorse for lattice based systems Gutzwiller (1963); Kanamori (1963); Hubbard (1963, 1964), including in the field of ultracold atoms in optical lattices Jaksch et al. (1998); Bloch (2005); Lewenstein et al. (2007); Dutta et al. (2015). One of the first topics of great interest to the ultracold gas community was the Mott-insulator to superfluid transition in the standard Bose-Hubbard model Fisher et al. (1989); Freericks and Monien (1994); Greiner et al. (2002); Zwerger (2003); Bakr et al. (2010). Ultracold gases in optical lattices allows for full control of the underlying periodicity and inherently produces no defects. In addition, several techniques also exist that can tune the strength of interaction and correlation processes such as Feshbach resonances Fedichev et al. (1996); Inouye et al. (1998); Theis et al. (2004); Dickerscheid et al. (2005) and laser assisted tunnelling Miyake et al. (2013); Wozniak et al. (2018). The standard Bose-Hubbard model is known to be a poor approximation for systems with strong, non-trivial long-range interactions, as can be realised in dipolar atomic species Menotti et al. (2008); Lucioni et al. (2018); Griesmaier (2007). Such additions to the standard Bose-Hubbard model are referred to as extended, or non-standard, Bose-Hubbard models, where long-range phenomena can significantly change the preoperties of the system Kremer et al. (2017). An example of a long-range interacting atomic gas is the case of dipolar atoms Aikawa et al. (2012); Lu et al. (2012), for which a long-range dipole-dipole interaction is present that decays as Trefzger et al. (2011).

For long-range interactions, non-standard terms can include a density-dependent tunnnelling, pair tunnelling, and/or inter-site interactions. It is known that introducing a density-dependent tunnelling changes the critical point of the Mott-insulator to superfluid transition Maik et al. (2013). In addition to changing ground state properties, this term also affects the dynamics of the system Łącki and Zakrzewski (2013). For dipolar interactions the additional non-standard pair tunnelling can even destroy the Mott-insulating domains and introduce new phases Sowiński et al. (2012a); Maik et al. (2013), including the pair superfluid. It is also known that the introduction of nearest-neighbour two-body interactions can induce density wave and supersolid ground states, which spontaneously break the translational symmetry of the lattice Sowiński et al. (2012b); Mosadeq and Asgari (2015); Rossini and Fazio (2012); Kawaki et al. (2017); Ohgoe et al. (2012).

In this work, we will consider the ground state phases of atoms with long-range interactions in one-dimensional optical lattices in detail by a Gutzwiller mean-field approach. There are a significant number of works in the current literature which consider the constrained density-dependent or extended Bose-Hubbard model Rossini and Fazio (2012); Kawaki et al. (2017); Ohgoe et al. (2012); Trefzger et al. (2011); Biedroń et al. (2018). By including all terms, the resulting phase diagrams can differ significantly even for modest parameter strengths. We will confirm the destruction of the Mott-insulating phase and the introduction of known supersolid, density wave and pair superfluid phases Jürgensen et al. (2015). By performing a detailed study of the ground state phase diagrams for various parameter regions, we find new staggered superfluid and supersolid phases, including sign staggered behaviour of the ordinary and pair superfluid and supersolid.

We will begin by defining the Bose-Hubbard model for a long-range interacting atomic gas, and discuss the derivation of the parameter strengths in terms of the well-known Wannier functions. In Sec. III, we will discuss the mean-field approach utilised in this work. The various phases encountered in the course of this work are defined and discussed in Sec. IV. We will then consider the phase diagrams for each non-standard Bose-Hubbard term being non-zero and the physically relevant case of combinations of non-zero additional terms in Sec. V.

## Ii Bose-Hubbard model

In this section we will define the one-dimensional Bose-Hubbard model for atoms with long-range interactions. The many-body Hamiltonian in second quantised form of an ultracold gas in an optical lattice described by is given by

(1) | |||

where is the chemical potential, are the bosonic annihilation (creation) field operators obeying the standard canonical commutation relations, is the two-body interaction potential, and is the mass.

We consider the case of long-range interacting atoms (e.g. dipolar atoms), which have an interaction potential of the form

(2) |

with a short-range interaction, usually of contact-type, of

(3) |

and a long-range interaction of the form

(4) |

where and are dimensionless scaling pre-factors and is the non-local spatial profile of the interaction. As an example, for dipolar atoms, the non-local spatial profile is of form .

For a periodic external potential (e.g. an optical lattice), the continuous field operators may be described by discrete mode excitations by virtue of Bloch’s theorem,

(5) |

where labels the band, is the quasi-momentum, are the bosonic particle destruction (creation) operators, and characterises the wave function. In the tight-binding limit the delocalised can be expressed in terms of localised and orthogonal Wannier functions, i.e.

(6) |

with corresponding to the lattice translation vector and being the Wannier function in the th band. Substituting the form in terms of Wannier functions into Eq. (5) gives

(7) |

The overlap and extent of the Wannier functions is determined by the external potential’s depth and is independent of the interaction strength and mechanism.

The usual tight-binding limit considers the depth of the lattice potential to be sufficient enough that the atoms are well localised to each lattice site, in an analogous way to electrons being tightly bound to atoms in solid state crystals Rosenberg (1988). In this limit, we can expand the wave function into the lowest set of Wannier functions. As the interaction and tunnelling terms usually decay substantially as the distance between sites increases, it is usually a good approximation to constrain the terms of the Hamiltonian to on-site and nearest-neighbour only. Using the Wannier representation, the general Bose-Hubbard model for long-range interacting atomic gases can be derived from Eq. (1) as,

(8) | ||||

where are labels of the lattice sites, is the number operator, and indicates nearest-neighbour summations.

In Hamiltonian (8) there are three terms which are contained in the standard Bose-Hubbard model; tunnelling with strength , two-body on-site interactions of strength , and a chemical potential of . However, there are more exotic terms which describe other possible two-body processes. The first of these terms is the two-body nearest-neighbour interaction of strength , which is vital when considering dipolar BECs. There is also a term which denotes the pair tunnelling process, i.e. when two atoms tunnel to a site together, which is of strength . The final term is the density-dependent tunnelling, which is a single-particle tunnelling process that depends on the density of atoms at each site involved in the tunnelling process. We will denote the strength of the density-dependent term as . Each of the five dynamic tunnelling and interaction terms in Hamiltonian (8) are illustrated in Fig. 1.

Each of the coefficients for the terms of Hamiltonian (8) are defined in terms of overlap integrals of the Wannier functions. The single-particle tunnelling strength is given by

(9) |

and the two-body on-site interaction by

(10) |

where

(11) |

stands for the short-range interaction integral and

(12) |

for the long-range interaction integral and are general labellings, with each being either or in the overlap integrals. We also get here the form of all three non-standard Bose-Hubbard terms, with the long-range two-body interaction given by

(13) |

the density-dependent single-particle tunnelling by

(14) |

and the two-body pair tunnelling by

(15) |

We consider on-site and nearest-neighbour terms only, as all higher order terms will be small. This means we consider only the cases of . All three non-standard Bose-Hubbard terms depend on the long-range interaction, and can, therefore, be important for the case of strong dipolar interactions. We note that the strength of the two-body terms follow the relation

(16) |

All parameters in Hamiltonian (8) are dependent on combinations of the external and two-body interaction potentials, as well as an overlap integral of the Wannier functions. However, in order to understand the effects of each process individually in this work, we will consider the parameters of Hamiltonian (8) to be independent variables. This serves to show the effect of each term on the ground state.

## Iii Self-consistent Gutzwiller mean-field

In order to study the ground state phases that are possible within the general Bose-Hubbard model, we will consider a Gutzwiller mean-field approach. For interacting problems, the task of exact diagonalisation becomes unfeasible due to the exponential increase in the dimension of the Hilbert space. The Gutzwiller approach was first developed in the 60’s for fermions Gutzwiller (1963, 1964, 1965) and relies on approximating the many-body wave function by a product of on-site only contributions. This mean-field method was later extended to the case of bosons Rokhsar and Kotliar (1991); Krauth et al. (1992) and applied to the Bose-Hubbard model of an ultracold bose gas in an optical lattice Jaksch et al. (1998). The Gutzwiller approach relies on the assumption that quantum correlations are small Lewenstein et al. (2012); Pethick and Smith (2008), which is valid in the limit of large particle numbers and/or weak correllations between lattice sites (i.e. weak tunnelling). In the case of infinite dimensions, this treatment becomes exact; however, for lower dimensional cases, and especially for one dimension, there can be significant quantum correlations. While we will consider a one-dimensional model, we will work with a large particle number per site and small tunnelling strengths. In this limit, correlations between sites should be small and the Gutzwiller approach valid. It is also known that the Gutzwiller approach provides qualitatively correct results in one dimension but can not be trusted for the quantitative prediction of critical points of phase transitions, often over-estimating numerical values due to the neglect of correlations Lewenstein et al. (2012); Zwerger (2003).

For bosonic atoms in an optical lattice, this mean-field approach is a suitable approximation to qualitatively capture the allowed phases since the majority of atoms will be in the condensed state. As such, operators may then be expressed as an average around some fluctuating operator,

(17) |

where denotes small deviations. Using Eq. (17), the nearest-neighbour summation terms in Hamiltonian (8) may then be decoupled to a problem which is on-site by linearising the fluctuation field. A natural extension of this mean-field is then to generalise the structure of the many-body wave function under these assumptions Zakrzewski (2005), which can be performed by utilising the Gutzwiller wave function,

(18) |

where is the maximum number of atoms allowed in each lattice site, denotes the size of the lattice, is the state of atoms in site , and are the coefficients which denote the mean-field wave function (commonly referred to as the Gutzwiller coefficients) and they are normalised such that

(19) |

The many-body wave function (18) is a product of on-site states such that we can rewrite the operators in Hamiltonian (8) which are over the nearest-neighbours in terms of on-site only operators. For the single particle tunnelling this is given by

(20) |

the nearest-neighbour interaction by

(21) |

the pair tunnelling by

(22) |

and the density-dependent tunnelling by

(23) | ||||

From taking these nearest-neighbour terms to on-site only terms, we can see that there are four distinct, independent expectation values that are required, , , , and . It comes as no surprise that these four expectations denote the four order parameters of Hamiltonian (8). The first order parameter characterises the single-atom transport properties which we will label as and is given by

(24) |

Next, there is an order parameter which denotes the density behaviour of a phase, which we will label as and is given by

(25) |

There is also an order parameter which defines the pair transport properties of a phase, which we will label as and is denoted by

(26) |

Finally, there is an order-parameter which defines the density-dependent transport properties, which we will label as and is given by

(27) |

Physically, the order parameters represent observables that can be measured in the system. Note, each order parameter is defined for each local lattice site, giving a total of values for each individual order parameter. However, due to the homogeneous nature of the lattice we consider, each of the order parameters will usually be the same, or 2-periodic.

With these relations, Hamiltonian (8) can then be written as a sum of on-site, mean-field Hamiltonians,

(28) |

with each on-site Hamiltonian being given by

(29) | ||||

where we have denoted nearest-neighbour summations by , i.e.

(30) |

In this form, the ground state phase diagrams can be determined using a variety of methods. We will consider the approach of a self-consistent loop. However, it is also possible to use an imaginary time propagation approach, and we have confirmed that our results are consistent with this.

, Phase | Abbreviation | ||||
---|---|---|---|---|---|

Mott Insulator | MI | ||||

Superfluid | SF | ||||

Density Wave | DW | ||||

Supersolid | SS | ||||

One-body Staggered Superfluid | OSSF | ||||

One-body Staggered Supersolid | OSSS | ||||

Pair Superfluid | PSF | ||||

Pair Supersolid | PSS | ||||

Staggered Phases | SP | ||||

Intermediate Staggered Superfluid | ISSF | ||||

Staggered Superfluid | SSF | ||||

Staggered Supersolid | SSS |

The self-consistent loop uses an exact diagonalisation scheme for the on-site problem. We initialise the loop by taking random and uniformly distributed order parameters for each site in the range . Then, the local Hamiltonians are diagonalised such that the order parameters can be redefined from the ground state. Therefore, each local order parameter is updated when each local Hamiltonian is solved. After diagonalising the local Hamiltonians and updating the corresponding order parameters, the loop is then repeated until the energy and order parameters have converged to a given accuracy. We will converge to a accuracy of . From our results, the mean-field approach is stable, with it being rare that the convergence gets stuck in local minima corresponding to excited states. All phases discussed in this work have been checked for a number of iterations of the random initial order parameters to ensure that the phase diagrams are reflecting the ground state properties of the system.

## Iv Characterisation of phases

Before discussing the full phase diagrams, we will first outline the individual phases that will appear and their relation to the four order parameters defined in Sec. III. We summarise all phases that are discussed in this work in Table. 1, where vectors define the corresponding order parameters for a particular lattice site, i.e . For the considered phases, there are at most two unique terms for each set of order parameters. This is represented in compact, periodic form as . In the standard Bose-Hubbard model, phases are defined by homogeneous order parameters. The translational symmetry in Hamiltonian (8) is broken by the two-site two-body terms, allowing for phases with two-site periodic order parameters.

For certain parameter regions, it is expected that we will observe the known Mott-insulator and superfluid phases. The Mott-insulator (MI) is defined by its fixed dynamics (no transport) and an integer valued uniform density across the lattice, i.e. and Fisher et al. (1989). However, the superfluid (SF) phase is given by its uniform non-zero transport property and uniform non-integer density, i.e. and . It is known that the introduction of nearest-neighbour interactions can result in density wave (DW) and supersolid (SS) phases Rossini and Fazio (2012); Kawaki et al. (2017); Ohgoe et al. (2012). The density wave phase is characterised by zero transport properties and a staggered (2-period) density, i.e. and , whereas the supersolid phase is defined by both staggered transport properties and staggered density, i.e. and .

Of course, the presence of the density-dependent and two-body pair tunnelling terms introduce more exotic phases. These phases are a result of the translational symmetry breaking two-site terms, including the density-dependent tunnelling, and the introduction of two-body dynamics by the pair tunnelling. Such transport staggered and pair superfluid phases have been previously observed for Hamiltonian (8) Maik et al. (2013). The one-body staggered superfluid (OSSF) is characterised by a sign staggering in the one-body transport properties, a constant two-body transport, and a constant non-integer density, i.e. , , and . Similar to the ordinary supersolid, the one-body staggered supersolid (OSSS) has staggered one- and two-body transport, and a staggered density, however, the staggering is no longer only characterised by a sign flip. Additionally, it is possible to observe the pair superfluid and supersolid (PSF and PSS), which have no single-particle transport properties and a non-zero two-body transport. The dynamics of these phases are therefore dominated by two-body processes.

In this paper, we observe a new unconventional set of phases, which we will label as staggered phases (SP). These phases are the most general of phase, with all four order parameters staggered. The region of staggered phases is mainly made up of two seperate phases with different symmetries in the density which exhibit an extended second-order (continuous) phase transition between them Landau and Lifshitz (1958). The staggered supersolid (SSS) phase has all four order parameters staggered, whereas the superfluid phase (SSF) has a constant density. Note, that as the transition between these two phases is of second order across a large region in parameter space, the labelling of a staggered superfluid and staggered supersolid phase is not clear in intermediate regions, and we can not rule out the possibility of multiple other phases existing. These regions of neither staggered superfluid or staggered supersolid phase could be a result of the second-order phase transition or the signature of more exotic phases. We expect that these staggered phases exhibited by the Gutzwiller mean-field could be related to the previously observed twisted complex staggered phases in DMRG studies Jürgensen et al. (2015); Lühmann (2016). However, an extension of the mean-field approach used here would be required to compare the staggered and twisted phases, as the twisted phases exhibit complex order parameters, which is outside the scope of this work. That is, due to the nature of the static Gutzwiller mean-field approach, where non-local quantum correlations neglected, care is required when considering initial states with phase differences between local lattice sites.

## V Ground State Phase Diagrams

In this section, mean-field phase diagrams are presented and discussed for various parameter regimes of Hamiltonian (8). To accurately resolve the phase boundaries, we use a grid of at least points for each phase diagram. We will work in units of the two-body on-site interaction strength and consider one-dimensional equally spaced lattices of sites with periodic boundary conditions. The size of the truncated basis, denoted by the maximum number of particles per site , is selected so that the convergence of order parameters is constant with respect to the desired convergence precision when is increased. We have found that in the considered site lattice, a maximum particle number of is required to have machine precision in convergence.

First, we check that for the case of , the well-known Mott-insulator to superfluid transition is observed. We indeed see this behaviour in Fig. 2, with the distinct Mott-insulator lobes which are characterised by different integer uniform fillings of the lattice sites. The critical point of the Mott-insulator to superfluid transition is found to be in agreement with previous Gutzwiller mean-field approaches Zwerger (2003); Lewenstein et al. (2012), with a critical transition point of for the first order Mott-insulator. As expected from mean-field results this is an over-estimation of the true critical point Kashurnikov and Svistunov (1996); Kühner et al. (2000).

We also confirm that given a non-zero , we observe the known supersolid and density wave phases, as shown in Fig. 3. For sufficiently strong , the Mott-insulating phases are completely destroyed and replaced with the density wave phase. This transition to density wave phases makes sense, as the nearest-neighbour interactions can be reduced by generating an offset between the density of nearest-neighbours. Therefore, there is a symmetry breaking of the ground state in order to reduce its energy. It is also observed that the higher order density wave phases exist over a consistent size of parameters, i.e. the area of the density phase is similar for different orders of the density wave. With increasing the superfluid phase is also replaced but with a supersolid. The supersolid phase exists because of the same symmetry arguments already invoked for the density wave, but starting from a superfluid.

A non-zero density-dependent tunnelling causes the destruction of the Mott-insulator phase, as we observe in Fig. 4. This is due to there being an incentive for the state to spread out its density, and favours the superfluid. That is, it is energetically favourable for the density of the ground state to be lower than that required for the Mott-insulator state. With , the magnitude of density-dependent tunnelling enlarges, thus reducing the stability of insulating phases and inflating the overall tunnelling rate.

If there is a non-zero pair tunnelling strength, then the non-trivial pair superfluid and staggered phases are observed, as seen in Fig. 5. The pair superfluid effectively replaces the Mott-insulating lobes at high enough chemical potential. This makes sense, as the system now has a non-zero pair tunnelling, and hence, the pair tunnelling order parameter can not remain zero at large chemical potential. The staggered phases arise mostly in the region of the phase diagram that usually consists of a superfluid. Due to the process of pair tunnelling, it is more likely for large that density and transport properties will clump together and favour some lattice sites over the others. This asymmetry in the density and transport properties results in the staggered phases, where all order parameters are in general staggered.

In Fig. 6a, we study a modification of the density dependent tunnelling process with the opposite sign to the linear tunnelling strength, corresponding to the presence of an attractive long-range interaction. This leads to a skewed structure of the phase diagram, with the Mott-insulating lobes looking more like elliptical structures for large . Furthermore, both a superfluid and one-body staggered superfluid phase are shown to exist. The one-body staggered superfluid phase arises in regions where , i.e. where the attractive non-local interactions are dominating the process, and hence a breaking of the translational symmetry is energetically favourable.

Combinations of several processes are then tuned within Figs. 6b and 6c in order to display the long-range (crystalline order) counterparts of the unconventional one-body staggered superfluid and pair superfluid phases. In Fig. 6b, both neighbour interactions and opposite sign density-dependent tunnelling are considered simultaneously. This leads to the familiar inclusion of supersolid and density wave phases that were already considered. However, we also observe a small parameter region where the one-body staggered supersolid phase exists. This phase is due to the staggering of the transport order parameters being favourable but the tunnelling being strong enough and the chemical potential small enough to not yet favour the one-body staggered superfluid phase. To access this particular phase diagram in a physical set-up, one would require the presence of both repulsive and attractive interactions simultaneously (with suitable tuning). Finally, in Fig. 6c we again demonstrate density wave and supersolid phases for when both neighbour interactions and pair tunnelling are present. We also observe the non-trivial phases of the pair superfluid and supersolid due to the favouring of pair transport for large numbers of atoms (large ), with no single-particle transport present in the ground states.

From the considered phase diagrams, we know that negligible is required in order to stabilise the long-range, symmetry breaking staggered phases. In particular, we must have to observe the staggered phases which, in a purely dipolar setup, requires carefully balanced repulsive and attractive long-range interactions as the relation given in Eq. (16) is not satisfied. However, synthetic many-body processes induced by light-matter interactions Caballero-Benitez and Mekhov (2016); Dogra et al. (2016); Niederle et al. (2016); Landig et al. (2016); Flottat et al. (2017) could alternatively be used to control the allowed phase transitions with greater freedom in an experimental scenario.

It is also worthwhile to consider the staggered phase regions in more detail. We will focus here on the staggered phases appearing in Fig. 5c, which has strong pair tunnelling and no long-range interactions or density-dependent tunnelling. In Fig. 7, we consider the staggered phase and label regions of staggered supersolid, staggered superfluid, and an intermediate staggered superfluid (ISSF). All transitions between each staggered phase is second-order and in Fig. 7a, we label regions where certain phases are dominant. The intermediate staggered superfluid is characterised by a constant density but with a pair transport which is staggered between different values, i.e. . This peculiar property of the intermediate phases are shown in the pair and density modulation order parameter plots of Fig. 7b and Fig. 7c respectively. The staggered superfluid is defined when both the pair and density modulation is equal to zero, whereas the staggered supersolid has finite modulation. In the intermediate case, there is finite pair modulation but zero density modulation.

It would be natural to think that the intermediate staggered superfluid is only a signature of the second-order phase transition. Regardless, it is interesting that in the second-order transition between the staggered superfluid and supersolid the density and pair order parameters appear to change their symmetry properties on different scales, resulting in the dominance of this third intermediate phase. However, it is observed that for moderately large and , the dominant phase is the intermediate staggered superfluid. This phase is of a superfluid nature with a constant non-integer density and a non-zero staggered transport property. The symmetry breaking of the pair tunnelling is a result of the two-body nature of the pair tunnelling process. This phase is similar to the staggered superfluid but it is of a more general type, i.e. all staggered superfluids can also be classified as being in the intermediate phase but not all intermediate staggered superfluids are of the staggered superfluid type.

## Vi Conclusions

In this work, we have considered the ground state phases of long-range interacting bosonic ultracold gases in optical lattices. This model is the most general of Bose-Hubbard models for two-body interactions. In addition to the standard local terms, this model includes pair tunnelling, long-range two-body interactions, and a density-dependent (or induced) tunnelling. We consider the ground state phases of this Bose-Hubbard model via a Gutzwiller local mean-field approach, which is valid in the limit of small quantum correlations. Therefore, it is necessary to consider large Hilbert spaces and small tunnelling strengths to be in the region where the mean-field is valid.

We confirm the presence of density wave and supersolid phases with intermediate and strong nearest-neighbour two-body interactions. In addition, we observe the known destruction of the Mott-insulating phase by a density-dependent tunnelling process. By considering the behaviour of the system with non-zero pair tunnelling, we observe an interesting mixed-phase region where all order parameters, or all but the density, are staggered. This region consists of a mixture between superfluid and supersolid phases, with second-order phase transitions between them over a large region of parameter space.

In summary, we have provided a detailed study of the ground state phases of the general Bose-Hubbard model for long-range interacting atoms. We observed the existence of new, 2-site periodic phases in the 1D limit. From the corresponding order parameters, we have observed not only unconventional sign staggering for both one-body and pair superfluids, but a rich, non-trivial structure in the symmetry transitions of the staggered phases. While these results are general, and not specific to an atomic species or set-up, we expect that the non-trivial ground states observed here could play an important role for dipolar atomic gases in optical lattices. In particular, the new staggered phases could be observed when dipolar gases are combined with synthetic techniques, e.g. light-matter interactions Caballero-Benitez and Mekhov (2016); Dogra et al. (2016); Niederle et al. (2016); Landig et al. (2016); Flottat et al. (2017), to induce competing long-range many-body interactions.

###### Acknowledgements.

D.J., N.W. and C.W.D. acknowledge support from EPSRC CM-CDT Grant No. EP/L015110/1. P.Ö. acknowledges support from EPSRC Grant No. EP/M024636/1.## References

- Gutzwiller (1963) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
- Kanamori (1963) J. Kanamori, Prog. Theor. Phys. 30, 275 (1963).
- Hubbard (1963) J. Hubbard, P. Roy. Soc. Lond. A Mat. 276, 238 (1963).
- Hubbard (1964) J. Hubbard, P. Roy. Soc. Lond. A Mat. 281, 401 (1964).
- Jaksch et al. (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
- Bloch (2005) I. Bloch, Nat. phys. 1, 23 (2005).
- Lewenstein et al. (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, and U. Sen, Adv. Phys. 56, 243 (2007).
- Dutta et al. (2015) O. Dutta, M. Gajda, P. Hauke, M. Lewenstein, D.-S. LÃ¼hmann, B. A. Malomed, T. Sowiński, and J. Zakrzewski, Rep. Prog. Phys. 78, 066001 (2015).
- Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
- Freericks and Monien (1994) J. K. Freericks and H. Monien, Europhys. Lett. (EPL) 26, 545 (1994).
- Greiner et al. (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, nature 415, 39 (2002).
- Zwerger (2003) W. Zwerger, J. Opt. B Quantum S. O. 5, S9 (2003).
- Bakr et al. (2010) W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Foelling, L. Pollet, and M. Greiner, Science 329, 547 (2010).
- Fedichev et al. (1996) P. O. Fedichev, Y. Kagan, G. V. Shlyapnikov, and J. T. M. Walraven, Phys. Rev. Lett. 77, 2913 (1996).
- Inouye et al. (1998) S. Inouye, M. Andrews, J. Stenger, H.-J. Miesner, D. Stamper-Kurn, and W. Ketterle, Nature 392, 151 (1998).
- Theis et al. (2004) M. Theis, G. Thalhammer, K. Winkler, M. Hellwig, G. Ruff, R. Grimm, and J. H. Denschlag, Phys. Rev. Lett. 93, 123001 (2004).
- Dickerscheid et al. (2005) D. B. M. Dickerscheid, U. Al Khawaja, D. van Oosten, and H. T. C. Stoof, Phys. Rev. A 71, 043604 (2005).
- Miyake et al. (2013) H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Phys. Rev. Lett. 111, 185302 (2013).
- Wozniak et al. (2018) D. Wozniak, F. M. Dobler, and A. Posazhennikova, Phys. Rev. A 98, 053628 (2018).
- Menotti et al. (2008) C. Menotti, M. Lewenstein, T. Lahaye, and T. Pfau, AIP Conf. Proc. 970, 332 (2008).
- Lucioni et al. (2018) E. Lucioni, L. Tanzi, A. Fregosi, J. Catani, S. Gozzini, M. Inguscio, A. Fioretti, C. Gabbanini, and G. Modugno, Phys. Rev. A 97, 060701 (2018).
- Griesmaier (2007) A. Griesmaier, J. Phys. B At. Mol. Opt. 40, R91 (2007).
- Kremer et al. (2017) M. Kremer, R. Sachdeva, A. Benseny, and T. Busch, Phys. Rev. A 96, 063611 (2017).
- Aikawa et al. (2012) K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 (2012).
- Lu et al. (2012) M. Lu, N. Q. Burdick, and B. L. Lev, Phys. Rev. Lett. 108, 215301 (2012).
- Trefzger et al. (2011) C. Trefzger, C. Menotti, B. Capogrosso-Sansone, and M. Lewenstein, J. Phys. B At. Mol. Opt. 44, 193001 (2011).
- Maik et al. (2013) M. Maik, P. Hauke, O. Dutta, M. Lewenstein, and J. Zakrzewski, New J. Phys. 15, 113041 (2013).
- Łącki and Zakrzewski (2013) M. Łącki and J. Zakrzewski, Phys. Rev. Lett. 110, 065301 (2013).
- Sowiński et al. (2012a) T. Sowiński, O. Dutta, P. Hauke, L. Tagliacozzo, and M. Lewenstein, Phys. Rev. Lett. 108, 115301 (2012a).
- Sowiński et al. (2012b) T. Sowiński, O. Dutta, P. Hauke, L. Tagliacozzo, and M. Lewenstein, Phys. Rev. Lett. 108, 115301 (2012b).
- Mosadeq and Asgari (2015) H. Mosadeq and R. Asgari, Phys. Rev. B 91, 085126 (2015).
- Rossini and Fazio (2012) D. Rossini and R. Fazio, New J. Phys. 14, 065012 (2012).
- Kawaki et al. (2017) K. Kawaki, Y. Kuno, and I. Ichinose, Phys. Rev. B 95, 195101 (2017).
- Ohgoe et al. (2012) T. Ohgoe, T. Suzuki, and N. Kawashima, Phys. Rev. B 86, 054520 (2012).
- Biedroń et al. (2018) K. Biedroń, M. Łącki, and J. Zakrzewski, Phys. Rev. B 97, 245102 (2018).
- Jürgensen et al. (2015) O. Jürgensen, K. Sengstock, and D. S. Lühmann, Sci. Rep. 5, 12912 (2015).
- Rosenberg (1988) H. M. Rosenberg, The Solid State, 3rd ed. (Oxford University Press, 1988) pp. 140–141.
- Gutzwiller (1964) M. C. Gutzwiller, Phys. Rev. 134, A923 (1964).
- Gutzwiller (1965) M. C. Gutzwiller, Phys. Rev. 137, A1726 (1965).
- Rokhsar and Kotliar (1991) D. S. Rokhsar and B. G. Kotliar, Phys. Rev. B 44, 10328 (1991).
- Krauth et al. (1992) W. Krauth, M. Caffarel, and J.-P. Bouchaud, Phys. Rev. B 45, 3137 (1992).
- Lewenstein et al. (2012) M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold Atoms in Optical Lattices: Simulating quantum many-body systems (Oxford University Press, 2012) pp. 69–72.
- Pethick and Smith (2008) C. J. Pethick and H. Smith, Bose–Einstein condensation in dilute gases (Cambridge university press, 2008) pp. 433–439.
- Zakrzewski (2005) J. Zakrzewski, Phys. Rev. A 71, 043601 (2005).
- Landau and Lifshitz (1958) L. Landau and E. Lifshitz, Statistical Physics (Course of Theoretical Physics vol 5), 1st ed. (Addison-Wesley publishing company, 1958) pp. 430–456.
- Lühmann (2016) D.-S. Lühmann, Phys. Rev. A 94, 011603 (2016).
- Kashurnikov and Svistunov (1996) V. A. Kashurnikov and B. V. Svistunov, Phys. Rev. B 53, 11776 (1996).
- Kühner et al. (2000) T. D. Kühner, S. R. White, and H. Monien, Phys. Rev. B 61, 12474 (2000).
- Caballero-Benitez and Mekhov (2016) S. F. Caballero-Benitez and I. B. Mekhov, New J. Phys. 18, 113010 (2016).
- Dogra et al. (2016) N. Dogra, F. Brennecke, S. D. Huber, and T. Donner, Phys. Rev. A 94, 023632 (2016).
- Niederle et al. (2016) A. E. Niederle, G. Morigi, and H. Rieger, Phys. Rev. A 94, 033607 (2016).
- Landig et al. (2016) R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner, and T. Esslinger, Nature 532, 476 (2016).
- Flottat et al. (2017) T. Flottat, L. d. F. de Parny, F. Hébert, V. G. Rousseau, and G. G. Batrouni, Phys. Rev. B 95, 144501 (2017).