Quantum phase transitions and the degree of nonidentity in the system with two different species of vector bosons

Experimental research of ultracold atoms in optical lattices have dramatically expanded the possibilities of a tunable simulation in quantum many-body physics [1–14]. Moreover, ultracold atoms open the path to the parameter range that is hardly possible or even impossible to achieve in the natural condensed matter systems [15–21].

The typical example is the system of vector bosons. This case corresponds to Bose–Hubbard model that is absent in the standard solid state theory and is the topic of intensive research last time [22–34], and [35–37] for review. Numerous striking effects, in particular, induced by the multispecie nature of boson system have been found recently including quantum phase transitions, many-body localization, and topological order, as well as the superfluidity and supersolidity of ultracold atomic systems [35–37]. Here, we focus on still unexplored physical phenomena in the systems with different species of vector bosons originating from the tunable interplay of spin degrees of freedom and of those identifying different sorts of atoms. Then the situation becomes quite intriguing: in addition to the standard parameters, there appear nontrivial ones related to the 'degree of atomic nonidentity': the difference of tunneling amplitudes and on-site interactions. We address here new physical effects including quantum phase transitions driven by atom nonidentity.

Vector cold atoms in optical lattices are characterized in general by the following parameters: intersite hopping amplitudes t_α, where α = 1, 2 labels different atoms, on-site ${U}_{\alpha ,\alpha ^{\prime} }$ Coulomb interactions, and spin channel interaction parameters U_s [31, 38].

In our recent paper [31], we have considered the simplest limiting case of nearly identical bosons in the Mott insulating state: U₁₂ ≃ U₁₁ ≃ U₂₂ = U₀ and t₁ ≃ t₂ ≪ U₀. This model differs from the case of perfectly identical bosons just by the absence of cross-tunneling term: tunneling with the change of boson identity was forbidden. It has been shown in [31] that such model can be reduced to the Kugel–Khomskii [39] type spin-pseudospin model, but in contrast to the usual version of the latter model, here spin is S = 1 and pseudospin-1/2 (where pseudospin labels different bosons). New effects in [31] have been found, in particular, that the ground state of the system always has the pseudospin domain structure while sign change of U_s switches the spin arrangement of the ground state within domains from ferro to aniferromagnetic (AFM).

Here we also consider two species of vector bosons in an optical lattice, but we get rid of the identity conditions U₁₂ ≃ U₁₁ ≃ U₂₂ = U₀, t₁ ≃ t₂ and trace the evolution of the system quantum state with the increase in the degree of atomic nonidentity (β = t₂/t₁, ξ₁₂ = U₁₂/U₁₁ and ξ₂₂ = U₂₂/U₁₁) starting from nearly identical atoms of [31].

Surprisingly, instead of some naive and predictable modification of our previous results [31], we find that the boson system evolves between six quantum states instead of two in [31] and the manifold of quantum phase transitions is so nontrivial that one can even compare the resulting phase diagrams with modern art paintings (we find some similarity with 'J Miró-like paintings' or 'K Malewicz-like' paintings depending on atom nonidentity), see figures 1–9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The rest of our paper is organized as follows: in section 2, we introduce the model Hamiltonian, then in section 3, we reduce the initial general Hamiltonian for vector bosons to the effective Hamiltonian appearing to be anisotropic spin-pseudospin model of the Kugel–Khomskii type [39]; in section 4, we investigate different possible configurations in spin and pseudospin spaces and find their energy; in section 5, we address the ground state energy of the system and the quantum phase transitions. In particular, in section 5, we illustrate the evolution of phase diagrams with the degree of atomic nonidentity, which looks like the transformation from the Joan Miró style artistic image to that of Kazimir Malewicz. Finally in the appendix, we present the analysis of several special limiting cases of the model Hamiltonian that relate our model system to some well-known results.

We consider two types of boson atoms with S = 1 in the optical lattice with sites labeled by index i. The corresponding creation operators ${c}_{i\alpha \sigma }^{\dagger }$ where σ = {−1, 0, 1} is the spin index and α = 1, 2 labels the type of boson.

The Hamiltonian includes three terms:

$$\begin{eqnarray}&&H={H}^{({U}_{0})}+{H}^{({U}_{s})}+{H}_{t}.\end{eqnarray} \tag{ 1 }$$

The interaction between bosons is given by two terms:

$$\begin{eqnarray}&&{H}^{({U}_{0})}=\displaystyle \sum _{i}{U}_{12}{n}_{i,1}{n}_{i,2}+\displaystyle \frac{1}{2}\displaystyle \sum _{i,\alpha =1,2}{U}_{\alpha \alpha }{n}_{i,\alpha }({n}_{i,\alpha }-1).\end{eqnarray} \tag{ 2 }$$

These terms corresponds to the repulsion between boson atoms at the same site [36]. Here, U₁₁, U₂₂, and U₁₂ are three interaction parameters and ${n}_{i\alpha }={\sum }_{\sigma }{c}_{i\alpha \sigma }^{\dagger }{c}_{i\alpha \sigma }$. Now we assume that the interaction constants can strongly differ from each other, in contrast to the case considered in [31].

The spin-dependent interaction term is taken in the standard form [40, 41]:

$$\begin{eqnarray}&&{H}_{i}^{({U}_{s})}={U}_{s}({{\bf{S}}}_{i}^{2}-2{n}_{i})/2,\end{eqnarray} \tag{ 3 }$$

where n_i = n_i,1 + n_i,2 is the total number of bosons at site i.

The hopping term

$$\begin{eqnarray}&&{H}_{t}=\displaystyle \sum _{\langle i,j\rangle }{h}_{t}\equiv \displaystyle \sum _{\langle i,j\rangle }{t}_{\alpha }({c}_{i\alpha \sigma }^{\dagger }{c}_{j\alpha \sigma }+{c}_{j\alpha \sigma }^{\dagger }{c}_{i\alpha \sigma }),\end{eqnarray} \tag{ 4 }$$

where $\langle i,j\rangle$ means the summation only over the nearest-neighbor sites and h_t is the hopping Hamiltonian for one link. As usual, the repeated indices imply summation.

Below we focus on the Mott insulating state, where cold atoms are localized at the sites of optical lattice with the number of bosons at each site equal to unity. We remind that in such case, the hopping terms (4) can be treated as perturbation compared to the interaction part of the Hamiltonian, ${H}^{({U}_{0})}+{H}^{({U}_{s})}$. Application of the perturbation theory reduces the initial general Hamiltonian to a simpler effective Hamiltonian written solely in terms of the spin and pseudospin (p-spin) operators related to the lattice sites with atom filling equal to one. These spin-1 S_i and p-spin-1/2 ${{\boldsymbol{ \mathcal T }}}_{i}$ operators are defined in a standard way [42]

$$\begin{eqnarray}&&{S}_{i}^{a}={c}_{i\alpha \sigma }^{\dagger }{s}_{\sigma \sigma ^{\prime} }^{a}{c}_{i\alpha \sigma ^{\prime} },\qquad {{ \mathcal T }}_{i}^{a}={c}_{i\alpha \sigma }^{\dagger }{\tau }_{\alpha \beta }^{a}{c}_{i\beta \sigma },\end{eqnarray} \tag{ 5 }$$

where a = x, y, z.

Below we outline the algorithm of transforming the initial Hamiltonian to the effective one.

3.1. Basis states

In what follows, when we consider the link $\langle i,j\rangle$ between the nearest-neighbor sites, we focus on the basis of possible states for two bosons with spins S₁ = 1 and S₂ = 1 at neighboring sites i = 1 and j = 2. We are interested in the case with single occupation, i.e. when one boson of either type is located at each lattice site, ${n}_{i1}+{n}_{i2}=1$. We can pass now to the basis of the eigenstates of the total spin squared ${{\bf{S}}}^{2}={({{\bf{S}}}_{1}+{{\bf{S}}}_{2})}^{2}$ and its z-projection ${S}^{z}={S}_{1}^{z}+{S}_{2}^{z}$. We designate these states as $| {SM}\rangle ,S=0,1,2$ and M = −S, ..., S. This basis can be written as follows

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{SM}}^{(f)}=| {\phi }_{S}^{(f)}\rangle | {SM}\rangle ,\end{eqnarray} \tag{ 6 }$$

where f = 1, ..., 4 enumerates the ways to distribute two types of bosons over two sites. The coordinate part $| {\phi }_{S}^{(f)}\rangle$ is given explicitly in [31].

Applying h_t, see equation (4), to the basis states (6), we obtain two kinds of intermediate (virtual) states. The first type will be realized for two identical bosons at one site (i or j), the second type is for two nonidentical bosons at one site. Intermediate energies depend on the spin and types of bosons. They are

$$\begin{eqnarray}{E}_{S=0}^{{aa}} & = & {U}_{11}-2{U}_{s},\quad {E}_{S=2}^{{aa}}={U}_{11}+{U}_{s},\\ {E}_{S=0}^{{bb}} & = & {U}_{22}-2{U}_{s},\quad {E}_{S=2}^{{bb}}={U}_{22}+{U}_{s},\\ {E}_{S=0}^{{ab}} & = & {U}_{12}-2{U}_{s},\quad {E}_{S=2}^{{ab}}={U}_{12}+{U}_{s},\\ {E}_{S=1}^{{ab}} & = & {U}_{12}-{U}_{s}.\end{eqnarray} \tag{ 7 }$$

There are no intermediate states corresponding to the total spin S = 1 for two identical a- or b-bosons due to the symmetry of the total wave function [40].

The energy of the virtual states with the double site occupancy is much larger than the energy of the states with the single site occupancy we are focusing at. To find corrections to the energy of the single occupancy states related to the hoppings, we need the second-order terms of the perturbation theory. So, further it will be convenient to work with the operator

$$\begin{eqnarray}&&h=-{h}_{\mathrm{eff}}={h}_{t}(1/{H}_{0}){h}_{t},\end{eqnarray} \tag{ 8 }$$

where ${H}_{0}={H}^{({U}_{0})}+{H}^{({U}_{s})}$. In the basis of states (6), the matrix of h can be presented in the following block form

$$\begin{eqnarray}h=\left(\begin{array}{cccc}{B}_{11} & 0 & 0 & 0\\ 0 & {B}_{22} & {B}_{23} & 0\\ 0 & {B}_{32} & {B}_{33} & 0\\ 0 & 0 & 0 & {B}_{44}\end{array}\right).\end{eqnarray} \tag{ 9 }$$

The matrix h here is in fact the block matrix, where each block is a 9 × 9 matrix. Blocks ${B}_{{kl}}=\langle {{\rm{\Phi }}}_{S^{\prime} M^{\prime} }^{k}| h| {{\rm{\Phi }}}_{{SM}}^{l}\rangle ,k,l=1,2,3,4$ are diagonal matrices. Their explicit forms are the following

$$\begin{eqnarray}{B}_{11}=4{t}_{1}^{2}\left(\begin{array}{cccc}\tfrac{1}{{E}_{S=0}^{{aa}}}{I}_{1} & & & \\ & 0\cdot {I}_{3} & & \\ & & \tfrac{1}{{E}_{S=2}^{{aa}}}{I}_{5} & \end{array}\right),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}{B}_{44}=4{t}_{2}^{2}\left(\begin{array}{cccc}\tfrac{1}{{E}_{S=0}^{{bb}}}{I}_{1} & & & \\ & 0\cdot {I}_{3} & & \\ & & \tfrac{1}{{E}_{S=2}^{{bb}}}{I}_{5} & \end{array}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}{B}_{22}={B}_{33}=({t}_{1}^{2}+{t}_{2}^{2})\left(\begin{array}{cccc}\displaystyle \frac{1}{{E}_{S=0}^{{ab}}}{I}_{1} & & & \\ & \tfrac{1}{{E}_{S=1}^{{ab}}}{I}_{3} & & \\ & & \tfrac{1}{{E}_{S=2}^{{ab}}}{I}_{5} & \end{array}\right),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{B}_{23}={B}_{32}=2{t}_{1}{t}_{2}\left(\begin{array}{cccc}\tfrac{1}{{E}_{S=0}^{{ab}}}{I}_{1} & & & \\ & \tfrac{1}{{E}_{S=1}^{{ab}}}{I}_{3} & & \\ & & \tfrac{1}{{E}_{S=2}^{{ab}}}{I}_{5} & \end{array}\right).\end{eqnarray} \tag{ 13 }$$

Here I_n, n = 1, 3, 5, are the identity n × n matrices. I₁ accounts for one state with S = 0, I₃ accounts for three states with S = 1, and I₅ accounts for five states with S = 2.

In what follows, we identify the single occupancy of site i with a- or b-boson by the pseudospin-1/2 states $| 1{\rangle }_{i}=| +{\rangle }_{i}$ and $| 2{\rangle }_{i}=| -{\rangle }_{i}$, respectively. It will be convenient to introduce explicitly two types of creation operators: ${a}_{{is}}^{\dagger }$ for bosons of type α = 1 and ${b}_{{is}}^{\dagger }$ for bosons of type α = 2.

For this purpose, we rewrite the pseudospin operator ${{ \mathcal T }}_{i}^{\gamma }$ at sites i, see equation (5), in the form

$$\begin{eqnarray}&&{{ \mathcal T }}_{i}^{\gamma }={a}_{{is}}^{\dagger }{\tau }_{11}^{\gamma }{a}_{{is}}+{a}_{{is}}^{\dagger }{\tau }_{12}^{\gamma }{b}_{{is}}+{b}_{{is}}^{\dagger }{\tau }_{21}^{\gamma }{a}_{{is}}+{b}_{{is}}^{\dagger }{\tau }_{22}^{\gamma }{b}_{{is}}.\end{eqnarray} \tag{ 14 }$$

We can rewrite, as usual, the set of ${{ \mathcal T }}^{\gamma }$ operators in the other equivalent form:

$$\begin{eqnarray}&&{T}_{i}^{+}={a}_{{is}}^{\dagger }{b}_{{is}},\quad {T}_{i}^{-}={b}_{{is}}^{\dagger }{a}_{{is}},\quad {T}_{i}^{z}=\displaystyle \frac{1}{2}({a}_{{is}}^{\dagger }{a}_{{is}}-{b}_{{is}}^{\dagger }{b}_{{is}}).\end{eqnarray} \tag{ 15 }$$

To describe the occupancy of sites i and j, we introduce the basis of pseudospin states $| \alpha \beta \rangle =| \alpha {\rangle }_{i}| \beta {\rangle }_{j}$. Then, we find the correspondence between two-boson orbital states (6) and pseudospin states $| \alpha \beta \rangle$. For example, $| {\phi }_{S}^{(1)}\rangle \longrightarrow | ++\rangle$.

In what follows, we map the matrix h, equation (9), onto an effective spin-pseudospin operator in the space $| \alpha {\rangle }_{i}| \beta {\rangle }_{j}| {SM}\rangle$. This operator will be given in terms of spin S = 1 operators S_i, S_j and pseudospin T = 1/2 operators T_i, T_j and it has the same structure as matrix h, equation (9).

Next, we introduce the projection operator Q_S and P_T onto the combination of states $| {SM}\rangle$ and $| {{TM}}_{T}\rangle$ corresponding to the total spin S = 0, 1, 2 and pseudospin T = 0, 1 at the link $\langle i,j\rangle$. The projectors in the spin space ${Q}_{S}={\sum }_{M=-S}^{S}| {SM}\rangle \langle {SM}|$ can be written as

$$\begin{eqnarray}{Q}_{0} & = & -\displaystyle \frac{1}{3}+\displaystyle \frac{1}{3}{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2},\\ {Q}_{1} & = & 1-\displaystyle \frac{1}{2}({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})-\displaystyle \frac{1}{2}{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2},\\ {Q}_{2} & = & \displaystyle \frac{1}{3}+\displaystyle \frac{1}{2}({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})+\displaystyle \frac{1}{6}{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2},\end{eqnarray} \tag{ 16 }$$

where Q₀ + Q₁ + Q₂ = 1.

Similarly, in the pseudospin space the projectors onto the singlet T = 0 and triplet T = 1 states are

$$\begin{eqnarray}&&{P}_{s}=\displaystyle \frac{1}{4}-{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j},\quad {P}_{t}=\displaystyle \frac{3}{4}+{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}.\end{eqnarray} \tag{ 17 }$$

It is also convenient to introduce the following projectors in the pseudospin space

$$\begin{eqnarray}{P}^{11} & = & | ++\rangle \langle ++| =\left(\displaystyle \frac{1}{2}+{T}_{i}^{z}\right)\left(\displaystyle \frac{1}{2}+{T}_{j}^{z}\right),\\ {P}^{22} & = & | +-\rangle \langle +-| =\left(\displaystyle \frac{1}{2}+{T}_{i}^{z}\right)\left(\displaystyle \frac{1}{2}-{T}_{j}^{z}\right),\\ {P}^{33} & = & | -+\rangle \langle -+| =\left(\displaystyle \frac{1}{2}-{T}_{i}^{z}\right)\left(\displaystyle \frac{1}{2}+{T}_{j}^{z}\right),\\ {P}^{44} & = & | --\rangle \langle --| =\left(\displaystyle \frac{1}{2}-{T}_{i}^{z}\right)\left(\displaystyle \frac{1}{2}-{T}_{j}^{z}\right),\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}{P}^{32} & = & | -+\rangle \langle +-| ={T}_{i}^{-}{T}_{j}^{+},\\ {P}^{23} & = & | +-\rangle \langle -+| ={T}_{i}^{+}{T}_{j}^{-}.\end{eqnarray} \tag{ 19 }$$

In addition, we use below the following identities:

$$\begin{eqnarray}{P}^{11}+{P}^{44} & = & \displaystyle \frac{1}{2}+2{T}_{i}^{z}{T}_{j}^{z},\\ {P}^{22}+{P}^{33} & = & \displaystyle \frac{1}{2}-2{T}_{i}^{z}{T}_{j}^{z},\\ {P}^{32}+{P}^{23} & = & 2{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}-2{T}_{i}^{z}{T}_{j}^{z}.\end{eqnarray} \tag{ 20 }$$

With the help of projectors (16)–(19) and identities (20), we rewrite the block matrix h in terms of spin and pseudospin operators as follows

$$\begin{eqnarray}h & = & 4{t}_{1}^{2}\,{P}^{11}\left[\displaystyle \frac{{Q}_{0}}{{E}_{0}^{{aa}}}+\displaystyle \frac{{Q}_{2}}{{E}_{2}^{{aa}}}\right]+4{t}_{2}^{2}\,{P}^{44}\left[\displaystyle \frac{{Q}_{0}}{{E}_{0}^{{bb}}}+\displaystyle \frac{{Q}_{2}}{{E}_{2}^{{bb}}}\right]+({t}_{1}^{2}+{t}_{2}^{2})[{P}^{22}+{P}^{33}]\left[\displaystyle \frac{{Q}_{0}}{{E}_{0}^{{ab}}}+\displaystyle \frac{{Q}_{1}}{{E}_{1}^{{ab}}}+\displaystyle \frac{{Q}_{2}}{{E}_{2}^{{ab}}}\right]\\ & & +2{t}_{1}{t}_{2}[{P}^{23}+{P}^{32}]\left[\displaystyle \frac{{Q}_{0}}{{E}_{0}^{{ab}}}+\displaystyle \frac{{Q}_{1}}{{E}_{1}^{{ab}}}+\displaystyle \frac{{Q}_{2}}{{E}_{2}^{{ab}}}\right].\end{eqnarray} \tag{ 21 }$$

Substituting the explicit form of pseudospin projectors, we can rewrite h in the form

$$\begin{eqnarray}h & = & \left\{\displaystyle \frac{{t}_{1}^{2}}{{E}_{0}^{{aa}}}+\displaystyle \frac{{t}_{2}^{2}}{{E}_{0}^{{bb}}}+2\left(\displaystyle \frac{{t}_{1}^{2}}{{E}_{0}^{{aa}}}-\displaystyle \frac{{t}_{2}^{2}}{{E}_{0}^{{bb}}}\right)({T}_{i}^{z}+{T}_{j}^{z})+4\left(\displaystyle \frac{{t}_{1}^{2}}{{E}_{0}^{{aa}}}+\displaystyle \frac{{t}_{2}^{2}}{{E}_{0}^{{bb}}}\right){T}_{i}^{z}{T}_{j}^{z}\right\}{Q}_{0}\\ & & +\left\{\displaystyle \frac{{t}_{1}^{2}}{{E}_{2}^{{aa}}}+\displaystyle \frac{{t}_{2}^{2}}{{E}_{2}^{{bb}}}+2\left(\displaystyle \frac{{t}_{1}^{2}}{{E}_{2}^{{aa}}}-\displaystyle \frac{{t}_{2}^{2}}{{E}_{2}^{{bb}}}\right)({T}_{i}^{z}+{T}_{j}^{z})+4\left(\displaystyle \frac{{t}_{1}^{2}}{{E}_{2}^{{aa}}}+\displaystyle \frac{{t}_{2}^{2}}{{E}_{2}^{{bb}}}\right){T}_{i}^{z}{T}_{j}^{z}\right\}{Q}_{2}\\ & & +\left\{\displaystyle \frac{1}{2}({t}_{1}^{2}+{t}_{2}^{2})-2{({t}_{1}+{t}_{2})}^{2}{T}_{i}^{z}{T}_{j}^{z}+4{t}_{1}{t}_{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}\right\}\left[\displaystyle \frac{{Q}_{0}}{{E}_{0}^{{ab}}}+\displaystyle \frac{{Q}_{1}}{{E}_{1}^{{ab}}}+\displaystyle \frac{{Q}_{2}}{{E}_{2}^{{ab}}}\right].\end{eqnarray} \tag{ 22 }$$

Finally, the effective Hamiltonian is written as the sum of h operators (22) over all $\langle i,j\rangle$ links

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}=-\sum _{\langle i,j\rangle }h.\end{eqnarray} \tag{ 23 }$$

This is the most general form of the effective Hamiltonian involving different Hubbard interaction parameters, different hopping amplitudes, and spin-dependent interaction. In the appendix, we show that in a number of limiting cases, this Hamiltonian can be simplified to some well-known forms.

To proceed with the calculation of the ground state energy, we rewrite the 'kernel' h of effective Hamiltonian (23) as follows

$$\begin{eqnarray}&&h=\displaystyle \frac{2{t}_{1}^{2}}{{U}_{11}}\left({R}_{0}{Q}_{0}+{R}_{2}{Q}_{2}+{R}_{1}\left[\displaystyle \frac{{Q}_{0}}{{E}_{0}^{{ab}}}+\displaystyle \frac{{Q}_{1}}{{E}_{1}^{{ab}}}+\displaystyle \frac{{Q}_{2}}{{E}_{2}^{{ab}}}\right]\right),\end{eqnarray} \tag{ 24 }$$

where the coefficients R₀, R₁ and R₂ are

$$\begin{eqnarray}&&{R}_{0}=\left[\displaystyle \frac{1}{2}\left(\displaystyle \frac{1}{{E}_{0}^{{aa}}}+\displaystyle \frac{{\beta }^{2}}{{E}_{0}^{{bb}}}\right)+\left(\displaystyle \frac{1}{{E}_{0}^{{aa}}}-\displaystyle \frac{{\beta }^{2}}{{E}_{0}^{{bb}}}\right)({T}_{i}^{z}+{T}_{j}^{z})+2\left(\displaystyle \frac{1}{{E}_{0}^{{aa}}}+\displaystyle \frac{{\beta }^{2}}{{E}_{0}^{{bb}}}\right)({T}_{i}^{z}{T}_{j}^{z})\right],\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{R}_{1}=\left[\displaystyle \frac{1}{4}(1+{\beta }^{2})-{(1+\beta )}^{2}({T}_{i}^{z}{T}_{j}^{z})+2\beta ({{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j})\right],\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{R}_{2}=\left[\displaystyle \frac{1}{2}\left(\displaystyle \frac{1}{{E}_{2}^{{aa}}}+\displaystyle \frac{{\beta }^{2}}{{E}_{2}^{{bb}}}\right)+\left(\displaystyle \frac{1}{{E}_{2}^{{aa}}}-\displaystyle \frac{{\beta }^{2}}{{E}_{2}^{{bb}}}\right)({T}_{i}^{z}+{T}_{j}^{z})+2\left(\displaystyle \frac{1}{{E}_{2}^{{aa}}}+\displaystyle \frac{{\beta }^{2}}{{E}_{2}^{{bb}}}\right)({T}_{i}^{z}{T}_{j}^{z})\right].\end{eqnarray} \tag{ 27 }$$

Here, we introduce the dimensionless parameter β = t₂/t₁ that characterizes the difference of the tunnel amplitudes for different boson species.

The energy E of the ground state can be formally written as the average of the effective Hamiltonian over the ground state wave function (we will find it later using the variational approach)

$$\begin{eqnarray}&&E=-\displaystyle \frac{\nu }{2}\langle h\rangle ,\end{eqnarray} \tag{ 28 }$$

where $\nu =2D$ is the number of nearest neighbors for the D-dimensional cubic lattice.

Within the mean-field approximation, we can neglect any correlations between spin and p-spin degrees of freedom. Then, for example, $\langle {R}_{0}{Q}_{0}\rangle \to \langle {R}_{0}\rangle \langle {Q}_{0}\rangle$, and

$$\begin{eqnarray}&&\langle h\rangle ={E}_{u}\left(\langle {R}_{0}\rangle \langle {Q}_{0}\rangle +\langle {R}_{2}\rangle \langle {Q}_{2}\rangle +\langle {R}_{1}\rangle \left[\displaystyle \frac{\langle {Q}_{0}\rangle }{{E}_{0}^{{ab}}}+\displaystyle \frac{\langle {Q}_{1}\rangle }{{E}_{1}^{{ab}}}+\displaystyle \frac{\langle {Q}_{2}\rangle }{{E}_{2}^{{ab}}}\right]\right),\end{eqnarray} \tag{ 29 }$$

where ${E}_{u}=\tfrac{2{t}_{1}^{2}}{{U}_{11}}$.

We take the trial wave function in the p-spin space as the 'mixed orbital state' on two-sublattices A and B

$$\begin{eqnarray}&&| \chi {\rangle }_{i}=\cos \theta | +{\rangle }_{i}+\sin \theta | -{\rangle }_{i},\quad i\in A,\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&| \chi {\rangle }_{j}=\cos \theta | +{\rangle }_{j}+\eta \sin \theta | -{\rangle }_{j},\quad j\in B,\quad \eta =\pm 1.\end{eqnarray} \tag{ 31 }$$

For θ = 0, it defines the p-spin ferromagnetic (FM) state

$$\begin{eqnarray}&&| {\uparrow }_{i}{\uparrow }_{j}{\uparrow }_{i}{\uparrow }_{j}...\rangle .\end{eqnarray} \tag{ 32 }$$

For $\theta =\tfrac{\pi }{2}$, it defines also the p-spin FM state

$$\begin{eqnarray}&&| {\downarrow }_{i}(-{\downarrow }_{j}){\downarrow }_{i}(-{\downarrow }_{j})\,...\rangle .\end{eqnarray} \tag{ 33 }$$

The FM and AFM p-spin states will occur at $\theta =\tfrac{\pi }{4}$ for η = +1 and η = −1, respectively. These states are polarized in x direction:

$$\begin{eqnarray}&&\displaystyle \frac{1}{\sqrt{2}}{(| \uparrow \rangle +| \downarrow \rangle )}_{i},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{\sqrt{2}}{(| \uparrow \rangle +\eta | \downarrow \rangle )}_{j}.\end{eqnarray} \tag{ 35 }$$

For other angles θ, the mixed orbital state is a some intermediate state between the FM and AFM ones. That is why, it is more correctly to refer to this mixed orbital state as the mixed orbital state with η = +1, or η = −1 rather than the FM and AFM states.

Now, we find averages over the trial mixed orbital state of the p-spin operators entering R_i, i = 0, 1, 2

$$\begin{eqnarray}\langle {T}_{i}^{z}+{T}_{j}^{z}\rangle =\left\{\begin{array}{cc}\cos 2\theta ,\quad & \eta =+1,\\ 0,\quad & \eta =-1,\end{array}\right.\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}\langle {T}_{i}^{z}{T}_{j}^{z}\rangle =\left\{\begin{array}{cc}\tfrac{1}{4}{\cos }^{2}2\theta ,\quad & \eta =+1,\\ -\tfrac{1}{4}{\cos }^{2}2\theta ,\quad & \eta =-1,\end{array}\right.\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}\langle {{\bf{T}}}_{i}{{\bf{T}}}_{j}\rangle =\left\{\begin{array}{cc}\tfrac{1}{4},\quad & \eta =+1,\\ -\tfrac{1}{4},\quad & \eta =-1.\end{array}\right.\end{eqnarray} \tag{ 38 }$$

Later on, we will minimize the energy E with respect to the angle θ.

In contrast to the p-spin space, the effective Hamiltonian in the spin space is isotropic. So, we do not expect any exotic states there. Hence, we take for trial wave functions in the spin space the usual FM, AFM, and nematic (NEM) states. It should be noted that NEM naturally appears for spin-1 isotropic Heisenberg Hamiltonian, it has zero average site spin, and it is discussed in detail in [40].

Then, the coefficients for projectors in the spin space are

$$\begin{eqnarray}\langle {Q}_{0}\rangle & = & \left\{\begin{array}{cc}0, & \quad \mathrm{FM},\\ \tfrac{1}{3}, & \quad \mathrm{AFM},\\ \tfrac{1}{3}, & \quad \mathrm{NEM},\end{array}\right.\qquad \langle {Q}_{1}\rangle =\left\{\begin{array}{cc}0, & \quad \mathrm{FM},\\ \tfrac{1}{2}, & \quad \mathrm{AFM},\\ 0, & \quad \mathrm{NEM},\end{array}\right.\\ \langle {Q}_{2}\rangle & = & \left\{\begin{array}{cc}1, & \quad \mathrm{FM},\\ \tfrac{1}{6}, & \quad \mathrm{AFM},\\ \tfrac{2}{3}, & \quad \mathrm{NEM}.\end{array}\right.\end{eqnarray} \tag{ 39 }$$

Finally, we have found all the averages and correlation functions entering equation (28) for the energy. We minimize numerically the energy E(θ, η), find the ground stats at different values of the parameters, and draw the corresponding phase diagrams.

5.1. General properties of the phase diagrams

We address here the ground state of the system and possible quantum phase transitions. Performing the energy minimization, we find the phase diagrams for different ranges of the parameters. The key parameters are β = t₂/t₁, λ = U_s/U₁₁, ξ₁₂ = U₁₂/U₁₁, and ξ₂₂ = U₂₂/U₁₁. All these parameters (except λ responsible for spin channel interaction) are related to the difference between the types of bosons.

There are six different phases: three phases have the FM pseudospin arrangement, while three others correspond to the AFM pseudospin state. Evolution of these phases is illustrated in figures 1–9; figure 1 is supplemented by the colorbar, where the correspondence between colors and phases is shown. Points at the colorbar correspond to the following orders in spin and pseudospin systems:

$$\begin{eqnarray}1 & \to & {E}_{\mathrm{FM}}^{{\rm{p}} \mbox{-} \mathrm{spin}\,\mathrm{FM}},\\ 2 & \to & {E}_{\mathrm{AFM}}^{{\rm{p}} \mbox{-} \mathrm{spin}\,\mathrm{FM}},\\ 3 & \to & {E}_{\mathrm{NEM}}^{{\rm{p}} \mbox{-} \mathrm{spin}\,\mathrm{FM}},\\ 4 & \to & {E}_{\mathrm{FM}}^{{\rm{p}} \mbox{-} \mathrm{spin}\,\mathrm{AFM}},\\ 5 & \to & {E}_{\mathrm{AFM}}^{{\rm{p}} \mbox{-} \mathrm{spin}\,\mathrm{AFM}},\\ 6 & \to & {E}_{\mathrm{NEM}}^{{\rm{p}} \mbox{-} \mathrm{spin}\,\mathrm{AFM}}.\end{eqnarray} \tag{ 40 }$$

Here, for example, color '1' corresponds to FM spin and pseudospin orders.

In figures 1–9, we demonstrate the evolution of the phase diagrams within the wide range of parameters β = t₂/t₁, λ = U_s/U₁₁, ξ₁₂ = U₁₂/U₁₁ and ξ₂₂ = U₂₂/U₁₁. The most interesting phase transition is that accompanied by the change of atom distribution over the optical lattice: p-spin FM ↔ p-spin AFM. Note that the cold colors correspond to FM p-spin, while the warm ones—to p-spin AFM. So the transitions with p-spin change can be found in the phase diagram at the lines, where cold colors change to warm ones.

Evolution of the phase diagrams for ξ₁₂ = U₁₂/U₁₁ = 0 with ξ₂₂ = U₂₂/U₁₁ on (λ, β)-plane is shown in figure 1. In this case, there is no Coulomb interaction between different boson species. We see that there is always quantum phase transition at the line λ = 0. This is true not only for figure 1, but also for all phase diagrams in figures 1–9. This phase transition is driven by the sign change of spin channel interaction U_s. This transition has been recently revealed in [31] for nearly identical vector bosons. Here, we show that this transition is quite robust with respect to the evolution of the degree of nonidentity.

One can also see in figure 1, that the 'left color' is always blue. It corresponds to FM spin ordering (with FM p-spin ordering). This situation is intuitively obvious: since the 'left' phase is determined by large negative λ—the interaction in the spin channel. In all the next figures, the 'left color' also corresponds to the FM spin ordering, sometimes with the AFM p-spin ordering (yellow).

Looking through the complete set of phase diagrams, one can notice the absence of exact symmetry with the respect to reflections λ → −λ and β → −β, though some traces are detectable. The β-symmetry is restored, when ξ₁₂ ≫ ξ₂₂.

Let us mention that the intuitive speculations useful, for example, in the case of simple Heisenberg model, can be misleading here, because the effective Hamiltonian is rather nontrivial.

We also underline the evolution of the artistic image of the phase diagrams. Namely, at small ξ₁₂ their style resembles the J Miró paintings, while at large ξ₁₂—those of K Malewicz.

5.2. Phase diagrams: specific features

To conclude, we have investigated the evolution of the quantum state of vector two species bosons in optical lattices with the 'degree of atomic nonidentity' that drives the cascade of quantum phase transitions. Surprisingly, this quite simple model exhibits rather rich and complicated manifold of quantum phase transitions between six different phases when interaction parameters of the model change in rather limited subspace. In fact, one can say that the system is 'unstable' with the respect to small variation of parameters.

Technically, we have reduced the initial general Hamiltonian for vector bosons, using the Mott insulating state natural small parameters, to the anisotropic spin-pseudospin model of the Kugel–Khomskii type that served as the effective Hamiltonian. This procedure gave a chance for analytical progress in investigations of the system ground state. The variational approach have been used to uncover the phase diagram of the system in hand. Finally, we have investigated also limiting cases of the effective Hamiltonian and demonstrated the relation of our rather complicated Hamiltonian to the well-known results. There are further perspectives to find new results in the two species fermion systems in optical lattices and hybrid boson-fermion ones looking at the 'degree of atomic nonidentity' as the driving parameter.

AM and AB are grateful to Wu-Ming Liu for interest to this work and hospitality in Beijing National Laboratory for Condensed Matter Physics. NC is grateful to A Pekovic for stimulating discussions at the initial stage of this work, to Laboratoire de Physique Théorique, Toulouse, where this work has been initiated, for the hospitality, and to CNRS.

This work was supported by the Russian Foundation for Basic Research (projects Nos. 16-02-00295, 16-02-00304, 17-52-53014, and 17-02-00323). We also express our gratitude to Government of the Russian Federation (project No. 05.Y09.21.0018) for support of the preparation of this paper revised version and the Computational Centers of Russian Academy of Sciences and National Research Center 'Kurchatov Institute' for providing the access to URAL, JSCC, and HPC facilities.

We remind that ${H}_{\mathrm{eff}}=-{\sum }_{\langle i,j\rangle }h$.

A.1. Equal Hubbard interaction parameters

Now, we take a look at more special cases. First we consider the case with equal Hubbard interaction parameters

$$\begin{eqnarray}&&{U}_{12}={U}_{11}={U}_{22}={U}_{0}.\end{eqnarray} \tag{ A.1 }$$

Then

$$\begin{eqnarray}&&{E}_{S}^{{ab}}={E}_{S}^{{aa}}={E}_{S}^{{bb}}={E}_{S},\qquad S=0,1,2,\end{eqnarray} \tag{ A.2 }$$

and as follows from equation (22), the effective matrix h reduces to the following form:

$$\begin{eqnarray}h & = & \displaystyle \frac{{Q}_{0}}{{E}_{0}}\left\{\displaystyle \frac{3}{2}({t}_{1}^{2}+{t}_{2}^{2})+2({t}_{1}^{2}-{t}_{2}^{2})({T}_{i}^{z}+{T}_{j}^{z})+2{({t}_{1}-{t}_{2})}^{2}\,{T}_{i}^{z}{T}_{j}^{z}+4{t}_{1}{t}_{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}\right\}\\ & & +\displaystyle \frac{{Q}_{2}}{{E}_{2}}\left\{\displaystyle \frac{3}{2}({t}_{1}^{2}+{t}_{2}^{2})+2({t}_{1}^{2}-{t}_{2}^{2})({T}_{i}^{z}+{T}_{j}^{z})+2{({t}_{1}-{t}_{2})}^{2}\,{T}_{i}^{z}{T}_{j}^{z}+4{t}_{1}{t}_{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}\right\}\\ & & +\displaystyle \frac{{Q}_{1}}{{E}_{1}}\left\{\displaystyle \frac{1}{2}({t}_{1}^{2}+{t}_{2}^{2})-2{({t}_{1}+{t}_{2})}^{2}\,{T}_{i}^{z}{T}_{j}^{z}+4{t}_{1}{t}_{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}\right\}.\end{eqnarray} \tag{ A.3 }$$

A.2. Equal Hubbard interaction parameters and equal hopping amplitudes

The effective matrix h can be simplified further if one considers equal hopping amplitudes, t₁ = t₂ = t. For this case, as it can be seen from equation (A.3), matrix h reduces to

$$\begin{eqnarray}&&h=\displaystyle \frac{{Q}_{0}}{{E}_{0}}\{3{t}^{2}+4{t}^{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}\}+\displaystyle \frac{{Q}_{2}}{{E}_{2}}\{3{t}^{2}+4{t}^{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}\}+\displaystyle \frac{{Q}_{1}}{{E}_{1}}\{{t}^{2}-8{t}^{2}\,{T}_{i}^{z}{T}_{j}^{z}+4{t}^{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}\}.\end{eqnarray} \tag{ A.4 }$$

The last expression can be simplified if we take into account that

$$\begin{eqnarray}&&{t}^{2}-8{t}^{2}\,{T}_{i}^{z}{T}_{j}^{z}+4{t}^{2}\,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}=4{t}^{2}\,{P}_{t0},\end{eqnarray} \tag{ A.5 }$$

where ${P}_{t0}={P}_{t}-({P}^{11}+{P}^{44})=[1/4-2{T}_{i}^{z}{T}_{j}^{z}+{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}]$ is the projector onto the pseudospin state $| T=1,{M}_{T}=0\rangle$, when the matrix h can be reduced to the form

$$\begin{eqnarray}&&h=4{t}^{2}\,\left\{\displaystyle \frac{{Q}_{0}}{{E}_{0}}{P}_{t}+\displaystyle \frac{{Q}_{2}}{{E}_{2}}{P}_{t}+\displaystyle \frac{{Q}_{1}}{{E}_{1}}{P}_{t0}\right\}.\end{eqnarray} \tag{ A.6 }$$

Note that due to the presence of the projector P_t0, the states with spin S = 1 will be automatically symmetric in the orbital space, while the antisymmetric combination of orbital states is automatically excluded from the effective Hamiltonian.

Next, we rewrite the Hamiltonian h in terms of the spin operators. We use the relation

$$\begin{eqnarray}&&-4{t}^{2}\,\left(\displaystyle \frac{{Q}_{0}}{{E}_{0}}+\displaystyle \frac{{Q}_{2}}{{E}_{2}}\right)=\left\{\displaystyle \frac{-4{t}^{2}}{3}\left(\displaystyle \frac{1}{{E}_{2}}-\displaystyle \frac{1}{{E}_{0}}\right)+\displaystyle \frac{-4{t}^{2}}{2{E}_{2}}\,{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+\displaystyle \frac{-4{t}^{2}}{3}\left(\displaystyle \frac{1}{{E}_{0}}+\displaystyle \frac{1}{2{E}_{2}}\right)\,{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2}\right\}\\ &&\quad =\,\epsilon +J\,{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+K\,{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2},\end{eqnarray} \tag{ A.7 }$$

where

$$\begin{eqnarray}&&J=\displaystyle \frac{-2{t}^{2}}{{E}_{2}},\quad K=\displaystyle \frac{-4{t}^{2}}{3}\left(\displaystyle \frac{1}{{E}_{0}}+\displaystyle \frac{1}{2{E}_{2}}\right),\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{-4{t}^{2}}{3}\left(\displaystyle \frac{1}{{E}_{2}}-\displaystyle \frac{1}{{E}_{0}}\right)=J-K.\end{eqnarray} \tag{ A.9 }$$

Parameters J, K, and can be expressed explicitly via the initial interaction constants U₀ and U_s as follows

$$\begin{eqnarray}&&J=\displaystyle \frac{-2{t}^{2}}{{U}_{0}+{U}_{s}},\quad K=\displaystyle \frac{-2{t}^{2}{U}_{0}}{({U}_{0}+{U}_{s})({U}_{0}-2{U}_{s})},\end{eqnarray} \tag{ A.10 }$$

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{4{t}^{2}{U}_{s}}{({U}_{0}+{U}_{s})({U}_{0}-2{U}_{s})}.\end{eqnarray} \tag{ A.11 }$$

Now, the Hamiltonian takes the form

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}=\sum _{\langle i,j\rangle }\left\{[\epsilon +J\,{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+K\,{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2}]{P}_{t}+\left(\displaystyle \frac{-4{t}^{2}}{{E}_{1}}\right){Q}_{1}{P}_{t0}\right\}.\end{eqnarray} \tag{ A.12 }$$

Finally, the relation ${P}_{t}={P}_{t0}+{P}^{11}+{P}^{44}$ and the definition of Q₁ (16) allows us to reduce the Hamiltonian to the form used in [31],

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}=\sum _{\langle i,j\rangle }\{[\epsilon +J\,{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+K\,{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2}]({P}^{11}+{P}^{44})+[\epsilon ^{\prime} +J^{\prime} \,{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+K^{\prime} \,{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2}]\,{P}_{t0}\},\end{eqnarray} \tag{ A.13 }$$

where we introduced

$$\begin{eqnarray}&&J^{\prime} =2{t}^{2}\left(\displaystyle \frac{1}{{E}_{1}}-\displaystyle \frac{1}{{E}_{2}}\right),\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}&&K^{\prime} =\displaystyle \frac{-2{t}^{2}}{3}\left(\displaystyle \frac{2}{{E}_{0}}+\displaystyle \frac{1}{{E}_{2}}-\displaystyle \frac{3}{{E}_{1}}\right),\end{eqnarray} \tag{ A.15 }$$

$$\begin{eqnarray}&&\epsilon ^{\prime} =\displaystyle \frac{-4{t}^{2}}{3}\left(\displaystyle \frac{1}{{E}_{2}}-\displaystyle \frac{1}{{E}_{0}}+\displaystyle \frac{3}{{E}_{1}}\right).\end{eqnarray} \tag{ A.16 }$$

A.2.1. One type of bosons

For the case of only one type of bosons, Hamiltonian (A.13) is equivalent to that of [43]. For single type of bosons, the pseudospin operators become c-numbers: ${T}_{i}^{z}=\pm 1/2,{{\bf{T}}}_{i}\cdot {{\bf{T}}}_{j}$ should be replaced by ${T}_{i}^{z}{T}_{j}^{z}=1/4$, and projectors ${P}^{11}+{P}^{44}=1,{P}_{t0}=0$. Then equation (A.13) is reduced to the Hamiltonian considered in [43]

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}=\sum _{\langle i,j\rangle }\{\epsilon +J\,{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+K\,{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2}\}.\end{eqnarray} \tag{ A.17 }$$

A.2.2. No spin-dependent interaction

If there is no spin-dependent interaction, i.e. U_s = 0, then E₀ = E₂ = E₁ and J = K = −2t²/U₀,  = 0, J' = K' = 0, and $\epsilon ^{\prime} =-4{t}^{2}/{U}_{0}=2J$. The effective Hamiltonian, equation (A.13), is reduced to

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}=J\sum _{\langle i,j\rangle }\{2{P}_{t0}+({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2})({P}^{11}+{P}^{44})\}.\end{eqnarray} \tag{ A.18 }$$

And at last, if in this case, there is only one type of bosons, then ${P}_{t0}=0,{P}^{11}+{P}^{44}=1$, and the effective Hamiltonian describing effective interaction between identical bosons has the FM character (J < 0)

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}=J\sum _{\langle i,j\rangle }\{{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+{({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j})}^{2}\},\quad J=\displaystyle \frac{-2{t}^{2}}{{U}_{0}}.\end{eqnarray} \tag{ A.19 }$$

The systems with identical bosons with odd and even number bosons per site were discussed in [40, 41].