Disorder, three body interaction and Bose glass phase in a spinor atomic gas in an optical lattice

We study the effects of disorder on the spin dependent interaction term of a spinor Bose Hubbard model with a three body interaction potential. The signature of the Bose glass (BG) phase is observed by computing the fraction of the lattice sites having finite superfluid (SF) order parameter and non integer occupation densities. We obtain the phase diagram both for the antiferromagnetic (AF) and ferromagnetic (F) cases via a percolation analysis. In the AF case, the BG phase intervenes between the odd-even Mott insulating (MI) lobes (for example, the lobes corresponding to occupation densities, n=1 and (n=2) but not between the even-odd MI lobes. In the ferromagnetic case, the presence of the BG phase is observed between all the MI lobes irrespective of them being even or odd. The BG phase almost destroys the first MI lobe while the MI phase looks more stable than the SF phase both in the AF and F cases due to the presence of the three body interactions.


Introduction
Ultracold atoms trapped in optical lattices offer a novel platform to explore the quantum phenomena where quantum fluctuations are responsible for phase transition from one macroscopic phase to another one. In magnetic trapping, ultracold atoms freeze their internal degrees of freedom and thereby behaves like scalar or spin-0 Bose gas while in optical trapping, they retains their spin degrees of freedom and show spinor like behaviour [1].
The properties of cold atoms in optical lattices first theoretically described by the well known Bose Hubbard model (BHM) developed by Jaksch et. al. [2] where the phase transition is entirely determined between the competition in tunneling and two body interaction strengths. Following the proposal by Jaksch et. al., Greiner et. al. [3] first experimentally realized the superfluid (SF) to Mott insulator (MI) transition of bosonic cold atoms loaded in optical lattices.
Several progress have been made to study the scalar or spinor Bose gas with different extensions of BHM in both homogeneous and inhomogeneous case from theoretical as well as experimental point view. Among them, one possible extension of the BHM include the study of three body interactions [4][5][6][7] on the SF to MI transition gained much attention after the successful experimental realization of three and higher body interaction strengths by Will et. al. [8].
In this work, we will study the effect of three body repulsive interaction strengths on spinor ultracold Bose gas in presence of disorder through mean field approach (MFA) [9] where an additional glassy type of phase known as Bose glass (BG) phase [10][11][12]14] intervene in between SF and MI due to the presence of disorder. Since MFA is quite inadequate to overcome the site dependencies in the Hamiltonian, numerical approach like quantum Monte Carlo (QMC) [10], stochastic mean field theory (SMFT) [11], probabilistic mean field perturbative (PMFA) approach [12,13] and strong coupling expansion [14] have been developed to study the BHM in presence of disorder.
Following the approach developed in [15], here we use the MFA to study the disordered BHM and define a new quantity which is the fraction of lattice sites having finite SF order parameter to remove the inhomogenities in the Hamiltonian and did a percolation analysis to find the extent of the different phases based on the appearance of SF percolating cluster.

Spinor Bose Hubbard Hamiltonian
The Bose Hubbard Hamiltonian (BHM) for spin-1 ultracold atoms in presence of random disorder with three body interaction term is [4][5][6]16], where t is the tunneling amplitude and µ is the chemical potential.â † iσ (â iσ ) is the boson creation (annihilation) operator at a site i with hyperfine spin state σ =+1, 0, -1 and the particle number operator isn i = σâ † iσâ iσ . The hyperfine spin operator at site i is given by S i =â † iσ F σσ ′â iσ ′ where F σσ ′ are the components of spin-1 matrices. W represents the three body interaction strength and U 0 and U 2 are spin independent and dependent two body interaction strengths respectively. U 0 and U 2 are related to the scattering lengths, a 0 and a 2 as U 2 /U 0 = (a 2 − a 0 )/(a 0 + 2a 2 ) and is called antiferromagnetic (AF) interaction if U 2 /U 0 > 0 and ferromagnetic (F) if U 2 /U 0 ≤ 0. ǫ i is the on-site disorder is introduced at site i, which is randomly chosen from a box distribution in the interval [−∆, ∆] where ∆ is the strength of the disorder. To decouple the hopping part in Eq.(1), we shall use mean field approximation as [17] where denotes the equilibrium value of an operator. Now defining the superfluid order parameter at a site i as, ψ iσ ≡ â † iσ ≡ â iσ and ψ i = ψ 2 iσ , the BHM can be written as a sum of single site Hamiltonians as, where φ iσ = j ψ jσ , j includes all nearest neighbors at the site i of a square lattice with z = 4, z being the coordination number. The local density, ρ i and compressibility, κ i can also be computed using, where ψ eq i is the equilibrium value of ψ i after self consistent ground state is obtained.

Results
In order to see the effect of disorder and three body interaction term, first we consider the atomic limit i.e t = 0. In the atomic limit t = 0, system will consist only the Mott insulating (MI) phase where there is a gap in the particle hole excitation spectrum. The energy gap, E g for a particular filling of the MI lobe n corresponds to the upper and lower boundaries of the MI region as µ − (n) < µ < µ + (n) [12]. In the AF case, this inequality leads to the following equations corresponding to the even MI lobes as Fig.3(a). It shows that both in the AF and F cases, the chemical potential width µ increases due to W/U 0 which suggests that the MI lobes becomes more stable compared to the SF phase. Interestingly the critical value of U 2 for the disappearance of the odd MI lobes in the AF case now depends on the three body interaction strengths which was 0.5 without W/U 0 .
We now turn on the disorder and diagonalize the mean field Hamiltonian on a 2D square lattice of size L = 30 self consistently to find the ground state energy and the equilibrium order parameter ψ eq i andψ = 1   (1))]. In the atomic limit, t = 0 for AF case, the chemical potential width µ for the even-odd transition is parallel to U 2 axis while for odd-even transition, it depends upon U 2 and disorder in U 2 will change the chemical potential width and hence the MI lobes [ Fig.3 (a)]. In F case, the BG phase intervenes in all the MI lobes since disorder in U 2 will affect the occupation densities [Figs.2 (a),(b)]. From the above discussion, although we have been able to identify the BG phase but in order to get a clear understanding, we have to determine the boundaries of these three phases and hence the phase diagram in presence of disorder. In MFA, all SF parameter varies randomly due to disorder, thus it is quite difficult to determine the location of different phases and hence the phase diagrams based on the averageψ andκ.
To remove the inhomegenities in the Hamiltonian, we define a new quantity, called χ as, χ = sites with ψ i = 0 and ρ i = integer total numbers of sites The variation of χ for AF case is shown in Fig.3(b). It shows that χ switches from 0 to 1 indicating a direct MI to SF phase transition in clean state while in presence of disorder, there is gradual increase in slope due to the intervening BG phase in between MI and SF phases.
Interestingly in AF case, χ only assumes 0 and 1 values corresponding to even-odd MI transition in both clean and disordered cases leaving the signature of BG phase [ Fig.3(b)(MI(2))]. While for the odd-even transition χ shows a gradual increase between 0 and 1, indicating the intervention of BG phase in between MI and SF phases as discussed earlier [ Fig.3(b)(MI(1))]. In F, χ takes all finite values signaling the appearance of BG phase for all MI lobes and shows similar properties if disorder is present in µ.
From the variation of χ, it is now quite clear that χ = 0 corresponds to MI and χ = 1 corresponds to SF phase but we have to determine if there is any critical value of χ i.e χ c below which the system is in BG phase above which it belongs to SF phase for 0 < χ < 1 region [ Fig.3(b)]. To answer the above concern, we will follow the percolation analysis to know if there is any SF percolating cluster exists. We have characterized the three phases based on the appearances of SF percolating cluster [15]. A SF percolating cluster means cluster of sites with non integer densities connected with their nearest neighbors that percolates throughout  the lattice, that is from left to right and/or top to bottom. In MI phase, all sites have integer densities thus no SF percolating cluster. While in SF phase, all sites have non integer densities and thus have at least one SF percolating cluster. But in the BG phase, some sites have non integer densities, while the rest will have integer densities. Hence the BG phase can be identified with having a SF cluster which does not percolate. Whether a SF cluster will percolate can be well understand from the following quantity, χ ∞ known as percolation probability, defined as the fraction of sites within the SF percolating cluster out of total no of occupied sites [18]. To find a SF percolating cluster, we shall use Hoshen-Kopelman(HK) algorithm to find a SF percolating cluster [19]. It is based on the special application of the union-find algorithm where it aims to assign a label to each cluster. For that we represent our 2D lattice as a matrix and label all occupied sites initially with -1 and unoccupied sites with 0. Our cluster index starts with 1. Corresponding to each occupied site, every time we check the neighbors at the top and left corner of the current site. If both sites are empty, we label a new cluster number that has not used so far. Else, if the site has one occupied neighbor, then we assign same cluster number to the current site. If both neighboring sites are occupied, then we set the smallest number cluster label of the occupied neighbors to use as the label for the current site. To link two clusters, we create a union between both labels and set the site as the lowest of the two labels. When we burn the lattice for a second time, we collect the unions and update the lattice [19]. χ ∞ is analogous to critical percolation probability p c which is 0.593 for infinite square lattice that is L = ∞ and the variation of χ ∞ with χ is shown in Fig.4(a). It shows that χ ∞ is 0 till χ < 0.524 means no SF percolating cluster exists while it is non zero above χ ≥ 0.524 where SF percolating cluster exists.
At this stage, we have now all information to determine the location of different phases and thereby the phase diagrams entirely based on χ and χ ∞ . In clean state, phase diagram for AF and F cases are obtained by either ψ or χ as shown in Figs.4(b) and 5(b). Although there is no change in the first MI lobe but higher order MI lobes becomes more stable and the transition from MI to SF phase increases due to W/U 0 . In the disordered case, we have obtained the phase diagrams by incorporating the χ = 0 for MI phase, for BG phase, 0 < χ < χ c = 0.43 corresponding to AF and 0 < χ < χ c = 0.53 corresponding to F cases below which χ ∞ = 0 and χ c ≤ χ ≤ 1 for the SF phase and shown in Figs.4 and 5. It is quite helpful to mention that the colors code in the legend are as follows-blue: MI, red: BG and yellow: SF phases respectively. In AF case, BG phase only appears in between odd-even MI phase and due to inclusion of three body interaction, the intervening BG phase becomes quite thin compared to the MI phase [ Fig.4(c)]. While in F case, the BG appears in between all MI lobes shows similar phase diagrams if disorder is present in µ [ Fig.5(c)]. We have also considered the disorder in µ  Fig.5(a), the BG can not destroy the MI lobes since the disorder strength needed to dismiss the respective MI lobes is now lower than the critical value ∆ c which depends on W/U 0 [12]. Similarly in F case, the first MI lobe vanishes due the BG phase as ∆ is higher than the critical disorder strength [ Fig.5(c)] [12].

Conclusion
We study the effect of spin dependent disorder in spinor BHM through percolation analysis in presence of three body interactions. The different phases in presence of disorder are observed by the SF order parameter and compressibility. The boundary for the MI, BG and SF phases are determined with the help of percolation analysis and finally obtain the phase diagrams based on χ. In AF case, the BG appears only for the odd-even MI transition and very thin compared with the MI lobes due to the presence of three body interaction. We have also studied the effect of disorder in ferromagnetic case and found that it has similar properties as that of disorder in the chemical potential µ.