---
abstract: 'We identify the quantum phases in a binary mixture of dipolar bosons in 2D optical lattices. Our study is motivated by the recent experimental realization of binary dipolar condensate mixture of Er-Dy \[Phys. Rev. Lett. 121, 213601 (2018)\]. We model the system using the extended two-species Bose-Hubbard model and calculate the ground state phase diagrams using mean field theory. For selected cases we also obtain analytical phase boundaries using the site decoupled mean field theory. For comparison we also examine the phase diagram of two-species Bose-Hubbard model. Our results show that the quantum phases with the long range intra-species interaction phase separate with no phase ordering. The introduction of the long range inter-species interaction modifies the quantum phases of the system. It leads to the emergence of phase separated quantum phases with phase ordering. The transition from the phase separated quantum phases without phase ordering to phase ordered ones breaks the inversion symmetry.'
author:
- Rukmani Bai
- Deepak Gaur
- Soumik Bandyopadhyay
- Hrushikesh Sable
- 'K. Suthar'
- 'D. Angom'
bibliography:
- 'ref.bib'
title: 'Segregated quantum phases of dipolar bosonic mixtures in two-dimensional optical lattices'
---

Introduction
============

The Bose-Hubbard model (BHM) [@hubbard_63; @fisher_89] describes the physics of ultracold bosonic atoms trapped in optical lattices [@jaksch_98]. The variation of the hopping term, equivalent to kinetic term in continuum models, in the BHM drives a quantum phase transition from the Mott insulator (MI) to the superfluid phase (SF) phase. And this transition has been experimentally observed [@greiner_02]. The inter-particle interaction in the BHM is onsite or contact in nature. The introduction of the nearest neighbour (NN) interaction in the BHM generates two more phases, density wave (DW) and supersolid (SS) phase. This model with the NN interactions is referred to as the extended BHM (eBHM) [@kuhner_00] and shows rich physics compared to the BHM. Such a model captures the physics of dipolar ultracold quantum gases in optical lattices [@baier_16; @bandyopadhyay_19]. A more complex system, ideal to model several condensed matter systems, is to fill the optical lattice with two species Bose-Einstein Condensate (TBEC). A TBEC could be condensate mixture of two different atomic species [@modugno_02; @lercher_11; @mccarron_11; @pasquiou_13; @wacker_15; @wang_16], two hyperfine states of an atom [@myatt_97; @hall_98; @stamper_98; @stenger_98; @maddaloni_00; @delannoy_01; @sadler_06; @mertes_07; @anderson_09; @tojo_10] or two different isotopes of an atomic species [@papp_08; @handel_11; @sugawa_11a]. It was experimentally first realized in the two hyperfine states $|F = 2, m_F = 2 \rangle$ and $|F = 1, m_F = -1 \rangle$ of $^{87}$Rb atom [@myatt_97]. The TBECs, in the weakly interacting continuum systems, have been used to investigate novel phenomena such as pattern formation [@sasaki_09; @gautam_10_1; @gautam_10_2; @ronen_08; @hoefer_11; @hamner_11; @de_14], phase separation [@ho_96; @ao_98; @gautam_11; @roy_15_2; @bandyopadhyay_17; @papp_08; @tojo_10; @mccarron_11; @wacker_15; @wang_16], nonlinear dynamical excitations [@gautam_12; @gautam_13; @roy_14_2; @kuopanportti_19; @mertes_07; @eto_16a; @eto_16], collective excitations [@roy_14_1; @suthar_15; @roy_15_1; @suthar_16; @suthar_17; @pal_17; @pal_18], Kibble-Zurek mechanism [@nicklas_15], and the production of dipolar molecules [@molony_14; @guo_16; @will_16]. The phase separation, among all the phenomena is a unique property of TBECs. In this work we study the TBECs trapped in the optical lattices, and can be described by the BHM with appropriate modifications. The experimental realization of TBECs in optical lattices are reported in Ref. [@catani_08; @gadway_10] and early theoretical studies are presented in Ref. [@altman_03; @chen_03; @kuklov_03; @kuklov_04]. A remarkable recent achievement related to TBECs is the experimental realization of TBEC with dipolar quantum mixtures of Er-Dy was reported in a recent work [@trautmann_18].

The physics of the two-species BHM (TBHM), the lattice counterpart of TBEC, in 1D has been investigated in detail [@mishra_07; @zhan_14; @wang_14]. And, there has been some works on 2D as well [@altman_03; @chen_03; @kuklov_03; @kuklov_04; @isacsson_05; @iskin_10; @anufriiev_16]. The phase diagram of TBHM shows different combinations of mixed MI-SF phases apart from the MI and SF phases. And, these have been investigated using quantum Monte-Carlo [@kuklov_03; @kuklov_04], mapping to spin systems [@altman_03], and with mean field theory [@chen_03; @isacsson_05; @iskin_10; @pai_12; @anufriiev_16]. These studies, except for Ref. [@pai_12], considered homogeneous systems. However, hitherto the phenomenon of phase separation in 2D TBHM is yet to be investigated in detail.

The quantum phases of TBHM in the phase separated domain, unlike in the TBECs, do not show segregation into two spatial domains. We attribute this to the lack of long range interaction in the TBHM. The simplest modification to include the effect of long range interaction is to add nearest neighbour interaction. The eBHM, as mentioned earlier, supports two more quantum phases density wave (DW) [@capogrosso_10; @flottat_17; @iskin_11] and supersolid (SS) phase [@ng_08; @iskin_11; @yamamoto_11; @boninsegni_12; @leonard_17]. The DW phase is an insulating phase similar to the MI phase but it has crystalline order or diagonal long range order. And, the SS phase is a compressible phase with both diagonal and off diagonal long range order. In a recent study of extended TBHM (eTBHM) [@wilson_16], it was shown that the SS phase exists for small value of NN interaction. In this work, the NN interaction was limited to either one of the species or between the species. We address this research gap by including all the possible intra- and inter-species NN interactions. Such a model is apt to describe the physics of dipolar Bose-Bose mixtures in optical lattices. An example of such a combination is the recently realized Er-Dy mixture [@trautmann_18]. An important result of our work is the possibility to realize compressible and incompressible quantum phases with spatial segregation. Such a phase could be instrumental in examining superfluid instabilities and other non-equilibrium properties in the lattice models of quantum liquids.

The remaining of the paper is organized into four sections. In the Section \[sec\_theory\] we describe the zero temperature Hamiltonian of the TBHM and discuss the Gutzwiller mean-field theory of the model. We then discuss the mean-field decoupling theory to calculate the compressible-incompressible phase boundaries analytically. This is followed by a brief discussion on the characterization of quantum phases. The phase diagrams of TBHM are discussed in the Section \[sec\_phdiag\]. The Section \[sec\_phdiag\_ext\] discussed the phase diagram of eTBHM. In particular, the miscible and immiscible phases. We end the paper with conclusion in Section \[sec\_conclude\].

Theory {#sec_theory}
======

TBHM Hamiltonian
----------------

At zero temperature, the TBHM Hamiltonian which describes the physics of a TBEC in 2D optical lattice is [@damski_03] $$\begin{aligned}
\hat{H}^{{\rm TBH}} &=& -\sum_{p, q, k}\bigg [ \Big( J_x^k 
              \hat{b}_{p+1, q}^{\dagger k}\hat{b}_{p, q}^k + {\rm H.c.}\Big)
              + \Big( J_y^k \hat{b}_{p, q+1}^{\dagger k}
              \hat{b}_{p, q}^k  \nonumber\\ 
              &&+ {\rm H.c.}\Big) - \frac{U_{kk}}{2}
              \hat{n}_{p, q}^k (\hat{n}_{p, q}^k-1) + 
              \tilde{\mu}^k_{p,q}\hat{n}_{p, q}^k\bigg]
              \nonumber\\
              &&+\sum_{p, q} U_{12}\hat{n}_{p, q}^1 \hat{n}_{p, q}^2,
\label{tbhm}  \end{aligned}$$ where $k = 1$,$2$ is the species index, $(p,q)$ are the lattice indices, $J_x^k$ ($J_y^k$) is the nearest neighbour (NN) hopping strength along $x$ ($y$) directions, $\hat{b}^{\dagger k}_{p,q}$ ($\hat{b}^k_{p,q}$) is the creation (annihilation) operator and $\hat{n}_{p, q}^k$ is the number operator at ($p,q$)th site. $U_{kk}$ is intra-species interaction strength, and $U_{12}$ is inter-species interaction strength between two species. Further, $\tilde{\mu}^k_{p,q} =\mu^k - \varepsilon_{p,q}^k$, is the local chemical potential at each site for the two species where $\varepsilon_{p,q}^k$ is envelop potential for the species. For a system of $K \times L$ lattices sites, the index along $x$ ($y$) has values $p=1,\ldots K$ ($q=1,\ldots L$). The unique feature of the TBECs is the phase separation and for continuum systems, the criterion for phase segregation is $U_{12}^2 > U_{11}U_{22}$ [@ho_96; @Trippenbach_00]. Otherwise, the TBEC is in the miscible phase. For the case of strongly interacting TBECs in optical lattices, described by the above Hamiltonian, we show the existence of different phases in both the miscible and immiscible domains.

To obtain the ground state of the TBHM Hamiltonian in Eq. (\[tbhm\]), we use single site Gutzwiller mean field (SGMF) theory [@rokhsar_91; @sheshadri_93; @bai_18; @pal_19; @bandyopadhyay_19; @kuldeep_19]. The starting point of this theory is to separate the operators into mean-field and fluctuation operator components as $\hat{b}_{p, q}^k = \phi_{p,q}^k + \delta \hat{b}_{p, q}^k\;$ and $\;\hat{b}^{\dagger k}_{p,q} = \phi_{p,q}^{k*} + 
\delta\hat{b}_{p,q}^{\dagger k}$. Then, the Hamiltonian in Eq. (\[tbhm\]) is reduced to sum of single-site mean-field Hamiltonian $$\begin{aligned}
\hat{h}_{p,q}^{{\rm TBH}} &= &- \sum_{k}\left[J_x^k 
                  \left(\hat{b}_{p+1, q}^{\dagger k}\phi_{p,q}^k 
                  + \phi^{k*}_{p + 1, q}\hat{b}_{p, q}^k \right ) + {\rm H.c.}
                  \right.  \nonumber\\
                  && +  J_y^k\left(\hat{b}_{p, q+1}^{\dagger k} \phi_{p,q}^k 
                  +  \phi^{k*}_{p, q+1}\hat{b}_{p, q}^k  \right) + {\rm H.c.}
                  \nonumber \\
                  &&-\left . \frac{U_{kk}}{2}\hat{n}_{p, q}^k
                  \left (\hat{n}_{p, q}^k-1\right) 
                  + \tilde{\mu}^k_{p,q}\hat{n}_{p, q}^k\right]
                  \nonumber\\
                  &&+ U_{12}\hat{n}_{p, q}^1 \hat{n}_{p, q}^2,
\label{ham_ss_tbec}\end{aligned}$$ where, $\phi_{p,q}^k$ ($\phi_{p,q}^{k*}$) is the SF order parameter. With this definition of the single-site mean field Hamiltonian, the total Hamiltonian of the system is $$\hat{H}^{{\rm TBH}} = \sum_{p,q}\hat{h}_{p,q}^{{\rm TBH}}.$$ For the details of the derivations see [@bai_18]. To get the ground state we digonalize the Hamiltonian in Eq. (\[ham\_ss\_tbec\]) at each site. And, for this we use the Gutzwiller ansatz, based on which the ground state at the ($p,q$)th site is [@anufriiev_16] $$|\psi\rangle_{p,q} = \sum_{n_1, n_2 }c^{(p,q)}_{n_1, n_2} 
                       |n_1, n_2\rangle_{p,q}.
  \label{gw_2s}$$ Here $|n_1, n_2\rangle$ is a Fock state, which is direct product of the $n_1$ and $n_2$ occupation number states of the 1st and 2nd species, respectively. The occupation number states $n_k\in [0,N_b -1]$ where $N_b$ is the total number of local Fock states used in the computation and $c_{n_1, n_2}^{p,q}$ are complex co-efficients with $\sum_{n_1,n_2} |c^{(p,q)}_{n_1,n_2}|^2$ = 1. From the ground state, we can compute the new SF order parameter of the two species as

$$\begin{aligned}
    \phi_{p, q}^1& =& _{p,q}\langle\psi|\hat{b}_{p, q}^1|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2}\sqrt{n_1} 
              {c^{(p,q)*}_{n_1-1, n_2}}c^{(p,q)}_{n_1,n_2},
              \label{gw_phi2s_1}\\
    \phi_{p, q}^2& =& _{p,q}\langle\psi|\hat{b}_{p, q}^2|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2}\sqrt{n_2} 
              {c^{(p,q)*}_{n_1, n_2-1}}c^{(p,q)}_{n_1,n_2}.
    \label{gw_phi2s_2}              
  \end{aligned}$$

Similarly, corresponding the lattice occupancies are

$$\begin{aligned}
   \rho_{p, q}^1& =& _{p,q}\langle\psi|\hat{n}_{p, q}^1|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2} n_1 |c^{(p,q)}_{n_1,n_2}|^2,
            \label{num2s_1}\\
   \rho_{p, q}^2& =& _{p,q}\langle\psi|\hat{n}_{p, q}^2|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2}n_2|c^{(p,q)}_{n_1,n_2}|^2. 
   \label{num2s_2}              
  \end{aligned}$$

Using the new SF order parameters, the ground state of the next lattice site is computed and this process is repeated till all the lattices sites are covered. One such sweep is identified as an iteration and we then, start the process again for the next iteration. The iterations are carried out till the convergence criterion $|\phi_{p,q}^{n-1} - \phi_{p,q}^n| \lesssim 10^{-12}$ is satisfied at the $n^{{\rm th}}$ iteration. In the present work, to determine the phase diagrams, we consider lattice system of size $10 \times 10$ and choose $N_b = 7$. That is, $K$ and $L$ are both 10. We find that the phase boundaries remain unchanged when the system size is augmented to $20\times 20$. We also use the augmented system size to validate key findings. In addition, we employ periodic boundary conditions to model an infinite-sized system.

Extended TBHM
-------------

The BHM with NN interaction, referred to as the extended BHM (eBHM), exhibits richer phase diagram than BHM and it has the novel feature of harbouring the SS phase. The phase diagram of the single species eBHM consists of DW, SS, MI and SF phases. Similarly, the extended TBHM (eTBHM) also exhibits these phases as well as miscible and segregated phases and the model Hamiltonian of the system is $$\begin{aligned}
\hat{H}^{{\rm ext}} &=& \hat{H}^{{\rm TBH}} 
                       + \sum_{p, q, k}  V_k \hat{n}^k_{p,q} \biggl [ 
                         \Big ( \hat{n}^k_{p-1,q} + \hat{n}^k_{p+1,q} 
                       + \hat{n}^k_{p,q -1}  \nonumber\\
                    && + \hat{n}^k_{p,q+1}\Big ) +  V_{12} \hat{n}_{p,q}^k
                         \Big(\hat{n}^{3-k}_{p-1,q} + \hat{n}^{3-k}_{p+1,q} 
                       + \hat{n}^{3-k}_{p,q -1} 
                                     \nonumber \\
                    && + \hat{n}^{3-k}_{p,q+1}\Big ) \biggr],
\label{bhm_ext}\end{aligned}$$ here $V_k$ is the intra-species NN interaction strength for both the species, and $V_{12}$ is the inter-species NN interaction strength. In the experiments the relative strengths of the NN interaction is controllable through the Feshbach resonances. Thus, to relate with the experimental observations and predict possible phases we vary the inter- and intra-species interaction strengths. We use SGMF theory to obtain the ground state of the system, then, in this method the total Hamiltonian is the sum of the single-site mean field Hamiltonian $$\begin{aligned}
\hat{h}_{p,q}^{{\rm ext}} &=& \hat{h}_{p,q}^{{\rm TBH}} 
                       + \sum_{k}  V_k \hat{n}^k_{p,q} \biggl [ 
                         \Big ( \hat{n}^k_{p-1,q} + \hat{n}^k_{p+1,q} 
                       + \hat{n}^k_{p,q -1}  \nonumber\\
                    && + \hat{n}^k_{p,q+1}\Big ) +  V_{12} \hat{n}_{p,q}^k
                         \Big(\hat{n}^{3-k}_{p-1,q} + \hat{n}^{3-k}_{p+1,q} 
                       + \hat{n}^{3-k}_{p,q -1} 
                                     \nonumber \\
                    && + \hat{n}^{3-k}_{p,q+1}\Big ) \biggr].
\label{bhm_ext_ss}\end{aligned}$$ For eTBHM also, we diagonalize the Hamiltonian at each site separately, and obtain the ground state. The NN interaction term contribute to the diagonal matrix element. From the single-site wavefunction the SF order parameter and lattice occupancy can be calculated from the expressions in Eqns. (\[gw\_phi2s\_1\]-\[gw\_phi2s\_2\]) and (\[num2s\_1\]-\[num2s\_2\]).

Mean-field decoupling theory
----------------------------

To calculate the phase boundaries between MI and SF phases analytically we use the site decoupled mean-field theory [@oosten_01; @iskin_09; @bandyopadhyay_19]. For this, we adapt perturbative analysis of the mean-field Hamiltonian in Eq. (\[ham\_ss\_tbec\]). It is important to note that the SF order parameter $\phi^{k}_{p,q} $ is zero in the MI phase, but nonzero in the SF phase. So, the vanishing of the SF order parameter $\phi^{k}_{p,q} \rightarrow 0^{+}$ marks the MI-SF phase boundary in the phase diagram. With this consideration, in the site decoupled mean-field theory, the interaction and the chemical potential terms constitute the unperturbed Hamiltonian $\hat{h}_{p,q,0}^{\rm TBH}$. From Eq. (\[ham\_ss\_tbec\]), $$\begin{aligned}
\hat{h}_{p,q,0}^{{\rm TBH}} 
            &= & \sum_{k}\left[ \frac{U_{kk}}{2}\hat{n}_{p, q}^k
                  \left (\hat{n}_{p, q}^k-1\right) 
                  - \tilde{\mu}^k_{p,q}\hat{n}_{p, q}^k\right]
                  \nonumber\\
                  &&+ U_{12}\hat{n}_{p, q}^1 \hat{n}_{p, q}^2,
\label{ham_ss_h0}\end{aligned}$$ which is diagonal with respect to the Fock basis states.Then, the hopping terms in Eq. (\[ham\_ss\_tbec\]) as the perturbation, $$\begin{aligned}
\hat{h}_{p,q,1}^{{\rm TBH}} &= &- \sum_{k}\Bigl [J_x^k 
                  \left(\hat{b}_{p+1, q}^{\dagger k}\phi_{p,q}^k 
                  + \phi^{k*}_{p + 1, q}\hat{b}_{p, q}^k \right ) + {\rm H.c.}
                   \nonumber\\
                  && +  J_y^k\left(\hat{b}_{p, q+1}^{\dagger k} \phi_{p,q}^k 
                  +  \phi^{k*}_{p, q+1}\hat{b}_{p, q}^k  \right) + {\rm H.c.}
                   \Bigr ],
\label{ham_ss_h1}\end{aligned}$$ with the SF order parameter $\phi^{k}_{p,q}$ as the perturbation parameter. Then, from the first-order perturbative correction to the ground state wavefunction (details given in Appendix \[appendix\_a\]), we have $$\!\!\!\!\!\phi^{k}_{p,q} = J\bar{\phi}^{k}_{p,q} 
    \left(\frac{n^{k}_{p,q} +1}{n^{k}_{p,q}U -\bar{\mu}^{k}_{p,q}}
        -\frac{n^k_{p,q}}{(n^k_{p,q} -1)U -\bar{\mu}^{k}_{p,q}}\right),
 \label{order_par_ana}$$ with $$\begin{aligned}
   \bar{\mu}^{k}_{p,q} &= &\tilde{\mu}^{k}_{p,q}-U_{12}n^{3-k}_{p,q},
                          \nonumber \\
   \bar{\phi}^{k}_{p,q} &= &\left(\phi^{k}_{p+1,q}+\phi^{k}_{p-1,q}
                           +\phi^{k}_{p,q+1}+\phi^{k}_{p,q-1}\right).
                           \nonumber\end{aligned}$$ For a homogeneous lattice system $\varepsilon^{k}_{p,q} = 0$. Then, in the MI phase the total density $\rho = \rho^{1}+\rho^{2}$ is integer commensurate and $\phi^{k}_{p,q}=0$. In the SF phase, the order parameter is non-zero and uniform, say $\phi^{k}_{p,q}=\varphi^{k}_{0}$. With these considerations, $\bar{\phi}^{k}_{p,q} = \bar{\phi}^{k} = 4\varphi^{k}_{0}$. Starting from the SF phase, at the SF-MI phase boundary $\varphi^{k}_{0}\rightarrow 0^{+}$. Considering this limit in Eq. (\[order\_par\_ana\]), we obtain the equation which defines the MI-SF phase boundary in terms of $J$ for a particular value of $\mu$.

### TBHM

For the MI phase with $\rho = 2$,in the miscible domain, atoms of both the species fill all the lattice sites. That is, $n^1_{p,q} =n^2_{p,q} =1$. The MI-SF phase boundary is, then, defined by $$\begin{aligned}
   \frac{1}{4J}= \frac{2}{U -\mu +U_{12}} +\,\frac{1}{\mu -U_{12}}.
   \label{tbhm_even}\end{aligned}$$ On the other hand for finite $U_{12}$, the system is in the immiscible domain for the MI phase with $\rho=1$. The density pattern has one atom at each lattice site chosen randomly from the two species. Thus, at a given lattice site $(p,q)$ we can have the occupancies as $n^{1}_{p,q} =1, n^{2}_{p,q} =0$ or $n^{1}_{p,q} =0, n^{2}_{p,q} =1$. In the perturbative analysis, without loss of generality, we consider neighbouring lattice sites are occupied by atoms of different species. This is also one realization of the energetically favourable configuration for $U_{12}<U$. Then, with the correction arising from $b^{\dagger 1}\phi^{1}$, the equation $$\begin{aligned}
   \frac{1}{4J}= \frac{2}{U_{12} -\mu} +\,\frac{1}{\mu},
   \label{tbhm_odd}\end{aligned}$$ defines the phase boundary of the MI lobe with $\rho=1$. Based on similar analysis, we can obtain the phase boundary of other MI lobes. For which we have to use Eq.(\[tbhm\_even\]) and (\[tbhm\_odd\]) for the even and odd integer values of $\rho$, respectively.

### Extended TBHM

The MI-SF boundary of the eTBHM can also be calculated similar to the TBHM. The expression of the order parameter is similar to Eq. (\[order\_par\_ana\]) but $\bar{\mu}^{k}$ is given by $$\begin{aligned}
   \bar{\mu}^{k}_{p,q} &= &\tilde{\mu}^{k}_{p,q}-U_{12}n^{3-k}_{p,q}
        -4V_k n^{k}_{p,q} -4V_{12} n^{3-k}_{p,q},\end{aligned}$$ For MI phase with $\rho=2n$, the occupancies are $n^{k}_{p,q}=n$. Further, assuming $V_1 =V_2$, the MI-SF boundary is given by $$\begin{aligned}
   \frac{1}{4J}= \frac{2}{U -\bar{\mu}} +\,\frac{1}{\bar{\mu}}.
 \label{exTBHM_MI2}\end{aligned}$$ with $\bar{\mu} = \mu -U_{12} -4V_1 -4V_{12}$. A similar analysis of MI phase with odd integer occupancies require specific approximations based on the density configuration. And, a general unified analysis in not applicable.

Characterization of the phases
------------------------------

  --------------------- -------- ---------- ---------------- ----------------
                        $\rho$   $\phi$     $\Delta\rho^k$   $\Delta\phi^k$
  Mott Insulator (MI)            $ = 0$     $ = 0$           $ = 0$
  Superfluid (SF)                $ \ne 0$   $ = 0$           $ = 0$
  Density Wave (DW)              $ = 0$     $ \ne 0$         $ = 0$
  Super Solid (SS)               $ \ne 0$   $ \ne 0$         $ \ne 0$
  --------------------- -------- ---------- ---------------- ----------------

  : Classification of different quantum phases with order parameters for our systems.[]{data-label="table"}

To identify different quantum phases of the system we compute the density contrast $\Delta\rho^k$, order parameter contrast $\Delta\phi^k$ and compressibility $\kappa^k$. To define $\Delta\rho^k$, divide the lattice site occupancies as $$n^k_{p,q}=\begin{cases}
                n^{k,A} & \text{if $(p,q) \in$ sublattice A}, \\
                n^{k,B} & \text{if $(p,q) \in$ sublattice B},
             \end{cases}
   \label{occ_def}$$ then, the density contrast of the $k$th species is $$\Delta\rho^k = n^{k,A}-n^{k,B}.
   \label{den_con}$$ The order parameter contrast is defined similarly as $$\Delta\phi^k = \phi^{k,A}-\phi^{k,B},
   \label{den_con}$$ where, $\phi^{k,A}$ and $\phi^{k,B}$, like in the case of density are the values of the order parameters at lattice sites with $(p,q)$ belonging to sublattices A and B, respectively. The compressibility of each species are calculated using the definition $\partial \mu^k/\partial \rho^k$.

The TBHM, like the single species BHM, shows two phases MI and SF phases. The MI phase is an incompressible phase with integer commensurate density $n^{k,A}=n^{k,B}\in\mathbb{N}$. And, incompressibility implies zero SF order parameter $\phi^{k,A}=\phi^{k,B}=0$. The SF phase, on the other hand is compressible. Hence, it has $n^{k,A}=n^{k,B}\in \mathbb{R}$, $\phi^{k,A}=\phi^{k,B}\in \mathbb{R}$ and $\kappa^k\neq0$. For these two phases, the density and SF order parameters are uniform, so the contrast order parameters $\Delta\rho^k$ and $\Delta\phi^k$ are zero. In the eTBHM, the NN interaction leads to the emergence of two more quantum phases, DW and SS. These two phases have non-uniform density and SF order parameters. As a result the distinguishing features of these phases are non-zero contrast order parameters. The DW phase has integer $n^k_{p,q}$ with $n^{k,A}\neq n^{k,B}$ and $\Delta \rho^k\in\mathbb{N}$. This phase has zero SF order parameter $\phi^{k,A}=\phi^{k,B}=0$ and hence, incompressible. The SS phase has real $n^k_{p,q}$ with $n^{k,A}\neq n^{k,B}$ and $\Delta \rho^k\in\mathbb{R}$. The SF order parameter in this phase is non-zero and non-uniform. This implies that $n^{k,A}\neq n^{k,B}$ and $\phi^{k,A}\neq\phi^{k,B}$. So, both the contrast order parameters are non-zero in the SS phase. For easy reference the properties of the different quantum phases are listed in Table. \[table\].

Phase diagram of TBHM {#sec_phdiag}
=====================

![Phase diagram of TBEC by varying the inter-species interaction strength $U_{12}$. Blue solid lines represent numerically obtained phase boundaries for the mean field hamiltonian. Filled dots marks analytically obtained phase boundaries between compressible and incompressible phases, obtained analytically by perturbative analysis of the mean-field hamiltonian. The odd occupancy MI lobes appear for non zero $U_{12}$ and enlarges with increasing $U_{12}$.[]{data-label="tbec_ph"}](tbc_ph.pdf){width="8.0cm"}

To compute the ground state wavefunction and determine the phase, we initialize the SF order parameter $\phi$. This, then, defines the Hamiltonian in Eq. (\[ham\_ss\_tbec\]) and Hamiltonian matrix elements are computed using Gutzwiller wavefunction in Eq. (\[gw\_2s\]). By digonalizing the Hamiltonian matrix for each site we obtain the ground state wavefunction. From the results, the MI-SF phase boundary is identified based on the SF order parameter and the lattice occupancy. For the incompressible MI phase, at each lattice site, $\phi$ is zero and $\rho$ is integer commensurate. For the SF phase, $\phi$ is nonzero and $\rho$ is real commensurate. The phase diagrams of BH Hamiltonian of TBEC in Eq. (\[ham\_ss\_tbec\]) for different values of $U_{12}$ are shown in Fig. \[tbec\_ph\].

For simplicity, we consider symmetric hopping $J_x^k = J_y^k = J$, equal chemical potential $\tilde{\mu}^1_{p,q} = \tilde{\mu}^2_{p,q} = \mu$ and identical intra-species interactions $U_{kk} = U$. We scale all the energies with $U$, and define the phase diagram in the $J/U-\mu/U$ plane. The phase diagram consists of a sequence of MI lobes having integer $\rho$. Without the inter-species interaction $U_{12} = 0$, as shown in Fig. \[tbec\_ph\](a), the phase diagram is equivalent to the case of single species. But, with twice the occupancy. That is the MI lobes have $\rho=2n$ with $\rho^1=\rho^2=n$ and $n\in\mathbb{N}$. So, the lowest MI lobe has $\rho=2$ and each lattice has one atom from each of the two species. As a result, the phase diagram is identical to single species. With the introduction of the inter-species interaction ($U_{12} \neq 0$) the half filled lobes like $\rho = 1$ emerge in the phase diagram with $0 < \rho^1 <1$, and then, $\rho^2 = 1-\rho^1$. This is discernible for $U_{12} =0.4 U$ from the Fig. \[tbec\_ph\](b). Based on the form of the interactions in the Hamiltonian of the system, the energies of system is degenerate for all the possible combinations of $\rho^1$ and $\rho^2$ in the allowed ranges. For example, with $U_{12}= 0.4U$ and for $\mu/U=0.2$, $J/U=0.01$ we observe $0.33\lesssim \rho^1\lesssim 0.7$. In the figure, the half filling lobe $\rho^k = 0.5$ and $\rho = 1$ at $J/U=0$ lies in the domain $0\leqslant\mu/U \leqslant 0.4$. In general, in the miscible domain, the half filling lobe $\rho =1$ at $J/U=0$ lies in the domain $0\leqslant\mu/U \leqslant U_{12}/U$. The other MI lobes with higher $\rho$ occur at the higher values $\mu/U$. In general, the MI lobes have $\rho=n$ with $n\in\mathbb{N}$ and $\rho^k = n/2$. Thus, for MI lobes with odd $n$ the average occupancy of each species is half integer.

With increasing $U_{12}$ the MI lobes with odd integer occupancies grows in size, but the size of MI lobes with even integer occupancies remain the same till $U_{12}=U_{kk}$ but shifts to higher $\mu/U$. This can be understood from Eq. (\[order\_par\_ana\]). The trend is discernible from the phase diagrams in Fig.\[tbec\_ph\](b-c). This, in the case of weakly interacting TBECs, is equivalent to a march towards phase separation [@ho_96; @ao_98; @gautam_11; @roy_15_2; @bandyopadhyay_17]. For $U_{12}>U_{kk}$, the criterion for phase separation, the size of the MI lobe $\rho =2$ is different. But, once the phase separation criterion is met, there is no change in the phase diagram with further increase in $U_{12}$. As an example the phase diagram for $U_{12}=1.2U$ is shown in Fig.\[tbec\_ph\](d). The lobes in this phase diagram are the same as in Fig.\[tbec\_ph\](a). The only difference is the occupancy is $\rho=n$ with $n\in\mathbb{N}$ and $\rho^k=n/2$. As a result, the density pattern of the lowest MI lobe ($\rho=1$) has one atom at each lattice site chosen randomly from the two species. The important point is, the MI lobes have the same sizes for $U_{12}=0$ and $U_{12}\geqslant U_{kk}$. But, the occupancy and hence the density patterns are different.

Phase diagram with long range interactions {#sec_phdiag_ext}
==========================================

$V_{12}=0$
----------

![Phase diagram of extended BHM for TBEC at the different inter-species interaction strength $U_{12}$ and for inter-species NN interaction $V_{12} = V_{21} = 0$, $V_1 = V_2 = 0.05U$. Brown lines forms the boundary of region comprising of incompressible phases (MI, DW). Around DW phase, SS phase exists and the boundary between SS and SF phases is represented by green lines. The SS region around the DW region enlarges with increasing $U_{12}$. In DW(n,0) phase both species have DW(n,0) pattern and in MI(1,1) phase both species have uniform unit occupancy. For $U_{12} > U$ the MI lobes are replaced by DW as seen in Fig(d).[]{data-label="ebhm_v12_0"}](ebh_v12p0_ph.pdf){width="8.0cm"}

The ground state of the extended BHM Hamiltonian in Eq. (\[bhm\_ext\_ss\]), like in the previous case, is obtained using the Gutzwiller ansatz. The long range interactions in the extended BHM introduce two more phases, DW and SS, in the phase diagram. To analyse and highlight the effect of long range intra- and inter-species interactions, we first consider the case of $V_{12} = 0$. And, set the intra-species NN interaction strength $V_k = 0.05U$. Then, vary the inter species onsite interaction strength $U_{12}$, which can be achieved in experiments through Feshbach resonance. The choice of low value of $V_k$ is based on the parameters realized in dipolar BEC experiments [@baier_16]. In these experiments, $V/\hbar$ is in the range $\sim 10-100$ Hz, whereas $U/\hbar$ has typical values in kHz. In addition, this choice of parameters has the unique possibility to study the MI-DW quantum phase transition by changing $U_{12}$ keeping $V_k$ fixed. This is to be contrasted with the eBHM, where the NN interaction strength $V\ge0.25U$ [@iskin_11; @kuldeep_19] marks the critical point for such quantum phase transitions. Like in the BHM case, we consider symmetric hopping $J_x^k = J_y^k = J$, identical chemical potential $\tilde{\mu}^k_{p,q} = \mu$ and $U_{kk} = U$. The phase diagram for $U_{12} =0$ is shown in Fig. \[ebhm\_v12\_0\] (a). It is identical to the phase diagram of single species extended BHM [@iskin_11; @kuldeep_19] and consists of the DW(1,0), MI(1,1), DW(2,1), MI(2,2), SS (green line) and SF phases. In the figure, the SS phase occurs as a thin envelope around the DW lobes. On increasing $V_k$ but keeping the other parameters fixed, the size of the DW lobes and the accompanying envelope of SS phase are enhanced. However, the MI lobes shrinks and disappear from the phase diagram. This is due to the higher energy cost of having commensurate occupancy due to the intra-species NN interaction. The same effect is reported in the single species eBHM [@iskin_11; @kuldeep_19].

The DW phases with $U_{12} =0$ are four-fold degenerate. Two of the states have $\Delta \rho^1=\Delta\rho^2$ and the other two have $\Delta \rho^1=-\Delta\rho^2$. For both set of states, one of the degenerate states is obtained by shifting both of the species by one lattice constant either along $x$ or $y$ direction. For the $\Delta \rho^1=\Delta\rho^2$ states, the occupancies of the two species at each lattice sites are the same $n^1_{p,q}=n^2_{p,q}$. From this the $\Delta \rho^1=-\Delta\rho^2$ states are obtained after translation of one of the species by one lattice constant either along $x$ or $y$ direction. Thus, in the latter we have $n^{1,A}=n^{2,B}$ and $n^{1,B}=n^{2,A}$. It is to be noted that the $\rho=1$ phase of the TBHM has the same average density as the DW(1,0). However, the two have different symmetries. The $\rho=1$ phase of the TBHM has atoms from the two species with random occupancies and has no diagonal long range order. But, the DW(1,0) has diagonal order arising from the non-zero $\Delta\rho^k$. As an example consider the DW(1,0) phase, the two degenerate states correspond to $\Delta\rho^1=\Delta\rho^2=1$ and $\Delta\rho^1=-\Delta\rho^2=1$. At higher $\mu$, the DW(2,1) intervenes the transition from MI(1,1) to MI(2,2) phase.

To study the effect of the inter-species interaction we increase $U_{12}$ retaining $V_{12}$ and $V_k$ fixed at $0$ and $0.05U$, respectively. The phase diagram corresponding to $U_{12} = 0.4U$ is shown in Fig. \[ebhm\_v12\_0\] (b). At finite $U_{12}$ the MI phase is energetically costly due to repulsion between atoms of the two-species co-existing on the same lattice site. So MI phase shifts to higher $\mu/U$ values with increasing $U_{12}$ which can be understood from Eq. (\[exTBHM\_MI2\]). As seen from the figure, the finite $U_{12}$ enhances the DW(1,0) lobe. The finite $U_{12}$ also lifts the degeneracy of the DW states, and the state with $n^1_{p,q}=n^2_{p,q}$ has higher energy. So, the density of the DW states with finite $U_{12}$ has $n^{1,A}=n^{2,B}$ and $n^{1,B}=n^{2,A}$.

The MI(1,1) lobe remains unchanged in size, but, it is shifted upward in the phase diagram. The shift is attributed to the increase in effective chemical potential arising from the interaction energy associated with finite $U_{12}$. Similar trend, enhancement of DW(1,0) lobe, occur in the case of $U_{12}=0$ on increasing $V_k$. In addition to the MI phase, the DW(2,1) and similar DW phase with non-zero $n^{k,A}$ and $n^{k,B}$ are also energetically disfavoured. However, the most important feature is the emergence of prominent SS phase envelope around each of the DW lobes. On increasing $U_{12}$ further, as seen from the Fig. \[ebhm\_v12\_0\] (c-d), the MI lobes are transformed into DW lobes. And, at higher $U_{12}$ only the DW(n,0) phase, with $n\in\mathbb{N}$, are present in the system. The domain of the SS phase also increases. Ultimately, the SS envelopes around each of DW lobes merge into single large SS domain and this is discernible in these figures.

$V_{12}>0$
----------

One of the phenomena unique to TBEC is the phase separation. This provides important insights to understand novel phenomena in nonlinear dynamics, pattern formation, quantum phase transitions in condensed matter systems, etc. [@gautam_10_1; @gautam_10_2; @gautam_11; @roy_15_2; @bandyopadhyay_17; @gautam_12; @gautam_13; @roy_14_2; @kuopanportti_19; @roy_14_1; @suthar_15; @roy_15_1; @suthar_16; @suthar_17; @pal_17; @pal_18] Phase separation of TBECs in the weakly interacting regime, as mentioned earlier, is well studied. This, however, is not the case for the strongly interacting two-species ultracold atoms in optical lattices.

![Phase diagram of extended BHM for TBEC at the different inter-species interaction strength $U_{12}$ and for inter-species NN interaction $V_{12} = V_{21} = 0.05$, $V_1 = V_2 = 0.05U$. The incompressible (MI, cDW) and compressible phase(SS, SF) regions are separated by brown lines. In cDW phase the two species occupy lattice sites randomly in such a way such that total density $\rho= \rho^1 +\rho^2$ have DW pattern. Around DW phase, SS phase exists and the boundary of the SS region is marked by green lines. (d) For $U_{12} = 1.2$ the DW and SF phases are phase separated.[]{data-label="ebhm_v12_0p05"}](ebh_v12p0p05_ph.pdf){width="8.0cm"}

As discussed earlier, in the TBHM we observe phase separation in the SF phase, where the density of the two species are spatially separated into two domains. The phase separated MI phases, on the other hand, has random filling of the two species and are not separated into two distinct domains. The inclusion of the NN interactions modify the density distribution of the MI phases in the phase separated domain. To study this, we solve the Eq. (\[bhm\_ext\_ss\]) with finite $V_{12}$ and keep it fixed to a value of $0.05U$. We, then, increase the inter-species interaction $U_{12}$ from the miscible domain $U_{12}^2 < U_{11}U_{22}$ to the immiscible domain $U_{12}^2 > U_{11}U_{22}$. The phase diagrams for selected values of $U_{12}$ are shown in the Fig. \[ebhm\_v12\_0p05\].

### Miscible phase

In the miscible domain, $U_{12}^2 < U_{11}U_{22}$, the phase diagram has lobes of incompressible quantum phases having $\rho=n$ with $n\in \mathbb{N}$. These lobes are similar to those in the TBHM. In the present case, however, the $\rho=n$ lobes are intervened by lobes of DW quantum phases with half-integer total average occupancies $\rho = (2m+1)/2$ with $m\in \{0,\mathbb{N}\}$. The total occupancy $n_{p,q} = n^1_{p,q} + n^2_{p,q}$ of these phases have diagonal long-range order. This is essentially induced by the non-zero inter-species NN interaction, $V_{12}>0$. The particle densities $n^k_{p,q}$, however, posses no diagonal long-range order. For this reason we refer to these as correlated DW (cDW) phases. This is to distinguish from the DW phases with $V_{12}=0$, in which case $n^k_{p,q}$ have diagonal long range order. Another important property of the cDW lobes is, these are surrounded by a thin envelope of SS phase. As an example, the phase diagram for $U_{12}=0.9U$ is shown in Fig. \[ebhm\_v12\_0p05\](a). In the figure, the cDW(1,0) has the lowest average occupancy $\rho=1/2$. One of the possible density distribution of this phase is $n^{k,A}=0$. And, at the other sub-lattice the occupancy is $n^{2,B}_{p,q} = 1 - n^{1,B}_{p,q}$. The values of $n^{1,B}_{p,q}$ is either 0 or 1 distributed randomly. And, the random distribution implies that there is no diagonal long range order. In other words, the lattice occupancies of the individual species are not structured but the total lattice occupancy is a structured quantum phase. Around the cDW phase, as $J/U$ is increased for fixed $\mu/U$, the quantum fluctuations drive a second order quantum phase transition from cDW to the SS phase. For the SS phase around the cDW(1,0) phase, the occupancies of the two sub-lattices are identical, and lie in the range $0\lesssim n^{1,A}_{p,q}=n^{2,A}_{p,q}\lesssim 0.25$ and $0.25\lesssim n^{1,B}_{p,q}=n^{2,B}_{p,q}\lesssim 0.50$. Hence, both the species have the same diagonal long-range orders. Here, the occupancies of the SS phase are defined over a finite range due to it’s finite compressibility. The SF order parameters, although different in value, follow similar trends $\phi^{1,A}_{p,q}=\phi^{2,A}_{p,q}$, $\phi^{1,B}_{p,q}=\phi^{2,B}_{p,q}$ and $\phi^{k,B}_{p,q}\neq\phi^{k,A}_{p,q}$. In short, the fluctuations drive the cDW(1,0) phase with random integer $n^{k,B}_{p,q}$ to identical occupancies. And, $n^{k,A}_{p,q}$ also acquire non-zero value. On increasing $J/U$ further, the quantum fluctuations drive another phase transition from SS to SF phase. In this transition, the diagonal long-range order is destroyed and translational invariance of the system is restored.

![Phase separation with side by side pattern of species occupancies, obtained with periodic boundary conditions along both $x$ and $y$ axes. The density distribution of the species over lattice sites is shown in Fig (a,b,c) for DW(2,1) phase, in Fig(d,e,f) for SS phase and in Fig(g,h,i) for SF phase.[]{data-label="dw_ss_sf_occ"}](dw_ss_sf_occupancy.pdf){width="8.0cm"}

The insulating phase with average occupancy $\rho=1$, has uniform total lattice occupancy $n_{p,q}=n^1_{p,q} + n^2_{p,q}=1$. And, the occupancies of the two species satisfy the condition $n^1_{p,q}=1-n^2_{p,q}$ with $n^2_{p,q}\in \{ 0,1\}$, where the values between the two possibilities are chosen at random. Thus, this phase like the conventional MI phase with integer commensurate integer occupancies, but in terms of the total occupancy $n_{p,q}$. Similar to cDW phase, we refer to this phase as correlated MI (cMI) phase. This implies that increasing the chemical potential or on adding more particles to the system, at a fixed but low $J/U$, the system starting from cDW(1,0) passes through SS, SF and then to the $\rho=1$ phase. At still higher $\mu$, the cDW(2,1) phase appears. The total occupancies of the two sub-lattices in this quantum phase are $n^A_{p,q} = n^{1,A}_{p,q} + n^{2,A}_{p,q} =2$ and $n^B_{p,q} = n^{1,B}_{p,q} + n^{2,B}_{p,q} =1$. This implies that both the species have same occupancies in the $A$ sub-lattice $n^{1,A}_{p,q} = n^{2,A}_{p,q} = 1$. And, it is equivalent to the DW(2,0) phase in the eTBHM with $V_{12}=0$. From this phase we obtain the cDW(2,1) phase by randomly adding one atom of either species at the $B$ sub-lattice sites. That is $n^{1,B}_{p,q}=1-n^{2,B}_{p,q}$ with $n^{2,B}_{p,q}\in \{ 0,1\}$, where the values between the two possibilities are chosen at random. So, effectively, the cDW(2,1) is a superposition of the DW(2,0) with the cDW(1,0). At higher $\mu$ the other lobes with increasing $\rho$ appear. And, these have similar occupancies and order parameter structure as the lobes with lower $\rho$. It is to be highlighted that the phase diagrams are different, qualitatively and quantitatively, from the two-species BHM where only one of the species is dipolar [@wilson_16].

![Phase separation with side by side pattern of species SF order parameter, obtained with periodic boundary conditions along both $x$ and $y$ axes. The SF order parameter at the lattice sites is shown in Fig (a,b,c) for SS phase and in Fig(d,e,f) for SF phase.[]{data-label="ss_sf_phi"}](ss_sf_phi.pdf){width="8.0cm"}

### Immiscible phase

The criterion for phase separation of the two species in the TBECs or weakly interacting domain is $U_{12}^2>U_{11}U_{22}$ [@ho_96; @ao_98]. And, as discussed earlier, at phase separation the atoms of different species do not occupy the same lattice site. This is the energetically favourable configuration. However, the local nature of the inter-particle interaction preserves the inversion symmetry and the species do not separate into two spatial domains. In the TBECs or weakly interacting domain, the contact interaction is sufficient to break the inversion symmetry and leads to formation of two spatial domains [@papp_08; @tojo_10; @mccarron_11; @wacker_15; @wang_16] at phase separation. The introduction of the long range inter-species interaction ($V_{12}>0$) in the eTBHM introduces the possibility to lower the energy of the density configurations which breaks inversion symmetry. Thus, there is phase ordering of the two species.

In the present case, for the parameters considered ( $U_{11}=U_{22}=U$) the phase separation criterion is equivalent to $U_{12}>U$. This choice of the parameters, as a representative case, capture the key qualitative and quantitative features of the eTBHM. More importantly, the long range nature of the $V_{12}$ introduce phase ordering at phase separation. As an example, the phase diagram for $U_{12}=1.2U$ is shown in Fig. \[ebhm\_v12\_0p05\](b). The structure of the insulating or incompressible and compressible phases are similar to the case of $U_{12}<U$ as shown in Fig. \[ebhm\_v12\_0p05\](a). But, there is one key difference, the cDW, SS and SF phases in Fig. \[ebhm\_v12\_0p05\](b) are phase separated. This is the combined effect of the onsite and long range inter-species interactions. And, this is indicated in the phase diagram with annotation PS (phase separated). But, the insulating phases with $\rho=1$ and $\rho=2$ are not phase separated. In $\rho=1$ phase, like in the case of $U_{12}<U$, each lattice site is singly occupied by an atom from the two species chosen randomly. If the phase separation is along one of the axes, say $x$-axis, the DW($n_A, n_B$) phase has occupancies $$n^k_{p,q}=
      \begin{cases}
          \Theta\left [(-1)^k (p-(K-1)/2)\right ]n_A 
                &  \text{for $(p,q) \in$ A}, \\
         \Theta\left [(-1)^k (p-(K-1)/2)\right ]n_B 
                & \text{for $(p,q) \in$  B},
      \end{cases}
   \label{dw_ps_stripe1}$$ here, $k$, as defined earlier, is the species index, $K$ is the size of the system along $x$-axis, and $n_A$ and $n_B$ are integers with $n_A\neq n_B$. The ground state is doubly degenerate as the above density configuration has the same energy when the species are interchanged. The occupancies of the phase separated SS and SF phases can also be defined in a similar way. However, in these two phases $n_A$ and $n_B$ are real. Furthermore, in the SS phase $n_A\neq n_B$ but in the SF phase $n_A= n_B$. The SF order parameters in the SS and SF phases are also defined in the same form. The presence of the Heaviside step functions in Eq. (\[dw\_ps\_stripe1\]) indicates inversion symmetry is broken. The Hamiltonian is, however, invariant under the inversion symmetry. Thus, the phase mixed to separation transition breaks the inversion symmetry spontaneously. And, the observed ground state is one of the degenerate configurations.

As an example, the phase diagram in the immiscible domain $U_{12} = 1.2 U$ is shown in Fig. \[ebhm\_v12\_0p05\] (b). In the phase diagram, the global features of the phase domains are qualitatively similar to phase diagram in the miscible parameter domain $U_{12} = 0.9 U$ shown in Fig. \[ebhm\_v12\_0p05\] (a). There is, however, an important difference. All the phases in the figure are phase separated and this is indicated in the phase diagram with PS. In the SF phase, phase separation occurs for all parameter domain. The density profiles of the DW(2,1) phase, and the SS and SF phases around it are shown in Fig.(\[dw\_ss\_sf\_occ\]). In the figure, consider the lattice sites with odd (even) values of $(p+q)$ as the $A$ ($B$) sublattice. And, for better representation of the density orders of the structured phases, we consider a system size of $10\times 10$. Then, from the density pattern in Fig.(\[dw\_ss\_sf\_occ\])(a) and (b), both the species have occupancies $n_A=2$ and $n_B=1$. And, as it is phase separated, from Eq. (\[dw\_ps\_stripe1\]), $$\!\!n^k_{p,q}=
     \begin{cases}
         \Theta\left [(-1)^k(p-\frac{9}{2}) \right ]2 &  
                 \text{for odd $(p+q)$}, \\
         \Theta\left [(-1)^k(p-\frac{9}{2}) \right ]1 &  
                 \text{for even $(p+q)$}.
      \end{cases}
   \label{dw21_sp12}$$ The density pattern shown in the figures Fig.(\[dw\_ss\_sf\_occ\])(a-c) correspond to the parameters $\mu/U=1.35 $ and $J/U=0.010$. The above occupancies of the species implies that each of the species are confined within a subsystem of $5\times 10$ lattice. The other species, as we apply periodic boundary conditions along both the directions, effectively provides a confining potential. This is better visualised when the system is mapped to a torus. Then, phase separation along one of the axes, divides the torus into two equal halves. And, each one occupied by one of the species. For such a configuration, there are two inter-species boundaries which segregate the two species. Thus, with a $10\times 10$ system size, the total length of the boundary is $20a$, where $a$ as defined earlier is the lattice constant. From the figure it is evident that other configuration is phase separation along the diagonal. This, however, is energetically not favourable as it has larger interface energy due to longer boundary $10(2+\sqrt{2})a$.

For the same value of chemical potential $\mu/U=1.35$, on increasing the hopping amplitude to $J/U=0.011$ we are in the SS phase domain. It is also phase separated and the lattice site occupancies has similar form as in Eq.(\[dw21\_sp12\]). The occupancies are real, have checker board order and are shown in Fig.(\[dw\_ss\_sf\_occ\])(d-f) . Another important point is, as seen from the figures, boundary effects are present in the SF order parameter. The reason is that the effective potential which segregates the two species is like a soft boundary condition. And, this is due to the long range interspecies interaction. The SS phase is a superfluid phase with diagonal long range order, and hence, has non-zero SF order parameter $\phi^k_{p,q}$. The SF order parameters of the two species are shown in Fig.\[ss\_sf\_phi\](a-c). The boundary effects are much prominent in these figures and at the boundaries, the deviations from the checker board order of $\phi^k_{p,q}$ are visible without ambiguity. It is to be mentioned here that the domain of the SS phase, for the parameters considered, is rather small. Despite this, SS quantum phase with phase segregation is a novel one and it deserves detailed investigations. On increasing $J/U$ further, we reach the SF phase, which is also phase separated. As an example, the occupancies and SF order parameters for $\mu/U=1.35$ and $J/U=0.015$ are shown in Fig.(\[dw\_ss\_sf\_occ\])(g-i) and Fig.\[ss\_sf\_phi\](d-f), respectively. In the SF phase, there is phase separation, but the occupancies and SF order parameter are uniform within the domains of each species. Thus, the average occupancies and lattice site occupancies are the same $n^1_{p,q, p\leqslant 5}=\rho^1$ and $n^2_{p,q, p>5}=\rho^2$. As we consider identical parameters for both the species $\rho=\rho^1=\rho^2$, where $\rho\in \mathbb{R}$ and $1\leqslant \rho\leqslant 2$. The values and range are also discernible from the figures. The key point from these case studies is that, for non-zero inter-species long range interaction and $U_{12}^2>U_{11}U_{22}$ the eTBHM has quantum phases which are phase separated.

Conclusions {#sec_conclude}
===========

In conclusion, we obtain the phase diagram of TBHM (two-species Bose-Hubbard model) and it’s extended version eTBHM with long range interactions in 2D optical lattices. Our finding are pertinent and timely in view of the recent experimental realization of the Er-Dy binary dipolar BEC mixture [@trautmann_18]. The phase diagram of the TBHM has the unique feature of additional MI lobes with average occupancies which are half integer. These lobes emerge due to the presence of the second species. And, the domain of these lobes are enhanced with the increase of the inter-species interaction strength. In the case of eTBHM, we obtain insulating phases with the non-overlapping density distributions even with $U_{12}^2<U_{11}U_{22}$, where the atoms of the two species are distributed across the system randomly. The non-overlapping densities is like phase separation, but in this work, we use phase separation to mean the configuration where the densities of the two species are segregated into two non-overlapping domains. One key finding of our study is that the DW-MI quantum phase transitions may occur by varying $U_{12}$ while keeping $V_{k}$ fixed. This is in contrast to the single species eBHM, where the NN interaction strength is required to be large to observe such quantum phase transitions. With finite inter-species NN interaction, we obtain the phase diagram in the miscible and immiscible regimes. Our novel result is that the DW, SS and SF phases in the eTBHM in the immiscible domain $U_{12}^2>U_{11}U_{22}$ are phase separated. And, they have side by side order.

Acknowledgements {#sec_acknowledge}
================

The results presented in the paper are based on computations using Vikram-100, the 100TFLOP HPC Cluster at Physical Research Laboratory, Ahmedabad, India. K.S. acknowledges the support of the National Science Centre, Poland via project 2016/21/B/ST2/01086.

Perturbation analysis of the TBH model {#appendix_a}
======================================

The unperturbed ground state at the lattice site $(p,q)$ has the form $|\psi\rangle_{p, q}^{(0)} =|n^1,n^2\rangle_{p, q}$. The energy of this unperturbed ground state is $$\begin{aligned}
    E^{(0)}_{n^{1}_{p,q},n^{2}_{p,q}} 
    &=& \frac{U}{2} \left[n^{1}_{p,q}(n^{1}_{p,q} -1)+n^{2}_{p,q}
    (n^{2}_{p,q} -1)\right]\nonumber\\
    &+& U_{12} n^{1}_{p,q} n^{2}_{p,q} -\mu^{1}_{p,q}n^{1}_{p,q} 
    -\mu^{2}_{p,q}n^{2}_{p,q},
 \label{unperturb_energy}\end{aligned}$$ where we have chosen $U_{11}=U_{22}=U$. Then, to the first order of the SF order parameter $\phi^{k}_{p,q}$ the perturbed ground state can be written as $$\begin{aligned}
    |\psi\rangle &=& |n^1,n^2\rangle_{p,q} \nonumber\\
            &+& \sum_{\substack{m^1,m^2\\ \neq n^1,n^2}} 
            \frac{_{p,q}\langle m^1,m^2| \hat{T} |n^1,n^2 
            \rangle_{p,q}}
            {E^0_{n^1_{p,q},n^2_{p,q}} -E^0_{m^1_{p,q},m^2_{p,q}}} 
            |m^1,m^2\rangle_{p.q},
 \label{perturb_gs}         \end{aligned}$$ where, considering uniform tunnelling strengths for both the species ($J_x^1 =J_x^2 =J_y^1 =J_y^2 =J$) and SF order parameters as real numbers $$\hat{T}_{p, q} = -J \left[ \bar{\phi}_{p, q}^1 \left(
            \hat{b}_{p, q}^{\dagger 1} +\hat{b}_{p, q}^1 \right)
            + \bar{\phi}_{p, q}^2 \left(\hat{b}_{p, q}^{\dagger 2} 
            + \hat{b}_{p, q}^2 \right) \right], 
\label{hop_perturb}$$ with $\bar{\phi}^{k}_{p,q} = \left(\phi^{k}_{p+1,q}+\phi^{k}_{p-1,q}
+\phi^{k}_{p,q+1}+\phi^{k}_{p,q-1}\right)$. Then, using Eqs. (\[unperturb\_energy\])–(\[hop\_perturb\]) the ground state can be calculated as $$\begin{aligned}
    |\psi\rangle_{p,q} =&& |n^1,n^2\rangle_{p,q} \nonumber\\
            +J \bar{\phi}^1_{p,q} 
            &&\left\{\frac{\sqrt{n^1_{p,q} +1}}
            {n^{1}_{p,q}U -\mu^{1}_{p,q} +U_{12}n^{2}_{p,q}}
            |n^{1} +1,n^{2}\rangle_{p,q} \right. \nonumber\\
            &&\left. - \frac{\sqrt{n^1_{p,q}}}{(n^1_{p,q}-1)U 
            -\mu^{1}_{p,q} +U_{12}n^2_{p,q}}
            |n^1 -1,n^2 \rangle_{p,q}\right\} \nonumber\\
            +J \bar{\phi}^2 &&\left\{\frac{\sqrt{n^2_{p,q} +1}}
            {n^2_{p,q}U -\mu^{2}_{p,q} +U_{12}n^1_{p,q}}
            |n^1,n^2 +1\rangle_{p,q} \right. \nonumber \\
            &&\left.-\frac{\sqrt{n^2_{p,q}}}{(n^2_{p,q}-1)U 
            -\mu^{2}_{p,q} +U_{12}n^1_{p,q}}
            |n^1,n^2 -1\rangle_{p,q}\right\}. \nonumber \\
 \label{perturb_psi}\end{aligned}$$ From this state, we obtain the SF order parameter $\phi^{k}_{p,q} = _{p,q}\langle\psi|\hat{b}^{k}_{p,q}|\psi\rangle_{p,q}$, which is mentioned in Eq. (\[order\_par\_ana\]).
