The Hubbard model usually serves as good approximation for a lattice system. It assumes only one Wannier function contributes in each lattice site. We note that the MCTDHX enables a modelling of the manybody state with an improved accuracy as compared to the Hubbardlatticemodels, because it uses a general, variationally optimized and not necessarily sitelocal basis set. This improved accuracy is of particular importance in cases where the singleband Hubbard description fails, for instance, due to the action of a bluedetuned cavity [3] or longrange dipolar interactions [5,6,4]. Direct comparisons between the MCTDHX and Hubbard models have been performed in Refs. [12,9,15,7,8,13,14,16,11,10].
The BoseHubbard Hamiltonian takes the following form in one dimension,
(4) 
In the limit , the system is in the superfluid phase, and the ground state is given by . In this state, the symmetry of the system is broken spontaneously, such that all the particles have the same phase. Particles in different sites are coordinated; they are condensed in the same singleparticle state. This condensation entails a Glauber onebody correlation between sites which is unity. In the limit , the system is in a Mott insulator phase where the ground state of the system is given by , where is the number of particles in the th site and satisfies the constraint . In this state, particles in different sites are disconnected and have no information on each other's behaviors. The phases of the particles in different sites lose coherence. Strong correlation effects build up and the onebody correlation between sites vanishes. The transition point given in units of depends on the structure of the lattice. In particular, it is roughly proportional to the filling factor and the inverse of the dimensionality [17].

For illustration, we consider weaklyinteracting () bosons in a doublewell potential with high barrier dw, where the two wells are Mott insulating from each other. According to the singleband BoseHubbard model, the energy of this system is minimzed when the bosons are evenly distributed into the two wells, with three bosons in each well. The manybody state can thus be written as , where and are the creation operators for the band in the left and right well, respectively. Although such a picture is well supported by the momentum space density distribution and the onebody correlation function [][cf. Fig. 2(e,f) in the main text], it contradicts with the twobody correlation function [] obtained from the simulations. In a Mott insulator state where between lattice sites vanishes, the diagonal term of is directly related to the number of particles in the corresponding site through [16]. If the singleband BoseHubbard model holds true, the diagonal term of would be uniformly since . This is not the case in the simulations when sufficiently many orbitals are used. The landscape of the correlation function is uneven, with values ranging from 0.4 to 0.8 [Fig. S1(b)]. On the contrary, if only orbitals are used in the simulations, indeed holds uniformly [Fig. S1(c)]. These results show that the higherorder () orbitals are indispensable for correctly capturing the particle correlations, and indicate that more than one mode has contributions in each well. The singleband BoseHubbard model is an oversimplification for the system. We have seen that although the third orbital is two order of magnitude lower than the second orbital , it has significant impact on the simulated particle correlations.
In summary, to decide whether the convergence in orbital number has been achieved, we should perform simulations with more orbitals and compare the convergence in different quantities of interest.