We consider a system of dipolar bosons confined to a 2D square optical lattice, as illustrated in figure 1, under the influence of both DDI and TBI. The lattice lies in the $x$-$y$ plane with a lattice constant $a$. The spheres denote atoms, while the arrows indicate the polarization direction, which is fixed within the $y$-$z$ plane. In the grand canonical ensemble, this system is described by an extended Bose-Hubbard model (eBHM) with TBI, and the Hamiltonian is where $p(q)$ is the lattice coordinate along the $x(y)$ direction, $\hat{b}_{p,q}^{\dagger} (\hat{b}_{p,q})$ is the bosonic creation (annihilation) operator, and $\hat{n}_{p,q}$ is the particle number operator for bosons at the lattice site. $J_x(J_y)$ is the tunneling amplitude between adjacent lattice sites along the $x(y)$ direction, and $V_x(V_y)$ is the DDI strength between nearest-neighbor atoms along the $x(y)$ direction, satisfy the relation $\ V_y = V_x (1-3 \sin^2 \theta) $ [13]. The polarization angle $\theta$ is defined as the angle between the atomic polarization direction and the $z$-axis. At $\theta = 0$, the dipole interactions in the $x$- and $y$-directions are identical. As $\theta$ gradually increases to $90^\circ$, the dipole interaction along the $x$-direction remains unchanged, while that along the $y$-direction changes from repulsive to attractive. Additionally, $U$ is the on-site interaction strength, $W$ is the TBI strength, and $\mu$ is the chemical potential.

To investigate the phase transitions of the system, we use the inhomogeneous Gutzwiller mean-field approximation to decouple the lattice [31, 37-39]. This method neglects inter-site spatial correlations and quantum fluctuations [39, 42], which results in an overestimation of ordered-phase stability and an underestimation of the temperature-driven NF phase. Nevertheless, it can capture key features of the eBHM and quantum phase transitions in dipolar quantum gases within optical lattices, and is expected to retain the overall qualitative features of the phase transitions both at zero and finite temperatures [9, 33, 43]. Within this framework, the bosonic creation and annihilation operators are expressed as $\hat{b}_{p,q}^{\dagger} = \phi_{p,q}^* + \delta \hat{b}_{p,q}^{\dagger}$ and $\hat{b}_{p,q} = \phi_{p,q} + \delta \hat{b}_{p,q}$, where $\phi_{p,q}^* = \langle \hat{b}_{p,q}^{\dagger} \rangle$ and $\phi_{p,q} = \langle \hat{b}_{p,q} \rangle$ are quantum mechanical expectation values for particle creation and annihilation at lattice site $(p,q)$, $\delta \hat{b}_{p,q}^{\dagger}$ and $\delta \hat{b}_{p,q}$ are fluctuation operator. Accordingly, the operator product of the creation and annihilation operators appearing in the hopping term can be written as and the DDI term is expressed as where ${n}_{p,q} = \langle \hat{n}_{p,q} \rangle$ is the average particle number at lattice site $(p,q)$, and lattice point$(p^{^{\prime}}, q^{^{\prime}})$ denotes a nearest-neighbor lattice sites of $(p,q)$. Using equations (2) and (3), the Hamiltonian in equation (1) can be written as where $\hat{H}_{p,q}^{\textrm{MF}}$ is the single-site mean-field Hamiltonian of site $(p,q)$, and is given by

Under the mean-field approximation, the coupling between adjacent lattice sites is mediated by the parameters $\phi_{p,q}$ and $n_{p,q}$. Consequently, the eigenstate of the entire lattice takes the form of a direct product of individual single-site states and can be written as where $\left| \psi \right\rangle_{p,q}$ denotes the ground-state wave function of the lattice site $(p,q)$, $N_b$ is the maximum number of bosons per site, and $\left| n \right\rangle_{p,q}$ corresponds to the Fock state of $n$ bosons in the occupation number representation at the lattice site $(p,q)$. The coefficient $c_n^{(p,q)}$ serves as the normalization factor for the lattice site $(p,q)$ and satisfies the condition $\sum_{n = 0}^{N_b}|c_n^{(p,q)}|^2 = 1$. Therefore, the expectation value of the annihilation operator at each lattice site $(p,q)$ is given by Similarly, the average particle number at site $(p,q)$ is The ground-state energy of the system can be determined numerically through a self-consistent procedure. First, a set of $(\phi_{p,q}, n_{p,q})$ is randomly assigned as initial values to each lattice site in the system. Following the diagonalization of the Hamiltonian in equation (5) within the particle number representation, the ground-state energy and corresponding wave function are subsequently calculated for each lattice site. Utilizing this ground-state wave function, updated expectation values for the annihilation operators and the particle number are derived via equations (7) and (8), yielding a new set $(\phi_{p,q}^{^{^{\prime}}}, n_{p,q}^{^{^{\prime}}})$. Then, these updated values $(\phi_{p,q}^{^{^{\prime}}}, n_{p,q}^{^{^{\prime}}})$ are iteratively reintroduced into equation (5) as initial conditions for further computation. The iterative procedure continues until convergence is achieved by monitoring the ground-state energy, thereby ensuring minimization of the system's total energy. Ultimately, this process yields the ground state of the system, along with the corresponding expectation values of the annihilation operator and the particle number at each lattice site.

At finite temperature, utilizing the complete eigenvalue spectrum $E^{l}_{p,q}$ and the corresponding eigenfunctions $\left| \psi\right\rangle^{l}_{p,q}$, the thermal average of $\phi_{p,q}$ in the mean-field approximation is and the thermal average of the particle number density $n_{p,q}$ is where $Z = \sum_{l = 0}^{N_b} \mathrm{e}^{-\beta E^{l}} $ is the single-site partition function, $\beta = 1/k_{\mathrm{B}}T$, $k_{\mathrm{B}}$ denotes the Boltzmann constant, and $T$ represents the temperature. In this study, temperature is defined in a dimensionless form, as energy is scaled relative to $U$, and the system is characterized by a uniform tunneling amplitude $J_x = J_y = J$.

The order parameter characterizes the degree of order or symmetry breaking in a system undergoing a phase transition. In our system, different quantum phases can be distinguished by combining the order parameter, the mean particle number $n$, the expectation value of the annihilation operator $\phi$, the density structure factor $S(\vec{k})$ [13] and the local density compressibility $\kappa$ [44]. The average occupancy per lattice site is The average expectation value of the annihilation operator is in which $N = L \times L$ denotes the size of the system. In the MI and DW phases, $\phi = 0$, $n_{p,q}$ takes integer values, and $n$ remains constant as the chemical potential is varied-indicating incompressibility. The SS and SF phases are compressible, with $\phi \gt 0$ and $n_{p,q}$ taking non-integer (real) values.

The quantum phase with periodic density modulation can be characterized by the density structure factor [13] where $\vec{k} \equiv \left( k_{x}, k_{y} \right)$ is a reciprocal lattice vector. Analysis of structure factors at characteristic wavevectors enables discrimination between checkerboard-type and striped-type phases (e.g. CBDW and CBSS phases exhibit a finite peak at $S(\pi, \pi) \gt 0$, whereas SDW and SSS phases are characterized by $S(\pi, 0) \gt 0$). Moreover, the DW phases are distinguished from the SS phases by the superfluid order parameter, with $\phi = 0$ for both CBDW and SDW phases, and $\phi \gt 0$ for CBSS and SSS phases.

At finite temperature, the NF phase can be distinguished from incompressible phases by examining the local variance of the system's density, which serves as a measure of local compressibility [44] In the MI and DW phases, $\kappa = 0$, while in the NF phase, $\kappa\neq 0$. The phase boundaries between different phases at zero temperature and at finite temperature can be determined through the above order parameters. It should be noted that at finite temperatures, the order parameter must be computed as a thermal average. The characteristic order parameters for each phase in the system under study are summarized in table 1.

Landau theory provides an effective method for studying phase transitions between compressible and incompressible quantum states [13, 40, 41]. We employ this approach to examine the phase transition of bosons in a 2D optical lattice with DDI and TBI. In line with the decoupling scheme outlined in section 2.2, the hopping and interaction terms are decoupled using the order parameters $\phi_{p,q}$ and $ n_{p,q}$. As the order parameter $\phi_{p,q}$ changes from zero to a finite value at the critical point, the hopping term can be treated as a perturbation. Through second-order perturbative analysis, the order parameter at the phase boundary between the incompressible and compressible phases satisfies the following relation where $\widetilde{\mu} = \mu - V_x \sum_{p^{^{\prime}}} \ n_{p^{^{\prime}},q}- V_y \sum_{q^{^{\prime}}} \ n_{p,q^{^{\prime}}} $ and $\overline{\phi}_{p,q} = \sum_{p^{^{\prime}},q^{^{\prime}}}$ ${\phi}_{p^{^{\prime}},q^{^{\prime}}}$, $(p^{^{\prime}}, q^{^{\prime}})$ denotes the nearest-neighbor lattice sites of $(p,q)$.

To determine the phase transition between the DW and SS phases, we divide the system into two nested sublattices, $A$ and $B$, which are coupled via nearest-neighbor hopping and interactions. Each sublattice is characterized by distinct occupation numbers, $n_A$ and $n_B$, and their corresponding annihilation operator averages, $\phi_A$ and $\phi_B$. The phase boundary separating the SDW and SSS phases is defined by the conditions: $\overline{\phi}_{p,q} = 2\left( \phi_{A} + \phi_{B} \right)$ and $\widetilde{\mu}_{A} = \mu -2 \left( V_{x}n_{B} + V_{y}n_{A} \right)$ for sites $\left( p, q \right)$ belonging to sublattice $A$, and $\overline{\phi}_{p,q} = 2\left( \phi_{A} + \phi_{B} \right)$ and $\widetilde{\mu}_{B} = \mu - 2\left( V_{x}n_{A} + V_{y}n_{B} \right)$ for sites $\left( p, q \right)$ in sublattice $B$. Consequently, the expectation values of the annihilation operators, $\phi_A$ and $\phi_B$, satisfy the following relations

By simultaneously solving the two equations above, the phase boundary between SDW and SSS phases (or between MI and SF phases) can be determined when $n_A \neq n_B$, the preceding equation defines the phase boundary separating the SDW and SSS phases in the presence of TBI. In the case $n_A = n_B$, it simplifies to the transition boundary between the MI and SF phases [13]. Through an analogous derivation, the phase boundary between CBDW and CBSS phases (or between MI and SF phases) is obtained as

In this section, we focus on the role of TBI and investigate its impact on the ground-state phases in different polarization directions. The distributions of $n_{p,q}$ and $\phi_{p,q}$ can reveal the characteristics of different phases in the system. In the CBDW and CBSS phases, $n_{p,q}$ exhibits a checkerboard spatial modulation, whereas in the SDW and SSS phases, it exhibits a stripe order. Figure 2(a) illustrates the CBDW(1,0) phase, where the particle density at all lattice sites is quantized to integer values and exhibits periodic modulation along both the $x$- and $y$-directions, while $\phi_{p,q}$ remains zero at all lattice sites. Figure 2(b) presents the CBSS phase, where the particle density remains a checkerboard modulated but adopts non-integer real values, and $\phi_{p,q}$ acquires finite real values at all lattice sites. The SDW(1,0) phase, shown in figure 2(c), exhibits a stripe-ordered density distribution along the $y$-direction, with integer filling at all sites and no modulation along the $x$-direction. Additionally, the order parameter $\phi_{p,q}$ also remains zero across the lattice. Figure 2(d) corresponds to the SSS phase, in which both $n_{p,q}$ and $\phi_{p,q}$ display stripe-like spatial modulation, with $\phi_{p,q}$ taking positive values at every site.

To further characterize the distinct phase transitions by means of the order parameters $n$, $\phi$, $S(\pi, \pi)$, and $S(\pi,0)$, figure 3 shows the variation of these order parameters with $J/U$ and $\mu/U$. As shown in figure 3(a), at $\mu = 0.5U$ and $V_x = 2V_y = 0.8U$, the system undergoes a transition from the CBDW phase to the CBSS phase with increasing $J/U$. In the CBDW phase, the order parameters are characterized by $\phi = 0$, commensurate filling ($n \in \mathbb{N}$), $S(\pi, \pi) \gt 0$, and $S(\pi,0) = 0$. For instance, the CBDW(1,0) phase is manifested as $n = 0.5$ and $S(\pi, \pi) = 0.25$. In the CBSS phase, $\phi$ becomes finite, $n$ monotonically increases with $J/U$, $S(\pi, \pi)$ remains positive, and $S(\pi,0)$ stays zero. Upon further increasing $J/U$, the system transitions into the SF phase, where $\phi \gt 0$, $n$ increases with $J/U$, and both $S(\pi, \pi)$ and $S(\pi,0)$ vanish. Figure 3(b) displays the behavior of the order parameters at $J = 0.01U$ and $V_x = 2V_y = 0.8U$. Here, the average particle density $n$ exhibits three plateaus at $n = 0.5$, $1.0$, and $1.5$, while $S(\pi, \pi) \gt 0$ and $S(\pi, 0) = 0$. These plateaus correspond to the CBDW(1,0), CBDW(2,0), and CBDW(3,0) phases, respectively. Within these parameter regimes, the system density remains constant against variations in the chemical potential, indicating the characteristic of an incompressible phase.

Figure 3(c) presents the variation of order parameters as a function of $J/U$ for $V_x = -2V_y = 0.2U$ and $\mu = 0.1U$. As $J/U$ increases, the system undergoes a phase transition from the SDW phase to the SSS phase. In contrast to the CBDW and CBSS phases, the stripe-ordered phases (SDW and SSS) are characterized by a vanishing order parameter $S(\pi, \pi) = 0$ and $S(\pi,0) \gt 0$. Specifically, the SDW phase is characterized by $\phi = 0$, $n \in \mathbb{N}$, $S(\pi, \pi) = 0$, and $S(\pi, 0) \gt 0$. For instance, the SDW(1,0) phase exhibits $n = 0.5$ and $S(\pi, 0) = 0.25$. Conversely, the SSS phase features $\phi \gt 0$, $n$ increases monotonically with $J/U$, $S(\pi, \pi) = 0$, and $S(\pi, 0) \gt 0$. With a further increase in $J/U$, the system transforms from the SSS phase to the SF phase. Figure 3(d) shows the variation of order parameters as a function of $\mu/U$ for $V_x = -2V_y = 0.2U$ and $J = 0.03U$. Increasing $\mu/U$ drives phase transitions among the SDW, MI, SF, and SSS phases, with the density $n$ exhibits four plateaus at $n$ = $0.5$, $1.0$, $2.0$, and $3.0$. In the $n = 0.5$ regime, the system exhibits SDW order, characterized by $S(\pi, \pi) = 0$ and $S(\pi,0) \gt 0$. At densities $n$ = 1.0, 2.0, 3.0, both $S(\pi, \pi)$ and $S(\pi,0)$ vanish, indicating a MI phase. In regions where $n$ increases continuously with $\mu/U$, the system enters an SSS phase when $S(\pi, \pi) = 0$ and $S(\pi,0) \gt 0$, and transitions to an SF phase with $S(\pi, \pi) = 0$ and $S(\pi,0) = 0$.

Based on the above information, we present the ground-state phase diagram in figure 4 for the polarization angle $\theta = 24.09^\circ$ and dipole interactions $V_x = 2V_y = 0.8U$ in the $\mu/U$-$J/U$ plane. At this polar angle, the DDI varies along the $x$ and $y$ directions, but both are repulsive. In the absence of TBI (figure 4(a)), the phase diagram exhibits three distinct quantum phases: CBDW, CBSS, and SF, consistent with earlier findings reported in [18]. For repulsive TBI, as shown in figures 4(b) and (c), CBDW lobes expand while the CBSS phase region contracts along the $\mu$-axis, leading to the collapse of the CBSS-SF boundary into the CBDW regime. This indicates that, under repulsive TBI, the CBDW phase exhibits enhanced stability relative to the CBSS phase. Notably, attractive TBI exerts a distinct influence on the phase diagram, substantially enhancing the stability of the CBSS phase compared to the $W = 0$ case. As illustrated in figure 4(d), the system exhibits a contraction of higher-density CBDW phases, such as CBDW(3,0) and CBDW(4,0), and a toward the emergence of CBDW(5,0). Concurrently, the CBSS phase displays markedly enhanced stability in the vicinity of the CBDW(3,0) regime.

We now discuss the ground-state phase diagram for $\theta = 45^\circ$ and dipole interactions $V_x = -2V_y = 0.2U$. At this polar angle, the DDI is repulsive along the $x$-direction and attractive along the $y$-direction. Comparison with the $\theta = 24.09^\circ$ case, figure 5 shows that the CBDW and CBSS phases transform into the SDW and SSS phases, respectively, and that the MI phase emerges between the SDW lobes at small hopping parameters. As shown in figures 5(a)-(c), increasing the repulsive TBI enlarges the MI (e.g. MI(2) and MI(3)) and SDW (e.g. SDW(2,1) and SDW(3,2)) lobes, while the SSS phase region shrinks. This is because repulsive TBI enhance the effective on-site interaction, thereby stabilizing the MI and SDW phases and suppressing the SSS phase. In contrast, weakly attractive TBI dramatically enhance the stability region of the SSS phase. As observed in figure 5(d), the SSS phase expands near the SDW phase, and the SSS to SF phase boundary shifts toward larger hopping. Compared with $W = 0$ case, the SDW phases (e.g. SDW(3,2) and SDW(4,3)) exhibit relatively weaker suppression of their hopping amplitude. For instance, the tip of the SDW(3,2) phase shifts from $J/U = 0.017$ to $0.015$, while that of the SDW(4,3) lobe shifts from $J/U = 0.013$ to $0.010$. This is because the attractive TBI partially screen the interatomic interactions, thereby suppressing the higher-density interstitial phases. A comparison between figures 5(a) and (d) further indicates that the tip of the SSS region (surrounding the SDW(3,2) lobe) shifts from $J/U$ = 0.021 to 0.023. The width of the SSS region also increases, such as at $J/U = 0.02$, the lower boundary changes from $\mu/U = 2.37$ to $2.24$, while the upper boundary shifts from $\mu/U = 2.6$ to $2.49$. The expansion of the SSS phase originates from attractive TBI suppressing interatomic repulsion, which increases the mobility of DW phase, thereby enhancing SSS phase characteristics. Comparing the phase boundaries between the incompressible and compressible phases obtained by numerical calculations and second-order perturbation theory in figures 4 and 5, we observe excellent agreement between the two approaches. It is notable that both repulsive and attractive TBI can be engineered and controlled among bosons [21, 23, 47, 48], thereby preventing rapid destruction from three-body losses and attractive collapse [24, 28, 45]. This enables the experimental observation of quantum phenomena arising under such repulsive and attractive TBI.

In this section, we turn our attention to the effects of thermal fluctuations on quantum phases stabilized by anisotropic DDI and TBI. Figure 6 presents the finite-temperature phase diagram at a fixed polarization angle of $\theta = 24.09^\circ$ and dipole interactions $V_x = 2V_y = 0.8U$. As shown in figures 6(a)-(c), at low temperature ($k_{\mathrm{B}}T = 0.02U$), all the zero-temperature phases are retained, with an additional NF phase emerging between the CBDW lobes (green shaded regions). The CBSS phase remains stable at finite temperatures in the presence of either repulsive (figure 6(b)) or attractive (figure 6(c)) TBI. Upon increasing temperature, the CBDW phase loses its crystalline order and fully transforms into the NF phase, whereas the CBSS phase remains stable, as shown in figures 6(d)-(f) at $k_{\mathrm{B}}T = 0.1U$. Notably, repulsive TBI substantially reduces the CBSS regions situated between the CBDW and SF phases, while the NF phase occupies an expanded region of the phase diagram. For instance, at fixed chemical potential ($\mu/U = 3.5$) and TBI strength ($W = 0.5U$), the critical hopping of the NF phase extends from $J/U = 0.14$ to $0.17$. In contrast, under attractive TBI, the NF phase is suppressed and the CBSS phase is stabilized. At the chemical potential ($\mu/U = 3.5$) and TBI strength ($W = -0.1U$), the critical hopping value for the NF phase decreases from $J/U = 0.14$ to $0.11$. As the temperature is further increased to $k_{\mathrm{B}}T = 0.3U$ (figures 6(g)-(i)), the width of the NF region expands and shifts toward larger values of the hopping parameter, while the CBSS phase region shrinks. However, the CBSS phase remains stable, as the temperature at this point has not yet reached its critical melting value under the given dipole interaction strength.

Then the polarization angle was increased to $\theta = 45^\circ$, and the dipole interactions were set to $V_x = -2V_y = 0.2U$. Figure 7 shows the thermal melting of both the SDW and SSS phases. At low temperature ($k_{\mathrm{B}}T = 0.02U$), as shown in figures 7(a)-(c), an NF phase emerges in the parameter space between the MI and SDW lobes (green shaded regions). The SSS phase is thermally stable and occupies narrow regions adjacent to the SDW lobes. In contrast, for the repulsive TBI with $W = 0.5U$ (figure 7(b)), both the MI and SDW phase shift toward larger hopping values $J/U$, whereas the SSS phase regime contracts. Conversely, for the attractive TBI with $W = -0.1U$ (figure 7(c)), the MI and SDW phases contract toward smaller $J/U$ values, and the SSS phase region expands. This indicates that attractive TBI can stabilize the SSS phase at finite temperatures. At intermediate temperatures ($k_{\mathrm{B}}T = 0.13U$ ), the SDW and MI phases are completely suppressed and replaced by the NF phase, whereas certain SSS lobes remain thermally stable, as shown in the phase diagrams of figures 7(d)-(f). Under repulsive TBI, the NF region expands while the SSS region shrinks. In contrast, with attractive TBI, the NF region shrinks and the SSS region expands. For instance at $\mu/U = 3.5$, the NF region broadens from $J/U = 0.014$ to $0.027$ in figure 7(e), narrows from $J/U = 0.014$ to $0.012$ in figure 7(f). Above the critical temperature$(k_{\mathrm{B}}T_\mathrm{c} \gt 0.16U)$, the crystalline order of the SSS phase is destroyed, completely melting into the NF phase. At $k_{\mathrm{B}}T = 0.3U$, as shown in figures 7(g)-(i), only the NF and SF phases remain, as thermal fluctuations dominate over quantum fluctuations during the phase transition. A comparison with figure 6 shows that at $k_{\mathrm{B}}T = 0.3U$, the CBSS phase in a Bose gas with anisotropic DDI persists up to relatively high temperatures, while the SSS phase is suppressed only at much higher temperatures, resulting in a direct NF-SF transition. In contrast, for an isotropic dipolar Bose gas at the same temperature $k_{\mathrm{B}}T = 0.3U$, thermal fluctuations destroy the CBSS phase [33]. This indicates that the melting of the SS phase depends strongly on both the anisotropy and the strength of the dipolar interaction. Considering the experimental parameters for realizing a dipolar gas in an optical lattice [17, 46], the temperature can be estimated to satisfy $k_{\mathrm{B}}T \lt 0.4U$. For instance, in experiments with $ ^{168}\text{Er}$ atoms, the temperature corresponds to $k_{\mathrm{B}}T \approx 0.3U$, which lies in the nanokelvin range [49]. This ensures that all phases induced by the DDI and TBI can be accessed in ultracold experiments with dipolar bosons.