Bose-Hubbard model and beyond with MCTDH-X

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 MCTDH-X enables a modelling of the many-body state with an improved accuracy as compared to the Hubbard-lattice-models, because it uses a general, variationally optimized and not necessarily site-local basis set. This improved accuracy is of particular importance in cases where the single-band Hubbard description fails, for instance, due to the action of a blue-detuned cavity [3] or long-range dipolar interactions [5,6,4]. Direct comparisons between the MCTDH-X and Hubbard models have been performed in Refs. [12,9,15,7,8,13,14,16,11,10].

The Bose-Hubbard Hamiltonian takes the following form in one dimension,

$$\hat{H}_{\text{BH}} = -t\sum_{\langle i,j\rangle} (\hat{c}^\dagger_i \hat{c}_j+\text{H.c.})+ U\sum_i \hat{n}_i(\hat{n}_i-1),$$ (4)

where, $\hat{c}_i$ is the annihilation operator for a particle in the Wannier state $w_i(x)$ that is localized at site $i$ of the lattice and $\hat{n}_i=\hat{c}^\dagger_i\hat{c}_i$ is the corresponding number operator. In the Hubbard model, the superfluid to Mott-insulator transition is a result of the competition between the hopping strength $t\propto \int dx w^*_i(x) \partial_x^2w_{i+1}(x+a)$ between lattice sites with lattice constant $a$, and the on-site interaction $U\propto\int dx \vert w_i(x)\vert^4$ between particles. We consider a system with a fixed total number of particles $N$ and therefore neglect the chemical potential.

In the limit $t/U\rightarrow\infty$, the system is in the superfluid phase, and the ground state is given by $\vert\Psi\rangle \propto \left(\sum_i a^\dagger_i\right)^N\vert\rangle$. In this state, the $U(1)$ symmetry $\hat{c}_i\rightarrow e^{i\theta}\hat{c}_i$ 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 single-particle state. This condensation entails a Glauber one-body correlation between sites which is unity. In the limit $U/t\rightarrow\infty$, the system is in a Mott insulator phase where the ground state of the system is given by $\vert\Psi\rangle \propto \prod_i \left(a^\dagger_i\right)^{n_i}\vert\rangle$, where $n_i$ is the number of particles in the $i$-th site and satisfies the constraint $\sum_in_i=N$. 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 one-body correlation between sites vanishes. The transition point given in units of $t/U$ depends on the structure of the lattice. In particular, it is roughly proportional to the filling factor and the inverse of the dimensionality [17].

Figure S1: (a) The real-space density distribution for $N=2$ bosons with extremely strong interactions $g=500$ performed with different orbital numbers $M=$1, 2, 15. As the number of orbitals increases, the simulated density profile becomes closer and closer to the density profile of fermions. (b,c) The two-body correlation function inside one well $g^{(2)}(x>0,x'>0)$ of a double-well potential simulated with $N=6$ bosons and (b) $M=6$ orbitals (c) $M=2$ orbitals. More details are unveiled by the extra orbitals in panel (a), indicating physics beyond the single-band Hubbard model.

For illustration, we consider $N=6$ weakly-interacting ($g=1$) bosons in a double-well potential with high barrier $E_$dw$=20$, where the two wells are Mott insulating from each other. According to the single-band Bose-Hubbard 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 many-body state can thus be written as $\vert\Psi\rangle \propto \left(\hat{c}^\dagger_L\right)^3\left(\hat{c}_R^\dagger\right)^3\vert\rangle$, where $\hat{c}^\dagger_L$ and $\hat{c}^\dagger_R$ 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 one-body correlation function [$g^{(1)}$][cf. Fig. 2(e,f) in the main text], it contradicts with the two-body correlation function [$g^{(2)}$] obtained from the simulations. In a Mott insulator state where $g^{(1)}$ between lattice sites vanishes, the diagonal term of $g^{(2)}$ is directly related to the number of particles $n_i$ in the corresponding site through $g^{(2)}(x,x) = 1-1/n_i$ [16]. If the single-band Bose-Hubbard model holds true, the diagonal term of $g^{(2)}$ would be uniformly $g^{(2)}(x,x) = 2/3$ since $n_i=3$. 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 $M=2$ orbitals are used in the simulations, indeed $g^{(2)}(x,x) = 2/3$ holds uniformly [Fig. S1(c)]. These results show that the higher-order ($M>3$) orbitals are indispensable for correctly capturing the particle correlations, and indicate that more than one mode has contributions in each well. The single-band Bose-Hubbard model is an over-simplification for the system. We have seen that although the third orbital $\rho_3=4.432\times10^{-3}$ is two order of magnitude lower than the second orbital $\rho_2=0.493$, 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.