Quantum phases of dipolar bosons in multilayer optical lattice

We solve the model using the cluster mean-field theory with Gutzwiller ansatz  Fisher et al. (1989); Rokhsar and Kotliar (1991); Sheshadri et al. (1993); Bai et al. (2018); Pal et al. (2019); Bandyopadhyay et al. (2019); Suthar et al. (2020a); Malakar et al. (2020). For this, we subdivide the $\displaystyle K\times L\times M$ lattice system ($\displaystyle K$ and $\displaystyle L$ are the number of lattice sites along $\displaystyle x$ and $\displaystyle y$ directions) into $\displaystyle W$ small clusters of $\displaystyle k\times l\times m$ lattice sites, that is, $\displaystyle W=(K\times L\times M)/(k\times l\times m)$. Then, like in the site-decoupled mean-field theory, we can define a local Hamiltonian of the clusters Lühmann (2013b); Bai et al. (2018); Pal et al. (2019); Suthar et al. (2020a), and the total Hamiltonian is the sum of the cluster Hamiltonians.

In the present work, we limit ourselves to the case of bilayer ($\displaystyle M=2$) and trilayer ($\displaystyle M=3$) systems. The Hamiltonian of a cluster in the bilayer system is

	$\displaystyle\displaystyle\hat{H}_{C}$	$\displaystyle\displaystyle=$	$\displaystyle\displaystyle\sum_{r}\Bigg{\{}\sum_{p,q\in C}^{\prime}\Big{[}-J\Big{(}\hat{b}_{p+1,q,r}^{\dagger}\hat{b}_{p,q,r}+\hat{b}_{p,q+1,r}^{\dagger}\hat{b}_{p,q,r}$		(7)
			$\displaystyle\displaystyle+{\rm H.c.}\Big{)}+V\Big{(}\hat{n}_{p+1,q,r}\hat{n}_{p,q,r}+\hat{n}_{p,q+1,r}\hat{n}_{p,q,r}\Big{)}\Big{]}$
			$\displaystyle\displaystyle+\sum_{p,q\in\delta C}\Big{[}-J\Big{(}\phi_{p+1,q,r}^{*}\hat{b}_{p,q,r}+\phi_{p,q+1,r}^{*}\hat{b}_{p,q,r}+{\rm H.c.}\Big{)}$
			$\displaystyle\displaystyle+V\Big{(}\langle\hat{n}_{p+1,q,r}\rangle\hat{n}_{p,q,r}+\langle\hat{n}_{p,q+1,r}\rangle\hat{n}_{p,q,r}\Big{)}\Big{]}$
			$\displaystyle\displaystyle-\mu\sum_{p,q\in C}\hat{n}_{p,q,r}\Bigg{\}}-\sum_{p,q\in C}\Big{[}-J^{\prime}\Big{(}\hat{b}^{\dagger}_{p,q,1}\hat{b}_{p,q,2}+{\rm H.c.}\Big{)}$
			$\displaystyle\displaystyle+V^{\prime}\hat{n}_{p,q,1}\hat{n}_{p,q,2}\Big{]},$

where the prime in the summation denotes sum over lattice sites $\displaystyle(p,q)\in C$, such that, $\displaystyle(p+1,q)$ and $\displaystyle(p,q+1)\in C$, and $\displaystyle(p,q)\in\delta C$ denote the lattice sites at the boundary of the cluster $\displaystyle C$. The mean-field $\displaystyle\phi^{*}_{p+1,q,r}=\langle\hat{b}^{\dagger}_{p+1,q,r}\rangle$ and average occupancy $\displaystyle\langle\hat{n}_{p+1,q,r}\rangle$ with $\displaystyle(p+1,q)\notin C$ are computed at the boundary of the neighbouring cluster along $\displaystyle x$-direction, and are required to describe the intercluster hopping and NN interaction, respectively. Similarly, $\displaystyle\phi^{*}_{p,q+1,r}=\langle\hat{b}^{\dagger}_{p,q+1,r}\rangle$ and $\displaystyle\langle\hat{n}_{p,q+1,r}\rangle$ with $\displaystyle(p,q+1)\notin C$ are required to describe the intercluster hopping and NN interaction along $\displaystyle y$-direction between adjacent clusters. Thus, within a cluster the hopping and NN interaction terms are exact, and the intercluster hopping and NN interactions are accounted through the mean-fields and average occupancies, respectively. It is important to note that the interlayer hopping and NN interaction terms are exact in the cluster Hamiltonian. Now, like in the site-decoupled mean-field theory, the cluster Hamiltonian matrix can be calculated using the Fock basis states of the cluster $\displaystyle\{|n_{0}n_{1}...n_{m^{\prime}}\rangle\}$, where $\displaystyle m^{\prime}=(k\times l\times m)-1$. So, these basis states are the direct products of the local Fock states of the $\displaystyle(k\times l\times m)$ lattice sites within the cluster. By diagonalizing the Hamiltonian matrix, the ground state of the cluster can be obtained in the form

$$|\psi^{C}_{\alpha}\rangle=\sum_{n_{0}n_{1}...n_{m^{\prime}}}c^{(\alpha)}_{n_{0}n_{1}...n_{m^{\prime}}}|n_{0}n_{1}...n_{m^{\prime}}\rangle,$$

(8)

where $\displaystyle\alpha$ is the cluster index, and $\displaystyle c^{(\alpha)}_{n_{0}n_{1}...n_{m^{\prime}}}$ are the complex coefficients of the ground state $\displaystyle|\psi^{C}_{\alpha}\rangle$. Then, employing the Gutzwiller ansatz, the ground state of the entire lattice is

$$|\Psi_{\rm GW}\rangle=\prod_{\alpha=1}^{W}|\psi^{C}_{\alpha}\rangle=\prod_{\alpha=1}^{W}\sum_{n_{0}n_{1}...n_{m^{\prime}}}c^{(\alpha)}_{n_{0}n_{1}...n_{m^{\prime}}}|n_{0}n_{1}...n_{m^{\prime}}\rangle.$$

(9)

The local superfluid order parameter and average occupancy at the $\displaystyle(p,q)$ lattice site of the $\displaystyle r$th layer are

	$\displaystyle\displaystyle\phi_{p,q,r}$	$\displaystyle\displaystyle=$	$\displaystyle\displaystyle\langle\Psi_{\rm GW}|\hat{b}_{p,q,r}|\Psi_{\rm GW}\rangle,$
	$\displaystyle\displaystyle n_{p,q,r}$	$\displaystyle\displaystyle=$	$\displaystyle\displaystyle\langle\Psi_{\rm GW}|\hat{n}_{p,q,r}|\Psi_{\rm GW}\rangle.$		(10)

A relevant parameter of a quantum phase is the average occupancy per lattice site,

$$\rho=\frac{1}{(K\times L\times M)}\sum_{p=1,q=1,r=1}^{K,L,M}n_{p,q,r}.$$

(11)

In this work, we study quantum phases with the hard-core approximation. So, we have $\displaystyle\rho\leq 1$. An important quantity to probe the valence bond order in the bilayer system is the average pair density

$$\tilde{\rho}=\frac{1}{(K\times L)}\sum_{p,q=1}^{K,L}\langle\Psi_{\rm GW}|\hat{n}_{p,q,1}\hat{n}_{p,q,2}|\Psi_{\rm GW}\rangle.$$

(12)

For valence bond checkerboard solid phases $\displaystyle\tilde{\rho}$ is zero while $\displaystyle\rho$ is finite. The zero superfluid order parameter indicates the incompressibility of these solid phases and the checkerboard order in each layer is identified by the structure factor

$$S_{r}(\pi,\pi)=\frac{1}{K\times L}\sum_{\begin{subarray}{c}{p^{\prime},q^{\prime}}\\ p,q\end{subarray}}e^{{\rm i}\pi\{(p-p^{\prime})+(q-q^{\prime})\}}\langle\hat{n}_{p,q,r}\hat{n}_{p^{\prime},q^{\prime},r}\rangle.$$

(13)

In the phases with checkerboard order $\displaystyle S_{r}(\pi,\pi)$ is finite, while it is zero in the phases with uniform density distribution.

The starting point of our cluster mean field, hereafter, referred as cluster Gutzwiller mean field (CGMF) theory, is to choose an appropriate initial guess state $\displaystyle|\Psi_{\rm GW}\rangle$. From this we compute the initial $\displaystyle\phi_{p,q,r}$ and $\displaystyle n_{p,q,r}$. In general, we choose the same initial guess state $\displaystyle|\Psi^{C}_{\alpha}\rangle$ for all the $\displaystyle W$ clusters and consider $\displaystyle c^{(\alpha)}_{n_{0}n_{1}...n_{m^{\prime}}}=1/\sqrt{2^{m^{\prime}+1}}$. Then, using corresponding $\displaystyle\phi_{p,q,r}$ and $\displaystyle n_{p,q,r}$, we calculate the Hamiltonian matrix of a cluster given in Eq. (7), which we diagonalize Bai et al. (2018); Pal et al. (2019); Suthar et al. (2020a). We update the state $\displaystyle|\Psi_{\rm GW}\rangle$ using the new ground state of the cluster obtained from the diagonalization. Afterwards, we compute $\displaystyle\phi_{p,q,r}$ and $\displaystyle n_{p,q,r}$ using this updated $\displaystyle|\Psi_{\rm GW}\rangle$, and advance to the next cluster to repeat the same steps. We sweep the entire lattice system by continuing the procedure, and one such sweep constitutes an iteration. We continue the iterations until desired convergence in $\displaystyle\phi_{p,q,r}$ and $\displaystyle n_{p,q,r}$ is obtained. In our computations, we consider clusters ranging in size from $\displaystyle 1\times 1\times 2$ to $\displaystyle 4\times 4\times 2$ to tile lattice systems ranging in size from $\displaystyle 8\times 8\times 2$ to $\displaystyle 16\times 16\times 2$. To model an uniform infinite lattice system, we employ periodic boundary conditions in $\displaystyle\phi_{p,q,r}$ and $\displaystyle n_{p,q,r}$ along $\displaystyle x$ and $\displaystyle y$-directions. We also corroborate the stability of the obtained ground states with respect to different initial guess states having inhomogeneous distribution in $\displaystyle n_{p,q,r}$ and $\displaystyle\phi_{p,q,r}$. The initial guess states considered have checkerboard and random density patterns.

The system admits particle and hole vacuum states. And, these states correspond to $\displaystyle\rho=0$ and $\displaystyle\rho=1$, respectively. Using the dimer notation these states can be represented as

$$|\Psi\rangle_{\rm VAC}^{\rho=0}=\prod_{(p,q)}|00\rangle_{p,q},$$

(14)

and

$$|\Psi\rangle_{\rm VAC}^{\rho=1}=\prod_{(p,q)}|11\rangle_{p,q}.$$

(15)

Thus, in these states the system either has no atoms or has uniform distribution of one atom per lattice site throughout the system, respectively. The interlayer attractive NN interaction energetically favors the dimer state $\displaystyle|11\rangle_{p,q}$. In addition, the intralayer repulsive NN interaction can induce checkerboard density order. This interaction induced spatially periodic intralayer density modulation may be considered as two interpenetrating sublattices, $\displaystyle A$ and $\displaystyle B$. It is important to note that $\displaystyle n_{p,q,1}$ and $\displaystyle n_{p,q,2}$ have identical distributions or the checkerboard structure in both the layers are aligned. This is due to the attractive interlayer NN interaction, and is in contrast to the case when $\displaystyle V^{\prime}$ is repulsive. As mentioned earlier, the Hamiltonian of this system with $\displaystyle J^{\prime}=0$ is equivalent to the two species lattice Hamiltonian in 2D. In which case, $\displaystyle V^{\prime}$ is identified as the onsite interspecies interaction strength. For this two-species system with repulsive onsite interspecies interaction the checkerboard structure can have phases with spatially separated density or phase-separated Bai et al. (2020).

Due to the intralayer repulsive NN interaction, the bilayer system can exhibit a state of the dimers with solid or diagonal order, which is referred to as dimer checkerboard solid (DCBS). Then, considering the sublattice description of the checkerboard order, the DCBS state can be expressed as

$$|\Psi\rangle_{\rm DCBS}^{\rho=1/2}=\prod_{(p,q)\in A}|11\rangle_{p,q}\prod_{(p^{\prime},q^{\prime})\in B}|00\rangle_{p^{\prime},q^{\prime}}.$$

(16)

On the other hand, the interlayer hopping can stabilize the triplet state as in Eq. (5). As discussed before, this state with the checkerboard ordered density is the VCBS state. Depending on the average occupancy the VCBS states have the forms

$$|\Psi\rangle_{\rm VCBS}^{\rho=1/4}=\prod_{(p,q)\in A}|t_{0}\rangle_{p,q}\prod_{(p^{\prime},q^{\prime})\in B}|00\rangle_{p^{\prime},q^{\prime}},$$

(17)

and

$$|\Psi\rangle_{\rm VCBS}^{\rho=3/4}=\prod_{(p,q)\in A}|t_{0}\rangle_{p,q}\prod_{(p^{\prime},q^{\prime})\in B}|11\rangle_{p^{\prime},q^{\prime}}.$$

(18)

It is important to note that the states described in Eqs. (14)– (18) correspond to the incompressible phases of the bilayer system and the SF order parameter $\displaystyle\phi_{p,q,r}$ is zero in these phases. In addition, the checkerboard order in the DCBS and VCBS states are quantified by the nonvanishing $\displaystyle S_{r}(\pi,\pi)$. The average pair density $\displaystyle\tilde{\rho}$ given in Eq. (12) serves as an order parameter to distinguish the VCBS states from the DCBS state. Note that $\displaystyle\tilde{\rho}$ is zero in the VCBS states, but it is nonzero in the DCBS state. The system also exhibits supersolid and superfluid phases in which $\displaystyle\phi_{p,q,r}$ is finite. In the SF phase the system has uniform $\displaystyle n_{p,q,r}$ and $\displaystyle\phi_{p,q,r}$. In contrast, the supersolid phase has checkerboard order in $\displaystyle n_{p,q,r}$ and $\displaystyle\phi_{p,q,r}$. We show the schematic representation of the quantum phases discussed here, in Fig. 2. In the figure, the dimer in the DCBS phase is recognizable, and the blue shade over the bonds indicate the resonating structure of the VCBS states.

Like in the triplet state in the bilayer system, the dipolar atoms in a trilayer optical lattice can exhibit states of the form

$$|w_{0}\rangle_{p,q}=\frac{1}{\sqrt{3}}\left(|100\rangle_{p,q}+|010\rangle_{p,q}+|001\rangle_{p,q}\right),$$

(19)

and

$$|w_{1}\rangle_{p,q}=\frac{1}{\sqrt{3}}\left(|011\rangle_{p,q}+|101\rangle_{p,q}+|110\rangle_{p,q}\right).$$

(20)

These states have a resonating structure in one of the two sublattices, and are analogous to the VCBS states, given in Eq.(17) and Eq.(18), for bilayer optical lattice. They are stabilized by the interlayer hopping. They resemble the W-state of the three qubits that is studied in the context of the quantum information theory  Dür et al. (2000). The density of these states can vary as $\displaystyle\rho=1/6$, $\displaystyle 2/6$, $\displaystyle 4/6$ and $\displaystyle 5/6$ based on the filling in sublattice A and B. For the $\displaystyle|w_{0}\rangle$ state in sublattice A, the density of these states can be $\displaystyle\rho=1/6$ and $\displaystyle 4/6$ for the zero and unit filling in sublattice B, respectively. In a similar way, for the $\displaystyle|w_{1}\rangle$ state in sublattice A, the density can be $\displaystyle\rho=2/6$ and $\displaystyle 5/6$. On the other hand, the state corresponding to the density $\displaystyle\rho=3/6=1/2$ does not have a resonating structure unlike others, and is referred as the trimer checkerboard solid (TCBS) state. The form of this state, in terms of the sublattice description, is

$$|\Psi\rangle_{\rm TCBS}^{\rho=1/2}=\prod_{(p,q)\in A}|111\rangle_{p,q}\prod_{(p^{\prime},q^{\prime})\in B}|000\rangle_{p^{\prime},q^{\prime}}.$$

(21)

This is a generalization of the DCBS $\displaystyle\rho=1/2$ state, given in Eq.(16), for a bilayer optical lattice. It has a checkerboard pattern between the sublattice A and B, aligned between the three layers. The same occupancy between the three layers at a given lattice site is the result of the strong attractive interlayer interaction.

We calculate phase boundaries between the incompressible and compressible phases of the bilayer system using the mean-field theory. We do this by adapting the site-decoupling scheme, and perform perturbative analysis of the mean-field Hamiltonian. In this method we decompose the creation, annihilation and occupation number operators of a lattice site in terms of the local mean-fields and fluctuation operators as $\displaystyle\hat{b}_{p,q,r}=\phi_{p,q,r}+\delta\hat{b}_{p,q,r}$, $\displaystyle\hat{b}^{\dagger}_{p,q,r}=\phi^{*}_{p,q,r}+\delta\hat{b}^{\dagger}_{p,q,r}$, and $\displaystyle\hat{n}_{p,q,r}=n_{p,q,r}+\delta\hat{n}_{p,q,r}$ Then, the bilinear terms of the Hamiltonian in Eq. (1) are sum of linear terms of the local operators. In the perturbation analysis, we consider the interaction terms as the unperturbed Hamiltonian. These are diagonal with respect to the single site Fock basis states. The off diagonal hopping terms are considered as perturbations. The incompressible to compressible phase boundaries are, then, marked by the vanishing superfluid order parameter, that is $\displaystyle\phi_{p,q,r}\rightarrow 0^{+}$. The $\displaystyle\phi_{p,q,r}$ being small and associated with the off diagonal terms, it can be considered as the perturbation parameter.

To perform the first order perturbation analysis, we write the site-decoupled unperturbed Hamiltonian as

	$\displaystyle\displaystyle\hat{H}^{(0)}$	$\displaystyle\displaystyle=$	$\displaystyle\displaystyle\sum_{p,q,r}\hat{h}^{(0)}_{p,q,r}$		(22)
		$\displaystyle\displaystyle=$	$\displaystyle\displaystyle\sum_{p,q,r}\hat{n}_{p,q,r}\left(V^{\prime}n_{p,q,3-r}+V\overline{n}_{p,q,r}-\mu\right),$

where $\displaystyle\overline{n}_{p,q,r}=(n_{p+1,q,r}+n_{p-1,q,r}+n_{p,q+1,r}+n_{p,q-1,r})$, with $\displaystyle n_{p,q,r}=\langle\hat{n}_{p,q,r}\rangle$ representing the ground state occupancy at the site $\displaystyle(p,q,r)$. Now, consider the sublattice description of incompressible states with aligned checkerboard order, the unperturbed energy of the two different sublattices are

$$E^{A}_{n^{A}_{1}n^{A}_{2}}=V^{\prime}n^{A}_{1}n^{A}_{2}+\sum_{r=1}^{2}\left(4Vn^{A}_{r}n^{B}_{r}-\mu n^{A}_{r}\right),$$

(23)

for $\displaystyle(p,q)\in A$, and

$$E^{B}_{n^{B}_{1}n^{B}_{2}}=V^{\prime}n^{B}_{1}n^{B}_{2}+\sum_{r=1}^{2}\left(4Vn^{B}_{r}n^{A}_{r}-\mu n^{B}_{r}\right),$$

(24)

for $\displaystyle(p,q)\in B$. Here, $\displaystyle n^{A}_{r}$ and $\displaystyle n^{B}_{r}$ are the occupancies in the $\displaystyle r$th layer. It is important to note that $\displaystyle E^{A}_{00}=E^{B}_{00}=0$. As mentioned earlier, the hopping terms in the Hamiltonian are treated as the perturbations. Then, the site-decoupled perturbation Hamiltonians are

	$\displaystyle\displaystyle\hat{T}^{A}$	$\displaystyle\displaystyle=$	$\displaystyle\displaystyle-4J\sum_{r=1}^{2}\phi_{r}^{B}\left(\hat{b}^{A}_{r}+{\hat{b}^{A^{\dagger}}_{r}}\right)-J^{\prime}\sum_{r=1}^{2}\phi_{3-r}^{A}\left(\hat{b}^{A}_{r}+\hat{b}^{A^{\dagger}}_{r}\right),$
	$\displaystyle\displaystyle\hat{T}^{B}$	$\displaystyle\displaystyle=$	$\displaystyle\displaystyle-4J\sum_{r=1}^{2}\phi_{r}^{A}\left(\hat{b}^{B}_{r}+\hat{b}^{B^{\dagger}}_{r}\right)-J^{\prime}\sum_{r=1}^{2}\phi_{3-r}^{B}\left(\hat{b}^{B}_{r}+\hat{b}^{B^{\dagger}}_{r}\right),$

for the sublattice-$\displaystyle A$ and $\displaystyle B$, respectively. Now, to the first order in the superfluid order parameters, the perturbed ground state is

$$|\chi^{\alpha}\rangle=|n^{\alpha}_{1}n^{\alpha}_{2}\rangle+\sum_{\begin{subarray}{c}(m^{\alpha}_{1},m^{\alpha}_{2})\\ \neq(n^{\alpha}_{1},n^{\alpha}_{2})\end{subarray}}\frac{\langle m^{\alpha}_{1}m^{\alpha}_{2}|\hat{T}^{\alpha}|n^{\alpha}_{1}n^{\alpha}_{2}\rangle}{(E^{\alpha}_{n^{\alpha}_{1}n^{\alpha}_{2}}-E^{\alpha}_{m^{\alpha}_{1}m^{\alpha}_{2}})}|m^{\alpha}_{1}m^{\alpha}_{2}\rangle,$$

(25)

where $\displaystyle\alpha\in\{A,B\}$. Then, the superfluid order parameters are $\displaystyle\phi^{\alpha}_{\kappa}=\langle\chi^{\alpha}|\hat{b}^{\alpha}_{\kappa}|\chi^{\alpha}\rangle$. In our study, we have considered that the intralayer parameters are same for both the layers–the layers are identical. Therefore, $\displaystyle\phi^{A}_{1}=\phi^{A}_{2}=\varphi^{A}$ and $\displaystyle\phi^{B}_{1}=\phi^{B}_{2}=\varphi^{B}$.

We first obtain the phase boundary separating the DCBS phase from a compressible phase, which can either be superfluid or supersolid. The DCBS state in Eq. (16) is an eigenstate of $\displaystyle\hat{H}^{(0)}$ in Eq. (22), with occupancies $\displaystyle(n^{A}_{1},n^{A}_{2})=(1,1)$ and $\displaystyle(n^{B}_{1},n^{B}_{2})=(0,0)$. Then, the superfluid order parameters calculated with respect to the perturbed ground state in Eq. (33) are $\displaystyle(p,q)\in A$

$$\varphi^{A}=-\frac{4J\varphi^{B}}{(E^{A}_{11}-E^{A}_{01})+J^{\prime}},$$

(26)

and

$$\varphi^{B}=-\frac{4J\varphi^{A}}{(E^{B}_{00}-E^{B}_{10})+J^{\prime}}.$$

(27)

We solve Eqs. (26) and  (27) simultaneously, and then, we take the limit $\displaystyle\{\varphi^{A},\varphi^{B}\}\rightarrow 0^{+}$ to get

$$16J^{2}=\left[(E^{A}_{11}-E^{A}_{01})+J^{\prime}\right]\times\left[(E^{B}_{00}-E^{B}_{10})+J^{\prime}\right].$$

(28)

Now, substituting the values of $\displaystyle E^{A}_{11}$, $\displaystyle E^{A}_{01}$ from Eq. (23), and $\displaystyle E^{B}_{10}$ from Eq. (24), we obtain the DCBS phase boundary as a solution of

$$16J^{2}=(V^{\prime}-\mu+J^{\prime})(\mu-4V+J^{\prime}).$$

(29)

From this, the DCBS lobe in the plane of $\displaystyle J/V-\mu/V$ can be obtained for different values of $\displaystyle J^{\prime}$ and $\displaystyle V^{\prime}$.

Next, we calculate the equations of the phase boundaries separating the parameter domains of the VCBS states described in Eqs. (17) and  (18) from the compressible phases of the system. It is important to notice that these two states are eigenstates of $\displaystyle\hat{H}^{(0)}$ in Eq. (22). But, as emphasized earlier, the interlayer hopping is essential to stabilize the $\displaystyle|t_{0}\rangle_{p,q}$ triplet state of the VCBS states. And, this term is not present in the site-decoupled unperturbed Hamiltonian $\displaystyle\hat{H}^{(0)}$. So, in order to obtain the triplet state $\displaystyle|t_{0}\rangle$ as an eigenstate, we define local unperturbed Hamiltonian of the sublattice-$\displaystyle A$ as

	$\displaystyle\displaystyle\mathcal{H}^{(0)}_{A}$	$\displaystyle\displaystyle=$	$\displaystyle\displaystyle-J^{\prime}(\hat{b}^{\dagger}_{p,q,1}\hat{b}_{p,q,2}+{\rm H.c.})+V^{\prime}\hat{n}_{p,q,1}\hat{n}_{p,q,2}$		(30)
			$\displaystyle\displaystyle+\sum_{r=1}^{2}\left(V\hat{n}_{p,q,r}\overline{n}_{p,q,r}-\mu\hat{n}_{p,q,r}\right),$

for $\displaystyle(p,q)\in A$. Then, we can express the unperturbed ground state energy for the sublattice-A as

$$E^{A}_{t_{0}}=\langle t_{0}|\mathcal{H}^{(0)}_{A}|t_{0}\rangle=-J^{\prime}+\sum_{r=1}^{2}\left(4Vn^{A}_{r}n^{B}_{r}-\mu n^{A}_{r}\right).$$

(31)

This has a contribution from interlayer hopping process, and in contrast to the Eq. (23), it does not have interlayer interaction energy. This is because the pair expectation, as mentioned earlier, is zero with respect to $\displaystyle|t_{0}\rangle$. It is to be noted that in Eq. (31), $\displaystyle n^{A}_{1}=n^{A}_{2}=0.5$, which are the occupancies calculated with respect to the $\displaystyle|t_{0}\rangle$ state. On the other hand, unlike in the previous case, the perturbation Hamiltonian now contains only the intralayer hopping terms, that is,

$$\hat{T^{\prime}}^{A}=-4J\sum_{r=1}^{2}\phi_{r}^{B}(\hat{b}^{A}_{r}+\hat{b}^{A^{\dagger}}_{r}).$$

(32)

Then, the perturbed ground state is

$$|\chi^{A}_{t_{0}}\rangle=|t_{0}\rangle+\sum_{m^{A}_{1},m^{A}_{2}}\frac{\langle m^{A}_{1}m^{A}_{2}|\hat{T^{\prime}}^{A}|t_{0}\rangle}{(E^{A}_{t_{0}}-E^{A}_{m^{A}_{1}m^{A}_{2}})}|m^{A}_{1}m^{A}_{2}\rangle,$$

(33)

where $\displaystyle(m^{A}_{1},m^{A}_{2})\in\{(0,0),(1,1)\}$. Using this state the superfluid order parameters in the sublattice-$\displaystyle A$ can be calculated as $\displaystyle\phi^{A}_{r}=\langle\chi^{A}_{t_{0}}|\hat{b}^{A}_{r}|\chi^{A}_{t_{0}}\rangle$. Similar to the previous case $\displaystyle\phi^{A}_{1}=\phi^{A}_{2}=\varphi^{A}$. Then, the superfluid order parameter is

$$\varphi^{A}=-4J\phi^{B}\left[\frac{1}{(E^{A}_{t_{0}}-E^{A}_{00})}+\frac{1}{(E^{A}_{t_{0}}-E^{A}_{11})}\right].$$

(34)

For the VCBS state with $\displaystyle\rho=1/4$, we calculate the superfluid order parameter as given in Eq. (34) from the perturbative correction to $\displaystyle|t_{0}\rangle_{p,q}$ in sublattice-$\displaystyle A$. But, we perform the perturbation analysis of the site-decoupled Hamiltonian in sublattice-$\displaystyle B$. So, we obtain the superfluid order parameter as given in Eq. (27) from the perturbative correction to $\displaystyle|00\rangle_{p,q}$ state in sublattice-$\displaystyle B$. We solve the Eqs. (34) and  (27) simultaneously, and then take the limit $\displaystyle\{\varphi^{A},\varphi^{B}\}\rightarrow 0^{+}$ to obtain

	$\displaystyle\displaystyle\frac{1}{16J^{2}}=$		$\displaystyle\displaystyle\left[\frac{1}{(E^{A}_{t_{0}}-E^{A}_{00})}+\frac{1}{(E^{A}_{t_{0}}-E^{A}_{11})}\right]$		(35)
			$\displaystyle\displaystyle\times\left[\frac{1}{(E^{B}_{00}-E^{B}_{10})+J^{\prime}}\right].$

Note that, for the unperturbed ground state the occupancies $\displaystyle(n^{A}_{1},n^{A}_{2})=(0.5,0.5)$ and $\displaystyle(n^{B}_{1},n^{B}_{2})=(0,0)$. Now, by substituting the values of $\displaystyle E^{A}_{t_{0}}$ from Eq. (31), $\displaystyle E^{A}_{11}$ from Eq. (31), and $\displaystyle E^{B}_{10}$ from Eq. (31), we obtain the VCBS ($\displaystyle\rho=1/4$) phase boundary as a solution of

$$16J^{2}(2J^{\prime}+V^{\prime})=(\mu-2V+J^{\prime})(J^{\prime}+\mu)(\mu-J^{\prime}-V^{\prime}).$$

(36)

It is worth mentioning that $\displaystyle E^{A}_{11}$ from Eq. (31) is equal to $\displaystyle\langle 11|\mathcal{H}^{(0)}_{A}|11\rangle$. Similarly, for the VCBS state with $\displaystyle\rho=3/4$, we perform the perturbative analysis to state $\displaystyle|11\rangle_{p,q}$ in the sublattice-$\displaystyle B$. We obtain the expression of the superfluid order parameter $\displaystyle\varphi^{B}$, which is similar to the Eq. (27) with the superscripts $\displaystyle A$ and $\displaystyle B$ interchanged. We, then simultaneously solve this equation and Eq. (34), and take the limit $\displaystyle\{\varphi^{A},\varphi^{B}\}\rightarrow 0^{+}$ like in the previous case. Then, we obtain the VCBS ($\displaystyle\rho=3/4$) phase boundary as a solution of

	$\displaystyle\displaystyle 16J^{2}(2J^{\prime}+V^{\prime})=$		$\displaystyle\displaystyle(2V+V^{\prime}-\mu+J^{\prime})(\mu-4V+J^{\prime})$		(37)
			$\displaystyle\displaystyle\times(\mu-J^{\prime}-4V-V^{\prime}).$

From Eqs. (36) and (37) the VCBS lobes with $\displaystyle\rho=1/4$ and $\displaystyle 3/4$ can be obtained in the $\displaystyle J/V-\mu/V$ plane for different values of $\displaystyle J^{\prime}$ and $\displaystyle V^{\prime}$.

The above formalism can be generalized to the trilayer system and the details of the derivation are given in the Appendix B. The equation defining the phase boundary between the TCBS phase and the compressible phase is

$$16J^{2}=(2V^{\prime}-\mu+2J^{\prime})(\mu-4V+2J^{\prime}).$$

(38)

Based on this, the VCBS $\displaystyle(\rho=1/6)$ phase to compressible phase boundary is given by

$$16J^{2}(3V^{\prime}+\mu+8J^{\prime})=(3\mu-4V+6J^{\prime})(\mu+2J^{\prime})(\mu-V^{\prime}),$$

(39)

and that of VCBS $\displaystyle(\rho=2/6)$ is

$$16J^{2}(5V^{\prime}-\mu+8J^{\prime})=(3\mu-8V+6J^{\prime})(\mu-2V^{\prime}-2J^{\prime})(\mu-V^{\prime}).$$

(40)

Invoking the particle-hole symmetry of the model, we can write the phase boundaries between the VCBS $\displaystyle\rho=4/6$ and the compressible phase as

	$\displaystyle\displaystyle 16J^{2}$		$\displaystyle\displaystyle\left(3V^{\prime}-4V+\mu+8J^{\prime}\right)=\left(6J^{\prime}+4V+6V^{\prime}-3\mu\right)$		(41)
			$\displaystyle\displaystyle\times\left(4V-\mu-2J^{\prime}\right)\left(-\mu+V^{\prime}+4V\right),$

and that of VCBS $\displaystyle\rho=5/6$ as

	$\displaystyle\displaystyle 16J^{2}$		$\displaystyle\displaystyle\left(5V^{\prime}+4V-\mu+8J^{\prime}\right)=\left(6J^{\prime}+8V+6V^{\prime}-3\mu\right)$		(42)
			$\displaystyle\displaystyle\times\left(4V-\mu+2J^{\prime}+2V^{\prime}\right)\left(-\mu+V^{\prime}+4V\right).$

The phase diagram obtained based on these equations, shown in the results section, is in good agreement with the one obtained numerically.

This is a suitable choice to probe the interplay between different hopping and interaction energies theoretically. In terms of lattice constants, it corresponds to $\displaystyle a_{z}=\sqrt[3]{2}a$. We extend the theoretical insights gained from this case to the $\displaystyle|V^{\prime}|\neq V$ regime.