2 A. Ground-State Wave Functions where the pseudo-spin σ ≡ {↑ , ↓} labels the trapped hyper- fine states of a given species of bosons, or labels different types of bosons in a two-species mixture, t$_{ij,σ}$ is the tun- neling (or hopping) matrix between sites i and j , b † i,σ ( b$_{i,σ}$ ) is the boson creation (annihilation) and ̂ n$_{i,σ}$ = b † $_{i,σ}$b$_{i,σ}$ is the boson number operator at site i , U$_{σσ}$ ′ is the strength of the onsite boson-boson interaction between σ and σ ′ compo- nents, and µ$_{σ}$ is the chemical potential. In this manuscript, we consider a d -dimensional hypercubic lattice with M sites, for which we assume t$_{ij,σ}$ is a real symmetric matrix with el- ements t$_{ij,σ}$ = t$_{σ}$ ≥ 0 for i and j nearest neighbors and 0 otherwise. The lattice coordination number (or the number of nearest neighbors) for such lattices is z = 2 d . We take the intraspecies interactions to be repulsive ( { U$_{↑↑}$, U$_{↓↓}$ } > 0 ), but discuss both repulsive and attractive interspecies interaction U$_{↑↓}$ as long as U$_{↑↑}$U$_{↓↓}$ > U 2 $_{↑↓}$. This guarantees the stability of the mixture against collapse when U$_{↑↓}$ ≪ 0 , and against phase separation when U$_{↑↓}$ ≫ 0 . How- ever, when the interspecies interaction is sufficiently lar ge and attractive, we note that instead of a direct transition from the Mott insulator to a single particle superfluid phase, it is po ssi- ble to have a transition from the Mott insulator to a paired su- perfluid phase (superfluidity of composite bosons, i.e. Bose - Bose pairs) [7–16]. Therefore, one needs to consider both possibilities, as discussed next. The perturbation theory is performed with respect to the ground state of the system when t$_{↑}$ = t$_{↓}$ = 0 , and therefore we first need zeroth order wave functions of the Mott phase and of its defect states. To zeroth order in t$_{↑}$ and t$_{↓}$ , the Mott insulator wave function can be written as, | Ψ ins(0) Mott 〉 = 1 √ n$_{↑}$ ! n$_{↓}$ ! ∏ i ( b † i, $_{↑}$) $^{n$_{↑}$}$( b † i, $_{↓}$) $^{n$_{↓}$}$| 0 〉 , (2) where $^{〈}$̂ n$_{i,σ}$ 〉 = n$_{σ}$ is an integer number corresponding to the ground-state occupancy of the pseudo-spin σ bosons, 〈· · · 〉 is the thermal average, and | 0 〉 is the vacuum state. On the other hand, the wave functions of the defect states are determined by degenerate perturbation theory. The reason for that lies in the fact that when exactly one extra elementary particle or hole is added to the Mott phase, it could go to any of the M lattice sites, since all of those states share the same energy when t$_{↑}$ = t$_{↓}$ = 0 . Therefore, the initial degeneracy of the defect states is of order M . When the elementary excitations involve a single- σ -particle (exactly one extra pseudo-spin σ boson) or a single- σ -hole (exactly one less pseudo-spin σ boson), this degeneracy is lifted at first order in t$_{↑}$ and t$_{↓}$ . The treatment for this case is very similar to the single-species BH model [19, 24], and the wave functions (to zeroth order in t$_{↑}$ and t$_{↓}$ ) for the single- σ - particle and single- σ -hole defect states turn out to be | Ψ s σ p(0) def 〉 = 1 √$_{n}$$_{σ}$ + 1 ∑ i f s σ p i b † i,σ | Ψ ins(0) Mott 〉 , (3) | Ψ s σ h(0) def 〉 = 1 √ n$_{σ}$ ∑ i f s σ h i b$_{i,σ}$ | Ψ ins(0) Mott 〉 , (4) where f s σ p i = f s σ h i is the eigenvector of the hopping matrix t$_{ij,σ}$ with the highest eigenvalue (which is zt$_{σ}$ with z = 2 d ) such that ∑ j t$_{ij,σ}$f s σ p j = zt$_{σ}$f s σ p i . The normalization condi- tion requires that ∑ i | f s σ p i | 2 = 1 . Notice that we choose the highest eigenvalue of t$_{ij,σ}$ because the hopping matrix enters the Hamiltonian as − t$_{ij,σ}$ , and we ultimately want the lowest- energy states. We use the many-body version of Rayleigh-Schr¨odinger perturbation theory in the kinetic energy term to perform the expansion (in powers of t$_{↑}$ and t$_{↓}$ ) for the different energies needed to carry out our analysis. The strong-coupling expan- sion technique was previously used to discuss the phase di- agram of the single-species BH model [19–21, 23], extended BH model [24], and of the hardcore BH model with a superlat- tice [25], and its results showed an excellent agreement with Monte Carlo simulations [23, 25]. Motivated by the success of this technique with these models, here we apply it to the two-species BH model. However, when the elementary excitations involve two par- ticles (exactly one extra boson of each species) or two holes (exactly one less boson of each species), the degeneracy is lifted at second order in t$_{↑}$ and t$_{↓}$ . Such elementary excita- tions occur when U$_{↑↓}$ is sufficiently large and attractive [26], and the wave functions (to zeroth order in t$_{↑}$ and t$_{↓}$ ) for the two-particle and two-hole defect states can be written as To determine the phase boundary separating the incom- pressible Mott phase from the compressible superfluid phase within the strong-coupling expansion method, one needs the energy of the Mott phase and of its ‘defect’ states (those states which have exactly one extra elementary particle or hole about the ground state) as a function of t$_{↑}$ and t$_{↓}$ . At the point where the energy of the incompressible state becomes equal to its defect state, the system becomes compressible, assum- ing that the compressibility approaches zero continuously at the phase boundary. Here, we remark that this technique can- not be used to calculate the phase boundary between two com- pressible phases. where f tp i = f th i turns out to be the eigenvector of the t$_{ij,}$$_{↑}$t$_{ij,}$$_{↓}$ matrix with the highest eigenvalue (which is zt$_{↑}$t$_{↓}$ with z = 2 d ) such that ∑ j t$_{ij,}$$_{↑}$t$_{ij,}$$_{↓}$f tp j = zt$_{↑}$t$_{↓}$f tp i . Since the elementary excitations involve two particles or two holes, the III. STRONG-COUPLING EXPANSION | Ψ tp(0) def 〉 = 1 √ ( n$_{↑}$ + 1)( n$_{↓}$ + 1) ∑ i f tp i b † i, $_{↑}$b † i, ↓ | Ψ ins(0) Mott 〉 , (5) th(0) 1 th ins(0) | Ψ th(0) def 〉 = 1 n$_{↑}$n$_{↓}$ ∑ i f th i b$_{i,}$$_{↑}$b$_{i,}$$_{↓}$ | Ψ ins(0) Mott 〉 , (6)