4 This expression is exact for infinite-dimensional hypercub ic lattices, and it recovers the known result for the single species BH model when n$_{−}$$_{σ}$ = 0 or U$_{↑↓}$ = 0 [1]. In the d → ∞ limit (while keeping dt$_{σ}$ constant), we checked that our strong-coupling perturbation results given in Eqs. (8) and (9) agree with this exact solution when the latter is expanded out to third order in t$_{↑}$ and t$_{↓}$ , providing an independent check of the algebra. Equation (10) also shows that, for infinite-dimensional lattices, the M ott lobes are separated by U$_{↑↓}$n$_{−}$$_{σ}$ , but their shapes and critical points (the latter are obtained by setting µ par σ = µ hol σ ) are independent of U$_{↑↓}$ . This is not the case for finite-dimensional lattices as can b e clearly seen from our results. It is also important to mention here that both the shapes and critical points are independent of the sign of U$_{↑↓}$ in finite dimensions (at the third-order presented here) as can be seen in Eqs. (8) and (9). However, when the elementary excitations involve two particles or two holes (which occurs when U$_{↑↓}$ is sufficiently large and attractive [26]), the degeneracy is lifted at second order in t$_{↑}$ and t$_{↓}$ . A lengthy but straightforward calculation leads to the energy of the two-particle defect state up to third order in t$_{↑}$ and t$_{↓}$ as E tp def = E ins $_{Mott}$+ U$_{↑↓}$ ( n$_{↑}$ + n$_{↓}$ + 1) + ∑ σ ( U$_{σσ}$n$_{σ}$ − µ$_{σ}$ ) + 2( n$_{↑}$ + 1)( n$_{↓}$ + 1) U$_{↑↓}$ zt$_{↑}$t$_{↓}$ + ∑ σ [ ( n$_{σ}$ + 1) 2 U$_{↑↓}$ − n$_{σ}$ ( n$_{σ}$ + 2) 2 U$_{σσ}$ + U$_{↑↓}$ + 2 n$_{σ}$ ( n$_{σ}$ + 1) U$_{σσ}$ ] zt 2 $_{σ}$+ O ( t $^{4}$) . (11) E th def = E ins Mott − U$_{↑↓}$ ( n$_{↑}$ + n$_{↓}$ − 1) − ∑ σ [ U$_{σσ}$ ( n$_{σ}$ − 1) − µ$_{σ}$ ] + 2 n$_{↑}$n$_{↓}$ U$_{↑↓}$ zt$_{↑}$t$_{↓}$ + ∑ σ [ n 2 σ U$_{↑↓}$ − ( n 2 σ − 1) 2 U$_{σσ}$ + U$_{↑↓}$ + 2 n$_{σ}$ ( n$_{σ}$ + 1) U$_{σσ}$ ] zt 2 $_{σ}$+ O ( t $^{4}$) , (12) Here, we assume { U$_{σσ}$, | U$_{↑↓}$ | , 2 U$_{σσ}$ + U$_{↑↓}$ } ≫ t$_{σ}$ . Equation (11) is valid for all d -dimensional hypercubic lattices, where the odd-order terms in t$_{↑}$ and t$_{↓}$ vanish [27]. To third order in t$_{↑}$ and t$_{↓}$ , we obtain a similar expression for the energy of the two-hole defect state given by which is also valid for all d -dimensional hypercubic lattices, where the odd-order terms in t$_{↑}$ and t$_{↓}$ vanish [27]. Here, we again assume { U$_{σσ}$, | U$_{↑↓}$ | , 2 U$_{σσ}$ + U$_{↑↓}$ } ≫ t$_{σ}$ . Since the single- σ -particle and single- σ -hole defect states have cor- rections to first order in the hopping, while the two-particl e and two-hole defect states have corrections to second order in the hopping, the slopes of the Mott lobes are finite as { t$_{↑}$, t$_{↓}$ } → 0 in the former case, but they vanish in the lat- ter case. Hence, the shape of the insulating lobes are expected to be very different for two-particle or two-hole excitations. In addition, the chemical-potential widths ( µ$_{σ}$ ) of all Mott lobes are U$_{σσ}$ in the former case, but they [ ( µ$_{↑}$ + µ$_{↓}$ ) / 2 ] are U$_{↑↓}$ + ( U$_{↑↑}$ + U$_{↓↓}$ ) / 2 in the latter. state becomes equal to its defect state, the system becomes compressible, assuming that the compressibility approaches zero continuously at the phase boundary. While E ins Mott and its defects E s σ p def , E s σ h def , E tp def and E th def depend on the lattice size M , their difference do not. Therefore, the chemical po- tentials that determine the particle and hole branches are inde- pendent of M at the phase boundaries. This indicates that the numerical Monte Carlo simulations should not have a strong dependence on M . It is known that the third-order strong-coupling expansion is not very accurate near the tip of the Mott lobes, as t$_{↑}$ and t$_{↓}$ are not very small there [19, 24]. For this reason, an extrapola- tion technique is highly desirable to determine more accurate phase diagrams. Therefore, having discussed the third-order strong-coupling expansion for a general two-species Bose- Bose mixtures with arbitary hoppings t$_{σ}$ , interactions U$_{σσ}$ ′ , densities n$_{σ}$ , and chemical potentials µ$_{σ}$ , next we show how to develop a scaling theory. We note that in the limit when t$_{↑}$ = t$_{↓}$ = t , U$_{↑↑}$ = U$_{↓↓}$ = U$_{0}$ , U$_{↑↓}$ = U $^{′}$, n$_{↑}$ = n$_{↓}$ = n$_{0}$ , µ$_{↑}$ = µ$_{↓}$ = µ , and z = 2 (or d = 1 ), Eq. (12) is in complete agreement with Eq. (3) of Ref. [11], providing an independent check of the algebra. In addition, in the limit when t$_{↑}$ = t$_{↓}$ = J , U$_{↑↑}$ = U$_{↓↓}$ = U , U$_{↑↓}$ = W ≈ − U , n$_{↑}$ = n$_{↓}$ = m , and µ$_{↑}$ = µ$_{↓}$ = µ , Eqs. (11) and (12) reduce to those given in Ref. [12] (after setting U$_{N N}$ = 0 there). However, the terms that are propor- tional to t$_{↑}$t$_{↓}$ are not included in their definitions of the two- particle and two-hole excitation gaps. We also checked the accuracy of Eqs. (11) and (12) via exact small-cluster (two- site) calculations with one particle of each species. We would also like to remark in passing that the energy dif- ference between the Mott phase and its defect states determine the phase boundary of the particle and hole branches. This is because at the point where the energy of the incompressible In this section, we propose a chemical potential extrapo- lation technique based on scaling theory to extrapolate our third-order power-series expansion into a functional form that is appropriate for the entire Mott lobes. It is known that the critical point at the tip of the lobes has the scaling behavior of a ( d +1 )-dimensional XY model, and therefore the lobes have IV. EXTRAPOLATION TECHNIQUE