Quantum phase diagram and time-of-flight absorption pictures of an ultracold Bose system in a square optical superlattice

The physics of ultracold atomic systems in optical lattices has been boosted as one of the major and fascinating fields during the past decade [1,2], and has received tremendous attention on account of the novelty and various potential applications [3,4], ranging from condensed-matter physics [5,6], quantum information processing [7–9], even to cosmological problems [10].

While many experiments with ultracold atoms have been performed in homogeneous simple cubic lattices, more and more experimental and theoretical efforts have been devoted to complex systems, including ultracold bosons in non-rectangular optical lattices [11–13], systems with long-range interactions [14,15], or with multi components [16]. Among them, an interesting topic is ultracold Bose systems in optical superlattices [17,18]. Due to the complexity of the lattice structures, the quantum phase transitions and corresponding phase diagrams of these systems are expected to be more complex, and hence need to be investigated analytically.

Actually, beside the mean-field theory [19], which underestimates the phase boundaries when compared with Monte Carlo results [20], and the strong coupling expansion method [21], which overestimates the phase boundaries, an alternative field-theoretical method, the so-called effective potential theory [22], has been developed in recent years. By utilizing it, the quantum phase diagrams of ultracold Bose systems in triangular, hexagonal, and Kagomé optical lattices have been investigated [23], the relative difference between our analytical results and numerical results [24,25] is less than $10\%$.

Meanwhile, momentum distribution functions, which can be revealed via time-of-flight (TOF) technique [1,3,11,26], should be calculated as well in order to compare the theoretical results with experimental data. Thanks to the re-summed Green's function method [27] and cumulant expansion [28], we have calculated out analytically the TOF absorption pictures for ultracold Bose systems in triangular optical lattices for various lattice depth [29], exhibiting a good agreement with experimental data [11].

In this letter, we will investigate analytically the quantum phase transitions of an ultracold Bose system in a square superlattice. In fact, Wang et al. [30] have developed a generalized effective potential theory and calculated out the phase diagram of an ultracold superlattice Bose system up to second order. With the help of the process chain technique [31], by the use of the Kato formulation [32] of perturbation, we will push the limit of the perturbative calculation of quantum phase boundaries of such superlattice ultracold Bose system up to the 8th order, and further goes to the order of infinity via the technique of linear fit extrapolation. Moreover, by making use of the Green's function method, we will calculate the TOF absorption pictures of the superlattice ultracold Bose system.

A system of spinless bosons trapped in a square superlattice is described by [30,33]

$$\begin{eqnarray} \hat{H}_{SL}&= -t\sum_{\left\langle i\in A,j\in B\right\rangle }\left(\hat{a}_{i}^{\dagger}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\hat{a}_{i}\right) +\frac{U}{2}\sum_{i\in A,B}\hat{n}_{i}\left(\hat{n}_{i}-1\right) \nonumber\\ &\quad -\,\left(\mu+\Delta\mu\right)\sum_{i\in A}\hat{n}_{i}-\mu\sum_{i\in B}\hat{n}_{i}, \end{eqnarray} \tag{ 1 }$$

where t is the nearest-neighbor hopping parameter, U denotes the on-site repulsion between two atoms. $\hat{a}_{i}^{\dagger}\ (\hat{a}_{i})$ represents the boson creation (annihilation) operator at site i, $\Delta \mu$ stands for the difference of the chemical potentials on sublattices A and B and without loss of generality, we set $\Delta \mu>0$.

The subtle balance among the hopping parameter t, on-site repulsion U, and $\Delta \mu$ introduces rich phase structures to this system. Even when the hopping parameter t = 0, the competition between U and $\Delta \mu$ will lead to different quantum phases at zero temperature, these phases correspond to different regions of the chemical potential.

Before we discuss the quantum phase boundaries of the system, the ground state of the site-diagonal part of the above Hamiltonian

$$\begin{equation} \hat{H}_{0}=\frac{U}{2}\sum_{i\in A,B}\hat{n}_{i}\left(\hat{n}_{i}-1\right)-(\mu+\Delta \mu)\sum_{i\in A}\hat{n}_{i}-\mu\sum_{i\in B}\hat{n}_{i} \end{equation} \tag{ 2 }$$

needs to be determined. Since all terms in this part of Hamiltonian are site-diagonal, its eigenstates can be denoted in terms of the occupation number of each sublattice, i.e. $|n_A,n_B\rangle$, the corresponding eigenenergies read

$$\begin{eqnarray} E_{n_{A},n_{B}}&= \frac{U}{2}n_{A}(n_{A}-1)+\frac{U}{2}n_{B}(n_{B}-1)\nonumber\\ &\quad -\mu(n_{A}+n_{B})-\Delta \mu n_A. \end{eqnarray} \tag{ 3 }$$

Mott insulator phases correspond to the cases when $n_A=n_B$, while charge-density-wave (CDW) phases when $n_A\neq n_B$. As is known, the condition for the existence of a Mott phase $|n,n\rangle$ is $E_{n,n}\leq E_{n+1,n}$ and at the same time $E_{n,n}\leq E_{n,n-1}$, from eq. (3), we see that this leads to

$$\begin{equation} U(n-1)\leq \mu \leq Un -\Delta \mu. \end{equation} \tag{ 4 }$$

It is easy to recognize that there is no Mott phase when $\Delta \mu>U$, and the ground state for t = 0 at any case would be CDW phases. Hence, without loss of generality, we would focus on the case of $0<\Delta \mu <U$ in the rest of the paper. In this case, the occupation number of these two sublattices in possible CDW phases may only be $n_B=n_A-1$, and the corresponding regions of μ read

$$\begin{equation} U(n-1)-\Delta \mu\leq \mu \leq U(n-1). \end{equation} \tag{ 5 }$$

With this information in hand, we can now go further to tune on the hopping parameter to see the quantum phase transitions from the incompressible Mott states or CDW states to superfluid phase. In order to investigate the quantum phase transitions of the superlattice Bose system and determine the corresponding phase boundaries analytically, similar to previous work [22,23], we would like to use the generalized effective potential method [30] under the framework of the Ginzburg-Landau field theory, i.e. we add, for the moment, additional source terms with strength $\textbf{J}=(\textbf{J}_{A}, \textbf{J}_{B})^{T}$ into the above superlattice Bose-Hubbard Hamiltonian as follows:

$$\begin{eqnarray} \hat{H}_{SL}\left(\textbf{J},\textbf{J}^{\dagger}\right)&= \hat{H}_{SL}+\sum_{j\epsilon A}\left(J_{A}\hat{a}_{j}^{\dagger}+J_{A}^{*}\hat{a}_{j}\right)\nonumber\\ &\quad +\,\sum_{j\epsilon B}\left(J_{B}\hat{a}_{j}^{\dagger}+J_{B}^{*}\hat{a}_{j}\right) \end{eqnarray} \tag{ 6 }$$

and treat the hopping term and the additional external source terms as perturbations.

It is quite straightforward to get the grand-canonical free energy of the system as power series of both the hopping parameter t and the external source $\textbf{J}\ (\textbf{J}^{\dagger})$ via Taylor expansion. Since the unperturbed ground states are either Mott phases or CDW states, J and $\textbf{J}^{\dagger}$ can only appear in pair, after re-grouping the terms in the free energy with respect to J and $\textbf{J}^{\dagger}$, the free energy reads

$$\begin{eqnarray} F\left(\textbf{J},\textbf{J}^{\dagger},t\right)&= N_{s}\left[F_{0}\left(t\right)+\textbf{J}^{\dagger}{\mathbb C}_{2}\left(t\right)\textbf{J}\right.\nonumber\\ &\quad +\,\left.\textbf{J}^{\dagger}\textbf{J}^{\dagger}{\mathbb C}_{4}\left(t\right)\textbf{J}\textbf{J}+\cdots\right] \end{eqnarray} \tag{ 7 }$$

with the expansion coefficient tensor ${\mathbb C}_2(t)=\Bigl(\begin{array}{cc} c_{2AA} & c_{2AB}\\ c_{2BA} & c_{2BB} \end{array}\Bigr)$ being $c_{2AA}=\sum_{n=0}^{\infty}(-t)^{n}\alpha_{2AA}^{(n)}$, $c_{2BB}=\sum_{n=0}^{\infty}(-t)^{n}\alpha_{2BB}^{(n)}$, $c_{2AB}=\sum_{n=0}^{\infty}(-t)^{n}\alpha_{2AB}^{(n)}$, N_s is the total number of the lattice sites.

Due to the nature of the superlattice structure of the system, the superfluid order parameter takes the form of the vector $\mathbf{\Psi}=(\psi_A, \psi_B)^T\ (\mathbf{\Psi}^{\dagger}=(\psi_A^*, \psi_B^*))$ whose components are defined as $\psi_m=\langle a_m\rangle\ (\psi_m^*=\langle a_m^{\dagger}\rangle)\ (m=A, B)$ which can be calculated from the free energy by [34] $\psi_{m}=\frac{1}{N_{S}}\frac{\partial F}{\partial J_{m}^{*}}\ (\psi_{m}^{*}=\frac{1}{N_{S}}\frac{\partial F}{\partial J_{m}})\ (m=A,B)$.

A Legendre transformation of the free energy $F(\textbf{J},\textbf{J}^{\dagger},t)$ leads to an effective potential $\Gamma(\mathbf{\Psi},\mathbf{\Psi} ^{\dagger},t)=F/N_{s}-\mathbf{\Psi}^{\dagger}\textbf{J}-\textbf{J}^{\dagger} \mathbf{\Psi}$ which is expressed in terms of $\mathbf{\Psi}\ (\mathbf{\Psi}^{\dagger})$ as

$$\begin{equation} \Gamma\left(\!\mathbf{\Psi},\mathbf{\Psi}^{\dagger},t\right)=F_{0}(t)+\mathbf{\Psi}^{\dagger} {\mathbb A}_{2}(t)\mathbf{\Psi}+\mathbf{\Psi}^{\dagger}\mathbf{\Psi}^{\dagger}{\mathbb A}_{4}\left(t\right)\mathbf{\Psi\Psi}+\cdots, \end{equation} \tag{ 8 }$$

this takes exactly the form of the Landau $\phi^4$ theory. According to the Landau theory, $\det {\mathbb A}_{2}(t)=0$ determines the phase boundaries between the uncompressed Mott or CDW states and superfluid phases. It is not difficult to find out that ${\mathbb A}_{2}(t)=-{\mathbb C}_{2}^{-1}(t)$, hence, the critical value of t can be found by looking for the radius of convergence of the determinant of ${\mathbb C}_{2}(t)$:

$$\begin{equation} \det {\mathbb C}_{2}(t)=c_{2AA}c_{2BB}-c_{2AB}c_{2BA}\equiv\sum_{k=0}^{\infty}D_2^{(2k)}t^{2k}. \end{equation} \tag{ 9 }$$

There is no odd order term of t because $\alpha_{2AA(BB)}$ has only even terms, while $\alpha_{2AB(BA)}$ has only odd terms. The convergence radius may be found via the so-called d'Alembert's ratio test, and the n-th–order approximation of the phase boundary reads

$$\begin{equation} t^{2}_c=\frac{D_2^{(n-2)}}{D_2^{(n)}}=\frac{\alpha_{2AA}^{\left(n-2\right)}\alpha_{2BB}^{\left(0\right)}-\alpha_{2AB}^{\left(n-3\right)}\alpha_{2BA}^{\left(1\right)}+ \cdots-\alpha_{2AB}^{\left(1\right)}\alpha_{2BA}^{\left(n-3\right)}+\alpha_{2AA}^{\left(0\right)}\alpha_{2BB}^{\left(n-2\right)}}{\alpha_{2AA}^{\left(n\right)} \alpha_{2BB}^{\left(0\right)}-\alpha_{2AB}^{\left(n-1\right)}\alpha_{2BA}^{\left(1\right)}+\cdots-\alpha_{2AB}^{\left(1\right)}\alpha_{2BA}^{\left(n-1\right)} +\alpha_{2AA}^{\left(0\right)}\alpha_{2BB}^{\left(n\right)}}. \end{equation} \tag{ 10 }$$

From the above discussion, we see that in order to find the critical value of t to determine the phase boundaries, the coefficients of $\alpha^{(n)}_{2}$ have to be calculated out, this can only be achieved perturbatively. Since the unperturbed state is a Mott state or a CDW state, every nonzero contribution to $\alpha^{(n)}_{2}$ includes exactly one creation operator at site i (associated with J_i) and one annihilation operator at site j (associated with $J^*_j$) as well as n nearest-neighbor hopping processes (associated with t) connecting site i and j.

Actually, when calculating many-body problems, a perturbation theory as the Rayleigh-Schrödinger perturbation expansion [35] may be applied, with the help of Kato's formulation of perturbation series [32] and process chain technique [31], results with very high-order correction may in principle be reached. Let us denote the ground state of $\hat{H}_0$ by $|m\rangle$. In general, when $|m\rangle$ is subject to some perturbation V, the eigenenergy of the total Hamiltonian $\hat{H}=\hat{H}_0+V$ can be perturbatively expressed as $E_m=E_m^{(0)}+\sum_nE_m^{(n)}$ with the n-th–order correction given by the trace [32]

$$\begin{equation} E_{m}^{(n)}=\text{Tr}\left[{\displaystyle\sum_{\left\{ \alpha_{\ell}\right\} }S^{\alpha_{1}}V}S^{\alpha_{2}}VS^{\alpha_{3}}\ldots S^{\alpha_{n}}VS^{\alpha_{n+1}}\right], \end{equation} \tag{ 11 }$$

the operators $S^{\alpha}$ connecting each perturbation terms V read

$$\begin{equation} S^{\alpha}=\left\{ \begin{array}{c@{\quad}c} -\left|m\right\rangle \left\langle m\right|, & \text{for }\alpha=0,\\ {\displaystyle \sum_{i\neq m}\frac{\left|i\right\rangle \left\langle i\right|}{\left(E_{m}^{\left(0\right)}-E_{i}^{(0)}\right)^{\alpha}}}, & \text{for }\alpha>0 \end{array}\right. \end{equation} \tag{ 12 }$$

with $E^{(0)}_{m}$ and $E^{(0)}_{i}$ the unperturbed energies of the $\hat{H}_{0}$'s eigenstates $|m\rangle$ and $|i\rangle$. Since the eigenstates of $\hat{H}_0$ form an orthonormal basis, it is easy to prove that ${S^{0}S^{\alpha}=0}$, for, ${\alpha > 0}$, and $S^{\alpha}S^{\beta}=S^{\alpha+\beta}$ for $\alpha,\beta > 0$. It should be emphasized that all possible sequences of $\{\alpha_{\ell}\}$ with constraint $\sum_{\ell=1}^{n+1}\alpha_{\ell}=n-1$ must be taken into account in the calculation [31,32].

Each Kato term $\langle m|S^{\alpha_{1}}V S^{\alpha_{2}}VS^{\alpha_{3}}\ldots S^{\alpha_{n}}VS^{\alpha_{n+1}}|m\rangle$ in eq. (11) represents a sum of processes going from the ground state $|m\rangle$ over series of different intermediate states $|i\rangle$, caused by the perturbation V, then back to $|m\rangle$, i.e. the Kato terms are sums of the so-called process chains [25,31].

In our case, the perturbation terms are hopping terms and external source terms, and the unperturbed ground state $|m\rangle$ stands for the Mott or CDW state. All these process chains can be abstractly represented by diagrams [31], due to the linked-cluster theorem [35], only connected diagrams contribute. As an example, we show in fig. 1 topologically different diagrams for sixth-order perturbations consisting of exactly one source in (denoted by a dot $(\bullet)$ in the diagram), one source out (denoted by $(\times)$) and four hoppings between nearest-neighbor sites (denoted by arrows).

Zoom In Zoom Out Reset image size

By applying the above sketch process chain method to the calculation of all $\alpha^{(n)}_2$ needed, together with eq. (10), the phase boundaries between the Mott (or CDW) state and the superfluid phase can then be calculated. The quantum phase boundaries of the ultra-cold Bose system in a square superlattice with $\Delta \mu=0.5U$ are shown in fig. 2 up to the 8th order. The comparison to Monte Carlo numerical results [30] shows that the relative deviation of our 8th-order process chain results from the Monte Carlo results is less than $4\%$.

Zoom In Zoom Out Reset image size

A more accurate result may further be obtained via the linear fit extrapolation method [25,31]. When plotting the logarithm of the coefficients $D_2^{(n)}$ against the order of n, i.e. the number of tunneling processes in the perturbative calculation, for different ground states (CDW or MI), we find a good linear behavior, as shown in fig. 3. This indicates that, in a good approximation, the ratio between the coefficients of adjacent orders is constant, and would be exact in the case of infinite dimensionality [31]. However, due to the tunneling perturbation, the ratio would change slightly when increasing the order n. As is known, in order to get the exact convergence radius of a power series via d'Alembert's ratio test, the order n should be sent to infinite. Based on the data of order 2, 4, 6, and 8 obtained above, the extrapolation over $1/n$ can be carried out, and we can then obtain the critical value of t_c for different μ when $n\rightarrow \infty$. As shown in fig. 4, our result is nearly identical with the numerical result [30]. By selecting different orders when doing the extrapolation, we find that the relative error of our result is less than 1%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our previous work [29], with the help of re-summed Green's function method [27,29], we have calculated analytically the TOF pictures of a Bose-Hubbard system in a triangular optical lattice. A generalized Green's function method [36] has also been developed to tackle problems of a Bose-Hubbard system in a bipartite optical lattice. In this section, by utilizing the generalized Green's function method, we will investigate the TOF of the Bose-Hubbard model on a square superlattice.

As is known, the superlattice Bose-Hubbard Hamiltonian in eq. (1) is the single-band approximation of the many-body Hamiltonian $H=\int \text{d}\textbf{r}\psi^{\dagger}(\textbf{r})[-\frac{\hbar^2 \nabla^2}{2m}+V(\textbf{r})]\psi(\textbf{r})+\frac g2\int \text{d}\textbf{r}\psi^{\dagger}(\textbf{r})\psi^{\dagger}(\textbf{r})\psi(\textbf{r})\psi(\textbf{r})$ in the tight-binding limit at extremely low temperature [37], the optical lattice potential $V(\textbf{r})$ is [38]

$$\begin{eqnarray} V_{op}(\textbf{r})&= -V_{0}\left[\cos^{2}\left(\frac{2\pi x}{\lambda}\right)+\cos^{2}\left(\frac{2\pi y}{\lambda}\right)\right. \nonumber\\ &\quad \left.+\,2\cos\theta \cos\left(\frac{2\pi x}{\lambda}\right)\cos\left(\frac{2\pi y}{\lambda}\right)\right] \end{eqnarray} \tag{ 13 }$$

which creates the checkerboard pattern of the square optical superlattice with alternate deep and shallow wells and the relative well depth can be tuned in real time by changing the phase difference $\theta\ (0<\cos\theta <1)$ between the counterpropagating laser beams. The lattice constant $a=\frac{\lambda}{2}$ with λ being the wavelength of the laser.

In the harmonic approximation, we expand the optical lattice potential in the vicinity of sublattice A (the locations of deep wells, the sublattice which the site $(0,0)$ belongs to) and sublattice B (for instance the site $(a,0)$), respectively, and easily get the harmonic-oscillator frequencies $\omega_{A,B}=\sqrt{\frac{V_{0}2\pi^{2}(1\pm\cos\theta)}{ma^{2}}}$ for deep wells (+ sign) and shallow wells (− sign). Correspondingly, in unit of recoil energy $E_R=\pi^2\hbar^2/(2ma^2)$, the lowest-band Wannier functions for sublattices A and B read $w(\textbf{r}-\textbf{r}_{i})|_{i\in A,B}=[\tilde{V}_{0}(1\pm\cos\theta)]^{\frac{1}{4}} (\frac{\pi}{a^{2}})^{\frac{1}{2}} \exp[-\frac{\pi^{2}}{2}\sqrt{\tilde{V}_{0}(1\pm\cos\theta)}\frac{(x-x_{i})^{2}+(y-y_{i})^{2}}{a^{2}}]$ with ${\tilde{V}_{0}=\frac{V_{0}}{E_{R}}}$.

When expanding the field operator $\psi(\textbf{r})$ in the basis of the orthonormal Wannier function $w(\textbf{r}-\textbf{r}_{i})$, $\psi(\textbf{r})=\sum_{i}w(\textbf{r}-\textbf{r}_{i})\hat{a}_i$, we find that the hopping parameter t and the on-site interaction U are determined by [37] $t_{ij}=-\int \text{d}\textbf{r}w^{*}(\textbf{r}-\textbf{r}_{i})(-\frac{\hbar^{2}}{2\pi}\nabla^{2}-V_{op}(\textbf{r}))w(\textbf{r}-\textbf{r}_{j})$ and $U=g\int \text{d}\textbf{r}|w(\textbf{r}-\textbf{r}_{i})|^{4}$, respectively. Meanwhile, from the expression of the optical lattice potential in eq. (13), we see that $\Delta \mu=4V_0\cos\theta$.

The density distribution function in momentum space [39] reads $n(\textbf{k})=\int \text{d}\textbf{r}\text{d}\textbf{r}'e^{i\textbf{k}\cdot(\textbf{r}-\textbf{r}')} \langle\psi^{\dagger}(\textbf{r})\psi(\textbf{r}')\rangle$ which can in turn be expressed in terms of the ${\hat a}_\textbf{k}=\frac1{N_S}\times \sum_i\hat{a}_ie^{-i\textbf{k}\cdot\textbf{r}_i}$ as $n(\textbf{k})=N_{S}|w(\textbf{k})|^{2}\langle \hat{a}_\textbf{k}^{\dagger}\hat{a}_\textbf{k}\rangle$ [40], or in other words, $n(\textbf{k})=N_{S}|w(\textbf{k})|^{2}\lim_{\tau\downarrow 0}G(\tau |0,\textbf{k})$, where $G(\tau |0,\textbf{k})\equiv\langle\hat{T}_{\tau}[\hat{a}^{\dagger}(\tau)_\textbf{k}\hat{a}(0)_\textbf{k}]\rangle$ is the Fourier transformation of the one-particle Green's function $G(\tau', j'|\tau, j) \equiv\langle\hat{T}_{\tau}[\hat{a}^{\dagger}_{j'}(\tau')\hat{a}_j(\tau)]\rangle$. By treating the hopping term in eq. (1) as perturbation, in the Dirac picture, the Green's function reads $G(\tau',j'\mid \tau,j)=\,\frac{\text{Tr}\{{{e^{-\beta \hat{H}_{0}}}\hat{T}_{\tau}[\hat{a}^{\dagger}_{j'}(\tau')\hat{a}_{j}(\tau)\hat{u}(\beta,0)]}\}} {\{\text{Tr}{e^{-\beta \hat{H}_{0}}}\hat{u}(\beta,0)\}}$, where $\hat{O}(\tau)=e^{\tau\hat{H}_{0}}\hat{O}e^{- \tau\hat{H}_{0}}$, and $\hat{u}(\beta,0)=\hat{T}_{\tau}[\exp(\int^{\beta}_{0}\text{d}\tau\sum_{\langle i,j \rangle} t\hat{a}^{\dagger}_{i}(\tau)\hat{a}_{j}(\tau))]$ is the evolution operator (setting $\hbar=1$) [41]. After expanding the evolution operator perturbatively, it is not difficult to see that the expansion of the Green's function consists of terms as $\frac{1}{n!}\sum_{i_{1},j_{1},\ldots,i_{n},j_{n}}t_{i_{1},j_{1}}\cdots t_{i_{n},j_{n}}\times\int^{\beta}_{0}\text{d}\tau_{1} \cdots \int^{\beta}_{0}\text{d}\tau_{n} \langle\hat{T}_{\tau}[\hat{a}^{\dagger}_{j'}(\tau')\hat{a}_{j}(\tau)\hat{a}^{\dagger}_{i_{1}}(\tau_{1}) \hat{a}_{j_{1}}(\tau_{1})\cdots\hat{a}^{\dagger}_{i_{n}}(\tau_{n})\times\hat{a}_{j_{n}}(\tau_{n})]\rangle_{0}$ with $\langle\;\rangle_0$ being the average quantity with respect to $\hat{H}_0$, in other words, n-particle Green's functions $G_n^{(0)}=\langle\hat{T}_{\tau}[\hat{a}^{\dagger}_{i_{1}}(\tau_{1}) \hat{a}_{j_{1}}(\tau_{1})\cdots\hat{a}^{\dagger}_{i_{n}}(\tau_{n})\hat{a}_{j_{n}}(\tau_{n})]\rangle_{0}$ with respect to $\hat{H}_0$ need to be calculated out, and this can be achieved by expanding them in terms of cumulants [27–29] $C^{(0)}_m(\tau'_{1},\ldots,\tau'_{m}|\tau_{1},\ldots\tau_{m}) =\langle \hat{T}_{\tau}[\hat{a}^{\dagger}(\tau'_{1})\times\hat{a}(\tau_{1})\cdots \hat{a}^{\dagger}(\tau'_{m})\hat{a}(\tau_{m})]\rangle_{0}$. However, since we are dealing with a system on a superlattice, according to the generalized Green's function method [36], the cumulants $C^{(0)}_{mA}=\langle \hat{T}_{\tau}[\hat{a}^{\dagger}_A(\tau'_{1})\hat{a}_A(\tau_{1})\cdots \hat{a}_A^{\dagger}(\tau'_{m})\hat{a}_A(\tau_{m})]\rangle_{0}$ on sublattice A and cumulants $C^{(0)}_{mB}=\langle \hat{T}_{\tau}[\hat{a}_B^{\dagger}(\tau'_{1})\hat{a}_B(\tau_{1})\cdots \hat{a}_B^{\dagger}(\tau'_{m})\hat{a}_B(\tau_{m})]\rangle_{0}$ on sublattice B are different and have to be treated separately.

In practice, the Green's functions can only be calculated out perturbatively by selecting specific groups of terms in the above cumulant expansion expressions in a proper way. According to the re-summed Green's function method [27,29,36], to the lowest order, only terms consisting of first-order cumulants are picked up. In terms of Matsubara frequency, the lowest-order re-summed one-particle Green's function in momentum space reads $\tilde{G}(\omega_{m},\textbf{k})=\sum_{l=0}^{+\infty}\!C^{(0)}_{1A} (\omega_{m})^{l+1} C^{(0)}_{1B}(\omega_{m})^{l}(t(\textbf{k}))^{2l} +\sum_{l=0}^{+\infty}\!C^{(0)}_{1A}(\omega_{m})^{l}C^{(0)}_{1B} (\omega_{m})^{l+1}(t(\textbf{k}))^{2l}\!+2\!\sum_{l=0}^{+\infty}\!C^{(0)}_{1A} (\omega_{m})^{l+1}\!\times C^{(0)}_{1B}(\omega_{m})^{l+\!1}(t(\textbf{k}))^{2l+\!1}$ where $C^{(0)}_{1A(B)}(\omega_{m})\!=\!\int^{\beta}_{0}C^{(0)}_{1A(B)}(\tau)\!\times e^{i\omega_{m}\tau} \text{d}\tau$ and $t(\textbf{k})\!=2t[\cos(k_{x}a)+\cos(k_{y}a)]$. Since $C^{(0)}_{1}\!(\tau)\!=\langle\, \hat{T}_{\tau}[\hat{a}^{\dagger}(\tau) \hat{a}(0)]\,\rangle_{0}$, by counting the detailed information of the eigenstates of $\hat{H}_0$, in the zero temperature limit, a complicated yet straightforward calculation leads to ${C^{(0)}_{1A}(\omega_{m})=\frac{n_{A}+1}{E_{n_{A}+1,n_{B}}-E_{n_{A},n_{B}}+i\omega_{m}} -\frac{n_{A}}{E_{n_{A},n_{B}}-E_{n_{A}-1,n_{B}}+i\omega_{m}}}$ and $C^{(0)}_{1B}(\omega_{m})=\frac{n_{B}+1}{E_{n_{A},n_{B}+1}-E_{n_{A},n_{B}}+i\omega_{m}} -\frac{n_{B}}{E_{n_{A},n_{B}}-E_{n_{A},n_{B}-1}+i\omega_{m}}$ together with the expression of the unperturbed ground-state energy E_nA,nB in eq. (3), the Green's function $G(\tau|0,\textbf{k})$ can then be calculated analytically via $G(\tau|0,\textbf{k})=\frac1{2\pi}\int_{-\infty}^{+\infty}\text{d}\omega_m \tilde{G}(\omega_m,\textbf{k})e^{-\omega_m\tau}$.

However, in order to calculate the one-particle Green's function quantitatively and to plot the corresponding TOF absorption pictures, some basic facts and parameters of possible experiments need to be clarified. We may propose that the superlattice Bose-Hubbard model may be realized in experiment by trapping an ultracold ⁸⁷Rb Bose gas in a cubic optical lattice created by laser beams with wavelength $\lambda\approx850\ \text{nm}$, the corresponding s-wave scattering length of ⁸⁷Rb is estimated to be about $a_{s}\approx5.34\ \text{nm}$. The superlattice structure is achieved, according to eq. (13), by tuning the phase difference θ. In order to show the alternate Mott and CDW phases, according to the discussion in the second section, we choose θ properly so that $\Delta \mu/U=0.5$. In the experiment, in order to eliminate the influence of the third dimension to reveal the 2D properties of the system, the lattice depth in the third dimension $V_3=30 E_R$ would be a reasonable choice [11,42,43].

By taking all the above information into account, we calculate and plot the TOF absorption pictures of the superlattice ultracold Bose system released from the CDW(2,1) state and from the MI(1,1) state for different lattice depths in fig. 5 and fig. 6, respectively, it is clear to see that the phases transit from the superfluid phase to the Mott or CDW phase when increasing the lattice depth.

Zoom In Zoom Out Reset image size

From figs. 5 and 6, we see that the quantum phase transitions of the CDW(2,1) case start to appear with deeper lattice depth when compared with what in the MI(1,1) case, in other words, the quantum critical value of the hopping parameter t_c for the CDW(2,1) case is smaller than in the MI(1,1) case, this is in consistence with the quantum phase diagram in figs. 2 and 4. When compared with the TOF pictures of the Bose-Hubbard system in a homogeneous square lattice [1], there are extra peaks at $(\pm \pi, \pm \pi)$ in the present TOF pictures, this is due to the bipartite superlattice structure, and reflects the inhomogeneity of the superfluid phase.

Zoom In Zoom Out Reset image size

The authors gratefully acknowledge Axel Pelster and Tao Wang for their stimulating and fruitful discussions. This Work was supported by National Natural Science Foundation of China under Grant No. 11275119 and by PhD Programs Foundation of Ministry of Education of China under Grant No. 20123108110004.