The effect of disorder on the phase diagrams of hard-core lattice bosons with cavity-mediated long-range and nearest-neighbor interactions

We use quantum Monte Carlo simulations with the worm algorithm to study the phase diagram of a two-dimensional Bose-Hubbard model with cavity-mediated long-range interactions and uncorrelated disorder in the hard-core limit. Our study shows the system is in a supersolid phase at weak disorder and a disordered solid phase at stronger disorder. Due to long-range interactions, a large region of metastable states exists in both clean and disordered systems. By comparing the phase diagrams for both clean and disordered systems, we find that disorder suppresses metastable states and superfluidity. We compare these results with the phase diagram of the extended Bose-Hubbard model with nearest-neighbor interactions. Here, the supersolid phase does not exist even at weak disorder. We identify two kinds of glassy phases: a Bose glass phase and a disordered solid phase. The glassy phases intervene between the density-wave and superfluid phases as the Griffiths phase of the Bose-Hubbard model. The disordered solid phase intervenes between the density-wave and Bose glass phases since both have a finite structure factor.


I. INTRODUCTION
The study of adding disorder to the interacting manybody bosonic systems attracts enormous attention both experimentally and theoretically [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. Experimentally, ultracold atoms in optical lattices are a promising tool to study quantum phases and quantum phase transitions in strongly correlated quantum many-body systems. It provides an unique possibility of engineering matter with an unprecedented level of control over the parameters entering the Hamiltonian. On one hand, short-range interactions can be realized using Feshbach resonances, while long-range interactions have been studied using ultracold gases of particles with large magnetic or electronic dipole moments [15][16][17], polar molecules [18,19], atoms in Rydberg states [20][21][22], or cavity-mediated interactions [23,24]. On the other hand, disordered potential can be introduced artificially into the ultracold atomic gases in optical lattices. Speckle patterns is the most commonly used to produce the random potential [2,25,26]. Bichromatic lattices [27], introduction of localized atomic impurities [4], and holographic techniques which produce point-like disorder [28] are also used to engineer disorder experimentally.
The effect of disorder on the phases and phase transitions of quantum many-body systems triggered many theoretical studies [29][30][31]. The disordered Bose-Hubbard model (DBHM) allows one to study the interplay between disorder and interactions of ultracold bosons in optical lattices. In the phase diagram of the DBHM, the gapless Bose glass (BG) phase, characterized by a finite compressibility and absence of an offdiagonal long-range order, always intervenes as a Griffiths phase between the superfluid (SF) and Mott insulator (MI) phases [32,33]. The original DBHM focuses on short-range on-site interactions [34]. Recently long-range interactions starts attracting the focus of theoretical research. In the absence of disorder, with longrange interactions, the BHM exhibits a richer phase di-agram with additional density wave (DW) and supersolid (SS) phases [35][36][37][38]. The ground state phase diagram of the extended BHM with cavity-mediated longrange interactions has been investigated extensively with the help of mean-field theory [37,[39][40][41][42], Gutzwiller ansatz [38,43], quantum Monte Carlo [36,38,41], and Variational Monte-Carlo [44] methods in 1D, 2D, and 3D. The addition of disorder to the BHM with long-range interactions leads to additional phases. In our recent study, we found that in the DBMH, long-range interactions enhance the supersolid phase [45]. However, a study of the DBHM with cavity-mediated long-range interactions in the hard-core limit is still lacking. In the hard-core limit and without disorder, equilibrium phases of lattice bosons with cavity-mediated long-range interactions were investigated in [46,47] in 1D where the result shows that the checkerboard supersolid does not exist. While in 2D with nearest-neighbor interactions [48], the result shows that the checkerboard supersolid is unstable. However, in the presence of disorder and with hard-core limit, whether the supersolid phase exists or not is still unknown.
In this paper, we use quantum Monte Carlo simulations with the worm algorithm [49] to study the phase diagram of the two-dimensional disordered Bose-Hubbard model with cavity-mediated long-range and nearest-neighbor interactions in the hard-core limit. The paper is organized as follows: in section II, we introduce the Hamiltonian of the system of hard-core bosons with cavity-mediated long-range and nearest-neighbor interactions in the presence of disorder. In section III A, we present the phase diagrams of hard-core bosons in the two-dimensional lattice with cavity-mediated longrange interactions for both clean and disordered systems. We also study the phase diagram of hard-core bosons in the two-dimensional lattice with nearest-neighbor interactions for both clean and disordered systems in section III B. Section IV concludes this paper. We consider bosons trapped in an optical lattice with both short-range on-site and cavity-mediated long-range interactions in the presence of disorder in the hard-core limit. The bosons are trapped in a two-dimensional (2D) square lattice with linear size L and periodic boundary conditions. The hard-core limit corresponds to large onsite interactions where the occupation of two bosons on the same lattice site is suppressed. The system is described by the Hamiltonian [35][36][37]: Here the first term is the kinetic energy characterized by the hopping amplitude t. · · · denotes nearest neighboring sites, a † i (a i ) are the bosonic creation (annihilation) operators satisfying the usual bosonic commutation relations and n i = a † i a i is the particle number operator. The second term is the cavity-mediated long-range interaction with interaction strength U l , the summations i ∈ e and j ∈ o denote summing over even and odd lattice sites, respectively [37]. The third term is the chemical potential term with chemical potential µ shifted by the on-site random potential ε i , where ε i is uniformly distributed within the range [−∆, ∆]. ∆ is the disorder strength. The hard-core condition a †2 i = 0 implies that sites with more than one atom are energetically suppressed due to a large on-site interaction energy penalty. In this limit the usual on-site interaction term of the Bose-Hubbard model does not play any role. The maximum atom per site is one. The unit of energy and length are set to be the hopping amplitude t. For each µ/U l , t/U l , and ∆/t, we average over 200-400 realizations of disorder.
We also consider the extended Bose-Hubbard model with nearest-neighbor interactions in the presence of disorder in the hard-core limit. The Hamiltonian is written as: Here, the first term is the kinetic energy with hopping amplitude t. The second term is the repulsive interaction with interaction strength U nn between bosons on nearest neighboring sites. The third term is the random potential term coupled with the chemical potential term. For each µ/U nn , t/U nn , and ∆/t, we average over 1000-3000 realizations of disorder.

III. GROUND STATE PHASE DIAGRAMS
In this section, we present the ground state phase diagram of the extended BHM with cavity-mediated long-range interactions U l (model 1) in section III A and nearest-neighbor interactions U nn (model 2) in section III B in the hard-core limit for both clean and disordered systems, respectively. To obtain the phase diagram for cavity-mediated long-range interactions ( Fig. 1) and nearest-neighbor interactions (Fig. 4), we measure the superfluid stiffness, compressibility, and structure factor to separate those quantum phases. Different phases can be distinguished by the different combinations of those order parameters. Table I shows quantum phases and corresponding parameters associated to the phase we find in phase diagrams Fig. 1 and Fig. 4. TABLE I. Quantum phases and the corresponding parameters: superfluid stiffness ρ, structure factor S(π, π), and compressibility κ.
The superfluid (SF) phase is characterized by a finite superfluid stiffness, which is easily accessible in the QMC simulations by calculating the winding number Here, W is the winding number. d = 2 is the dimension of the system, L is the linear size of the system, and β is the inverse temperature. The density wave (DW) phase has a finite structure factor which is defined as S(k) = r,r exp [ik(r − r )] n r n r /N . Here, k is the reciprocal lattice vector and k = (π, π). The supersolid (SS) phase possesses both the diagonal long-range order and off-diagonal long-range order, and it is characterized by finite ρ and S(π, π). Both the disordered solid (DS) phase and the Bose glass (BG) phase are characterized by a finite compressibility, the difference between them is that DS phase has a finite structure factor. The compressibility measures the density fluctuations and it is defined as: κ = β( n 2 − n 2 ). The Mott insulator (MI) phase has zero superfluid stiffness, zero compressibility, and zero structure factor.
A. Cavity-mediated long-range interaction Figure 1 shows the ground state phase diagram of hard-core bosons trapped in a two-dimensional optical lattice with cavity-mediated long-range interactions for both clean ( Fig. 1 (a)) and disordered systems ( Fig. 1 (b)-(d)), respectively. The x-axis is t/U l and the y-axis is µ/U l where t is the hopping amplitude, U l is the strength of cavity-mediated long-range interactions, and µ is the chemical potential. Here we set the hopping amplitude t = 1. The phase boundary is determined by considering cuts through the x axis (t/U l ) and calculating the above three order parameters as a function of µ/U l . Finite-size scaling method is also used to get accurate transition points on the phase boundary. Figure 1 (a) shows the ground state phase diagram of the clean system. There are three phases in the phase diagram: a SF phase, a MI phase, and a DW phase. At small t/U l , there is a large region of metastable states (MS). A metastable state is a state that the system resides in but differs from its state of least energy for an extended period of time. Here, cavity-mediated long-range interactions tend to induce metastability. Experimentally, the MS can be probed by preparing the system in a well-defined state, providing extra energy, and observing which state the system relaxes to [51]. In the simulation, the MS is determined by starting simulations with two well-defined initial states (DW and MI) and measuring order parameters to determine the final state. If the final state depends on the initial state, the system is in the MS. Figure 1 (a) shows that the MS surrounds the DW phase for clean system at a small value of t/U l . The MS and DW persists until t/U l = 0.32 ± 0.02. When t/U l > 0.32, the system stays in the SF phase at a lower filling and MI at an integer filling n = 1. Here, the MI to SF phase transition is the second-order phase transition. Figure 1 (b) shows the ground state phase diagram at disorder strength ∆/t = 2. As disorder is added to the system, the region of MS shrinks. Disorder tends to localize bosons and suppress metastability. Besides the SF, MI and DW phases, there are two new phases emerge: a BG phase and a supersolid (SS) phase. The BG phase intervenes as a Griffiths phase between the MI and SF phases, and it is explained by the theory of inclusions [32,33]. The SS phase has both diagonal long-range order and off-diagonal long-range order and is characterized by a finite superfluid stiffness ρ and a finite structure S(π, π). Figure 2 (a) shows the structure factor S(π, π) and superfluid stiffness ρ as a function of µ/U l for different system sizes L = 8, 10, 12, and 14 (red dots, blue rectangles, green up triangles, and purple diamonds) at disorder strength ∆/t = 2 and t/U l = 0.2857. Between 0 < µ/U l < 0.25, the system is in the SS phase with both a finite superfluid stiffness and a finite structure factor. Figure 2 (b) shows the finite-size scaling of structure factor S(π, π), where we plot S(π, π)L 2β/ν as a function of µ/U l for system sizes L = 8, 10, 12, and 14 (the critical exponents 2β/ν = 1.0366 (8) correspond to the three-dimensional Ising universality class [52]). The crossing of different curves marks the transition point at µ/U l = 0.2 ± 0.05. The insert shows the data collapse result using ν = 0.67 and∆ c = (µ/U l ) c = 0.2 corresponding to the critical point extracted from main plot. Here, the SS goes to the SF phase via a second-order phase transition which belongs to the (2+1)-dimensional Ising type transition. Figure 3 shows the finite-size scaling of superfluid stiffness ρ and structure factor S(π, π) as a function of µ/U l for different system sizes at ∆/t = 2 and t/U l = 0.25. Dashed lines are the finite-size scaling result of superfluid stiffness. We plot ρL d+z−2 as a function of µ/U l at t/U l = 0.25 for a variety of system sizes. Here, z = 1 is the dynamic critical exponent and the inverse temperature β = L is used. The crossing of different curves marks the transition point at µ/U l = 0.325 ± 0.01. This shows at t/U l = 0.25, the DW to SS phase transition is the second-order phase transition and happens at µ/U l = 0.325 ± 0.01. Solid lines are the finite-size scaling result of structure factor. The crossing shows that the SS to SF phase transition happens at µ/U l = 0.368 ± 0.02. At t/U l = 0.25, the system is in the SS phase at 0.325 < µ/U l < 0.368. Figure 1 (c) shows the ground state phase diagram at disorder strength ∆/t = 5. As disorder increases, both the MS and SF phases shrink. This is due to the fact that disorder tends to localize bosons and destroy superfluidity and metastability. Interestingly, we do not find the SS phase but the DS phase. The DS is characterized by a finite compressibility and a finite structure factor but no superfluid stiffness. Figure 1 (d) shows the ground state phase diagram at disorder strength ∆/t = 10. At such a strong disorder, there is no SF phase anymore. By comparing phase diagrams in Figure 4, we can see that disorder tends to shrink the region of MS states and destroy superfluidity. The BG phase intervenes as a Griffiths phase between the MI and SF phases. The SS or DS intervenes between the DW and SF phases depending on the strength of disorder.

B. Nearest-neighbor interaction
In this subsection, we study the ground state phase diagram of hard-core bosons trapped in an optical lattice with nearest-neighbor interactions for both clean (Fig. 4  (a)) and disordered systems (Fig 4 (b)-(d)), respectively. The x-axis is t/U nn and the y-axis is µ/U nn , where t is the hopping amplitude, U nn is the strength of nearestneighbor interactions, and µ is the chemical potential. The hopping amplitude is set to t = 1. The phase boundary is determined by using system size L = 16, finite size scaling method is also used to get accurate transition points on the phase boundary. Figure 4 (a) shows the ground state phase diagram of the extended BHM with nearest-neighbor interactions for FIG. 1. Ground state phase diagrams of model 1 as a function of t/U l and µ/U l for clean system (a), and disordered system with disorder strength ∆/t = 2 (b), 5 (c), and 10 (d), respectively. The gray shadowed region represents the metastable states.
the clean system. There are three phases in the phase diagram: a SF phase, a MI phase, and a DW phase. Here, the MI to SF phase transition is the second-order phase transition at an integer filling while the DW to SF phase transition is the first-order phase transition at a half filling [48]. Figure 4 (b) shows the ground state phase diagram at disorder strength ∆/t = 3. As disorder is added to the system, the glassy phase (BG and DS) emerge. This can be explained by the theory of inclusions [32,33], which states that a compressible glassy phase surrounds the incompressible phase. There are two kinds of glassy phase, the BG phase and DS phase. The BG phase is characterized by a finite compressibility κ but zero structure factor S(π, π) while the DS phase is characterized by both a finite compressibility and a finite structure factor. The incompressible phase here is the DW phase. As t/U nn increases, the system goes from the DW to DS phase transition and then the DS to BG phase transition. The DS phase intervenes between the DW and BG phase since both have a finite structure factor. Finitesize scaling method is used to determine all critical points on the phase boundary. As shown in Figure 5, at fixed t/U nn = 0.2, as µ/U nn increases, the DW phase goes to the DS phase first. The main plot of Figure 5 (a) shows the finite-size scaling of S(π, π) at fixed disorder strength ∆/t = 3 and t/U nn = 0.2 for system sizes L = 12, 16, 20, and 24 (red circles, blue rectangles, green up triangles, and orange diamonds, respectively). The crossing of different curves marks the transition point at µ/U nn = 3.275 ± 0.02. Insert shows the data collapse result using ν = 0.67 and∆ c = (µ/U nn ) c = 3.275 corresponding to the critical point extracted from main plot, which shows the DS to BG phase transition belongs to the (2+1)-dimensional Ising type transition. Figure 5 (b) shows the scaling of superfluid stiffness ρL d+z−2 with z = 1, as a function of µ/U nn for t/U nn = 0.2 and L = 12, 16, 20, and 24. Here, z is the dynamic critical exponent and the inverse temperature β = L is used. The crossing of different curves marks the transition point at µ/U nn = 3.652 ± 0.02, which shows that at t/U nn = 0.2, the BG to SF phase transition is the second-order phase transition and happens at µ/U nn = 3.652 ± 0.02. increases, the DW phase is unstable and the system goes to the DS phase. Figure 6 (a) shows the compressibility κ as a function of t/U nn for different system sizes L = 16, 20, 24, and 30 (blue rectangles, green up triangles, orange diamonds, and purple down triangles, respectively). The compressibility κ stays zero until t/U nn = 0.18 and then becomes finite, which shows the DW to DS phase transition happens around t/U nn = 0.18 ± 0.01. As t/U nn increases further, the DS is unstable and the system enters the BG phase. Figure 6 (b) shows the finite-size scaling of S(π, π) for above system sizes. The crossing of different curves marks the transition point at t/U nn = 0.212±0.01. Insert shows the data collapse result using ν = 0.67 and∆ c = (t/U nn ) c = 0.212 corresponding to the critical point extracted from main plot. This shows the DS goes to the BG phase via a second-order phase transition which belongs to the (2+1)-dimensional Ising type transition. Finally, the BG goes to the SF as t/U nn increases further. Figure 6 (c) shows the finite-size scaling of superfluid stiffness ρL as a function of t/U nn for above system sizes. The crossing of different curves marks the transition point at t/U nn = 0.3241 ± 0.005. Figure 4 (d) show the ground state phase diagram at disorder strength ∆/t = 10. At such a strong disorder, there is no SF phase anymore. By comparing phase diagrams in Fig. 4, we can see that disorder tends to localize bosons and destroy superfluidity. Compared to the disordered BHM without nearest-neighbor interactions [8], we find that the DW phase can not go to the BG phase directly, there is the DS phase intervenes between them. the scaling of superfluid stiffness ρL and structure factor S(π, π)L 1.0366 as a function of µ/U l for different system sizes L = 8, 10, 12, and 14 (red dots, blue rectangles, green up triangles, and purple diamonds) at ∆/t = 2 and t/U l = 0.25. Solid lines represent structure factor and dashed lines represent superfluid stiffness. The crossing of different solid curves marks the transition point at µ/U l = 0.325 ± 0.01 for the density-wave to supersolid phase transition while the crossing of different dashed curves marks the transition point at µ/U l = 0.368 ± 0.02 for the supersolid to superfluid phase transition.

IV. CONCLUSION
In this paper, we use quantum Monte Carlo simulations with the worm algorithm to study the phase diagram of a two-dimensional Bose-Hubbard model with cavity-mediated long-range interactions for both clean and disordered systems in the hard-core limit. We find the SS phase at a weak disorder and DS phase at a stronger disorder. Due to the long-range interaction term, we find a large region of metastable states in both clean and disordered systems. Disorder suppresses metastable states and superfluidity. We also study the phase diagram of the extended Bose-Hubbard model with nearest-neighbor interactions for both clean and disordered systems in the hard-core limit. We find that there are two kinds of glassy phases: the BG phase and the DS phase. The glassy phase intervenes as a Griffiths phase between the DW and SF phases, which can be explained by the theory of inclusions. The DS phase intervenes between the DW and BG phase since both the DS and DW phases have a finite structure factor.