3 degenerate defect states cannot be connected by one hopping, but rather require two hoppings to be connected. Therefore, one expects the degeneracy to be lifted at least at second order in t$_{↑}$ and t$_{↓}$ , as discussed next. Next, we employ the many-body version of Rayleigh- Schr¨odinger perturbation theory in t$_{↑}$ and t$_{↓}$ with respect to the ground state of the system when t$_{↑}$ = t$_{↓}$ = 0 , and cal- culate the energy of the Mott phase and of its defect states. The energy of the Mott state is obtained via nondegenerate perturbation theory, and to third order in t$_{↑}$ and t$_{↓}$ it is given by This is an extensive quantity, i.e. E ins Mott is proportional to the number of lattice sites M . The odd-order terms in t$_{↑}$ and t$_{↓}$ vanish for the d -dimensional hypercubic lattices considered in this manuscript, which is simply because the Mott state given in Eq. (2) cannot be connected to itself by only one hopping, but rather requires two hoppings to be connected. Notice that Eq. (7) recovers the known result for the single-species BH model when one of the pseudo-spin components have vanish- ing filling, e.g. n$_{↓}$ = 0 [19, 24]. The calculation of the defect-state energies is more involved since it requires using degenerate perturbation theory. As mentioned above, when the elementary excitations involve a single- σ -particle or a single- σ -hole, the degeneracy is lifted at first order in t$_{↑}$ and t$_{↓}$ . A lengthy but straightforward calcula- tion leads to the energy of the single- σ -particle defect state up to third order in t$_{↑}$ and t$_{↓}$ as E ins Mott M = ∑ σ U$_{σσ}$ 2 n$_{σ}$ ( n$_{σ}$ − 1) + U$_{↑↓}$n$_{↑}$n$_{↓}$ − ∑ σ µ$_{σ}$n$_{σ}$ − σ n$_{σ}$ ( n$_{σ}$ + 1) zt 2 σ U$_{σσ}$ + O ( t $^{4}$) . (7) B. Ground-State Energies E s σ p def = E ins Mott + U$_{↑↓}$n$_{−}$$_{σ}$ + U$_{σσ}$n$_{σ}$ − µ$_{σ}$ − ( n$_{σ}$ + 1) zt$_{σ}$ − n$_{σ}$ [ n$_{σ}$ + 2 2 + ( n$_{σ}$ + 1)( z − 3) ] zt 2 σ U$_{σσ}$ − 2 n$_{−}$$_{σ}$ ( n$_{−}$$_{σ}$ + 1) U 2 ↑↓ U 2 − σ − σ − U 2 ↑↓ zt 2 − σ U$_{−}$$_{σ}$$_{−}$$_{σ}$ − n$_{σ}$ ( n$_{σ}$ + 1) [ n$_{σ}$ ( z − 1) 2 + ( n$_{σ}$ + 1)( z − 1)( z − 4) + ( n$_{σ}$ + 2)(3 z/ 4 − 1) ] zt 3 σ U 2 σσ − 4( n$_{σ}$ + 1) n$_{−}$$_{σ}$ ( n$_{−}$$_{σ}$ + 1) U 2 ↑↓ U 2 − σ − σ − U 2 ↑↓ ( z − 1 − U 2 − σ − σ U 2 − σ − σ − U 2 ↑↓ ) zt$_{σ}$t 2 − σ U 2 − σ − σ + O ( t $^{4}$) , (8) E s σ h def = E ins Mott − U$_{↑↓}$n$_{−}$$_{σ}$ − U$_{σσ}$ ( n$_{σ}$ − 1) + µ$_{σ}$ − n$_{σ}$zt$_{σ}$ − ( n$_{σ}$ + 1) [ n$_{σ}$ − 1 2 + n$_{σ}$ ( z − 3) ] zt 2 σ U$_{σσ}$ − 2 n$_{−}$$_{σ}$ ( n$_{−}$$_{σ}$ + 1) U 2 ↑↓ U 2 − σ − σ − U 2 ↑↓ zt 2 − σ U$_{−}$$_{σ}$$_{−}$$_{σ}$ − n$_{σ}$ ( n$_{σ}$ + 1) [ ( n$_{σ}$ + 1)( z − 1) 2 + n$_{σ}$ ( z − 1)( z − 4) + ( n$_{σ}$ − 1)(3 z/ 4 − 1) ] zt 3 σ U 2 σσ − 4 n$_{σ}$n$_{−}$$_{σ}$ ( n$_{−}$$_{σ}$ + 1) U 2 ↑↓ U 2 − σ − σ − U 2 ↑↓ ( z − 1 − U 2 − σ − σ U 2 − σ − σ − U 2 ↑↓ ) zt$_{σ}$t 2 − σ U 2 − σ − σ + O ( t $^{4}$) , (9) where ( − ↑ ) ≡↓ and vice versa. Here, we assume U$_{σσ}$ ≫ t$_{σ}$ and { U$_{−}$$_{σ}$$_{−}$$_{σ}$, | U$_{−}$$_{σ}$$_{−}$$_{σ}$ ± U$_{↑↓}$ |} ≫ t$_{−}$$_{σ}$ . Equation (8) is valid for all d -dimensional hypercubic lattices, and it recovers the known result for the single species BH model when n$_{−}$$_{σ}$ = 0 [19, 24]. Note that this expression also recovers the known result for the single species BH model when U$_{↑↓}$ = 0 , which provides an independent check of the algebra. To third order in t$_{↑}$ and t$_{↓}$ , we obtain a similar expression for the energy of the single- σ -hole defect state given by which is also valid for all d -dimensional hypercubic lattices, and it also recovers the known result for the single-species BH model when n$_{−}$$_{σ}$ = 0 or U$_{↑↓}$ = 0 [19, 24]. Here, we again assume U$_{σσ}$ ≫ t$_{σ}$ and { U$_{−}$$_{σ}$$_{−}$$_{σ}$, | U$_{−}$$_{σ}$$_{−}$$_{σ}$ ± U$_{↑↓}$ |} ≫ t$_{−}$$_{σ}$ . We also checked the accuracy of the second-order terms in Eqs. (8) and (9) via exact small-cluster (two-site) calculations with one σ and two − σ particles. We note that the mean-field phase boundary between the Mott ph ase and its single- σ -particle and single- σ -hole defect states can be calculated as µ par , hol σ = U$_{σσ}$ ( n$_{σ}$ − 1 / 2) + U$_{↑↓}$n$_{−}$$_{σ}$ − zt$_{σ}$/ 2 ± √ U 2 $_{σσ}$/ 4 − U$_{σσ}$ ( n$_{σ}$ + 1 / 2) zt$_{σ}$ + z $^{2}$t2 $_{σ}$/ 4 . (10)