{'fraction_non_alphanumeric': 0.03349009548239989, 'fraction_numerical': 0.014322359982898675, 'mean_word_length': 4.213595839524517, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Monte-Carlo study of a two-species boson Hubbard model

We consider a two-species hard-core boson Hubbard model for a supersolid, where the two types of bosons represent vacancies and interstitials doped into a commensurate crystal. The on-site inter-species interaction may create bound states of vacancies and interstitials facilitating vacancy condensation at lower energies than in a single-species model, as suggested in an earlier mean field study. Here we carry out quantum Monte Carlo simulation to study possible supersolid phases of the model, corresponding to superfluid phases of the vacancies or interstitials. At low temperatures, we find three distinct superfluid phases. The extent of the phases and the nature of the phase transitions are discussed in comparison to mean-field theory.


I. INTRODUCTION
A supersolid is a special type of solid with superfluid properties. It has a diagonal particle density long range order as in a usual crystal, and an off-diagonal long range order in particle density as in a superfluid. The simplest model for supersolid was proposed by Andreev and Lifshitz in 1969. 1 Their model was introduced to describe possible supersolid phase in Helium-4. In their model, vacancies or interstitials of solid Helium may exist in the ground state and condense due to the large quantum fluctuation of Helium atoms. The interaction between vacancy and interstital is neglected in their model.
In this paper we study a two-species boson Hubbard model, which is an extension of the Andreev-Lifshitz model to include the interaction between vacancy and interstitial. This two-species model was recently introduced by Dai, Ma, and Zhang, 2 motivated by the observation of non-classical rotational inertia moment in solid helium-4 reported by Kim and Chan. 3 They used a mean field theory to study the ground state of the model and the possibility of the supersolid phase. It was shown that the interaction of vacancies and interstitials may facilitate a supersolid phase. In this paper we use quantum Monte Carlo (QMC) simulations to study the possible supersolid and the finite temperature phase transition in the two-species boson model. The simulations support the qualitative conclusion obtained in the mean field theory that the vacancy-interstitial interaction may facilitate supersolidity. Using QMC, we calculate the phase diagram, the superfluid densities of bosons, and the specific heat of the system. The two-species boson model and our calculations may be useful to understand other boson problems such as bosons in optical lattices. 4 Before we present the model and our results, we briefly summarize the current situation in study of supersolid Helium-4. Because of its light mass and its bosonic nature, solid helium-4 has been a natural candidate for possible supersolid at low temperatures and high pressures. Theoretically, such a possibility was proposed by Andreev and Lifshitz 1 and by Chester. 5 Leggett further pre-dicted the non-classical rotational inertia moment of such a supersolid in a rotating experiment. 6,7 The interest of supersolid has been recently revived due to the observation of non-classical rotational inertia in solid helium-4. 3 By now, the non-classical inertia moment in solid helium has been confirmed by other groups. 8,9,10 However, it remains controversial if the phenomenon is related to the supersolidity and if the supersolid phase is a bulk equilibrium phenomenon. 11,12,13 On the theoretical side QMC simulations did not find a supersolid phase in Helium-4. 14,15,16 Furthermore, the vacancies or interstitials in helium are shown to attract to each other and to tend to have phase-separation, 16 indicating that the Andreev-Lifshitz model may not describe solid helium.

II. MODEL AND METHOD
We consider a two-species boson Hubbard model in a cubic lattice with z = 6 nearest neighbors: where a j is an annihilation operator of boson a at lattice site j, representing a vacancy, and b j an annihilation operator of boson representing an interstitial, in a vacuum representing a defect-free insulating crystal of bosonic atoms. n j,a = a † j a j and n j,b = b † j b j are the number operators for aand b-bosons, and ǫ a and ǫ b are site boson energies, respectively. We consider the interesting case ǫ a > 0, and ǫ b > 0. We assume both vacancy and interstital are hard-core bosons, so that the allowed values for n j,a and n j,b are either 0 or 1. An exciton is described by the state with both a vacancy and an interstitial at the same lattice site n j,a = n j,b = 1. The couplings t a and t b are the hopping integrals for boson a and b, respectively, and we assume t a and t b to be positive without loss of generality. U is the on-site attractive interaction between a vacancy and an interstitial. Note that the attractive interaction between a vacancy and an interstitial reflects the strong short range repulsion between two nearby atoms when an interstitial atom is added into the lattice.
Without interaction, at U = 0, the two-species model decouples into independent vacancy and interstitial models. The ground state of the a-boson (vacancy) model is superfluid (a vacancy supersolid) if zt a > ǫ a and an empty vacuum state (insulating solid) otherwise. Similarly, the ground state of the b-boson (interstitial) model is superfluid (an interstitial supersolid) if zt b > ǫ b and the empty vacuum state (an insulating solid) otherwise.
The attractive inter-species boson interaction U couples the two types of boson, and the problem cannot be solved analytically without approximation. This model was studied by using a mean field theory at zero temperature, 2 and a special limiting case with ǫ b → ∞ but a finite ǫ b − U was investigated by a modified spin wave theory. 4 The main effect of the attractive term is to facilitate vacancy or/and interstital condensation due to excitons or bound states.
In this paper we will use QMC methods to study the phase transition and superfluid properties of the model. We use a slightly extended version of the directed loop algorithms 17,18 of the ALPS project. 19 In the stochastic series expansion (SSE) representation used by this algorithm the superfulid density can be measured through the fluctuations of the spatial winding numbers. 20 In three dimensions for a simple cubic lattice the relationship is: where α = a, b refers to the type of boson, W α is the spatial winding number of the bosons in one direction, L is the linear size of the cubic lattice and T the temperature. In addition we consider the correlated winding numbers and define We performed simulations on lattices with up to 10 3 lattice sites. Larger lattices did not equilibrate using the directed loop algorithm due to the formation of bound states between a and b bosons. A two-worm algorithm such as developed for the one-dimensional case in Ref. 21 would be required to go to larger lattices, however we found that the sizes used here were sufficient to determine the nature of the phases.

III. SYMMETRIC CASE
We first consider the symmetric case with ǫ a = ǫ b = 1. The energy cost of an exciton, (both an aand a b-boson on the same lattice site), is ∆ = ǫ a + ǫ b − U . The system becomes a trivial exciton lattice if ∆ < 0, which we will not discuss. For the symmetric case, we choose U = 1.8 in our simulations to examine correlation effects. In Fig. 1, we summarize our result by showing the phase diagram in the parameter space of temperature T and boson hopping zt = zt a = zt b . The simulations are carried out at temperatures ranging from 0.1 to 0.2, which allow us to estimate a zero temperature phase boundary at zt ≈ 0.65, smaller than the mean field value of zt ≈ 0.80, and smaller than that of the non-interacting case of U = 0 at zt = 1. Hence, quantum fluctuations which are neglected in the mean field theory further favor the superfluid phase.
We examine the temperature dependences in more detail along line A in the phase diagram of Fig. 1. Figure  2 shows the superfluid densities ρ ± and the specific heat c V as functions of T along line A in Fig. 1 (zt = 0.88) at the temperature region from T = 0.03 to T = 0.4. As T decreases, ρ ± rise abruptly below T = 0.13 and saturate to ρ ± = 0.1, c V develops a clear peak around T = 0.13, and the peak becomes sharper as the size increases. Note that ρ + = ρ − within our error bars, indicating that there are no correlations and the two types of bosons condense independently with the same superfluid density ρ a = ρ b ≈ ρ + /2.
In Fig. 3, we show the the superfluid density for system sizes L = 4, 6, 8. Finite size scaling for a second order phase transition in the U (1) universality class implies that ρ s L is a constant at the transition temperature T c . The three curves in Fig. 3 indeed cross at a single point, from which we can estimate T c = 0.129 ± 0.002. Finally we investigate the dependence on the hopping amplitude zt at various temperatures. In Fig. 4 we plot the superfluid densities and specific heat at T = 0.15 (along line C). Superfluidity develops at around zt = 0.92 as we can see from both ρ + and c V . The superfluid density as functions of zt are plotted in Fig. 5 for different system sizes along lines B, C, and D. Each set of curves cross at one point, consistent with the expected scaling at a second order phase transition.

IV. NON-SYMMETRIC CASE
We now discuss the non-symmetric case, which is more interesting and possibly more relevant to physical systems since there is a lack of vacancy-interstitial symmetry. In all the simulations reported for the non-symmetric model, we consider ǫ a = 1, ǫ b = 4, and U = 4. We use smaller values of U than in the mean-field work of Ref.
2 since larger values of U cost too much CPU time.
Our phase diagram in the parameter space of zt a and zt b at T = 0.15, is summarized in Fig. 6(a). This temperature is low enough to observe the expected supersolid The phase boundaries labeled by red circles in Fig. 6(a) are obtained from simulations on systems with up to L = 4. Calculations on larger size systems in these parameter region require much more computational effort, and are only carried out for four selected points, labeled by blue squares on the boundaries in the figure, representing typical interesting cases of the three most interesting different phase-transitions in the parameter space.
For comparison, we show in Fig. 6(b) the result of mean field calculations for the same parameters considered. Note that the QMC predicts a larger parameter space for the supersolid phases than the mean field theory, indicating again that the quantum fluctuation neglected in the mean field theory but included in the QMC is in favor of the supersolid phase. In the remaining part of this section, we discuss the phase transition as a function of temperature and as a function of boson hopping integrals.
To study the temperature dependence, we choose two typical points A (zt a = 0.3 and zt b = 3.75) and B (zt a = 0.8 and zt b = 3.5) in the parameter space as indicated in Fig. 6(a). In Fig. 7, we show the superfluid density and specific heat as functions of the temperature for the system at point A. As the temperature decreases, ρ b starts to increase sharply at around T ≈ 0.92, while ρ a remains zero. This indicates that only b−bosons condense. A scaling analysis of ρ b L gives T c = 0.924 ± 0.002.
In Fig. 8 we show ρ a,b and c V for the system at the point B where there are two transitions. As the temperature is lowered, the system first undergoes a transition at T = 0.814 ± 0.002 from a normal-solid into the V-SF(B)phase, in which the b-bosons condense. As the temperature is further lowered, the system undergoes a second phase transition at T = 0.23 ± 0.01 within the supersolid state from the V-SF(B) phase to the VI-SF-phase where both the aand b−bosons condense. The critical temperatures have again been estimated by a scaling analysis similar to the symmetric case. Note that below the lower transition point, ρ b further increases, due to the attractive interaction with the a-bosons which effectively increase the chemical potential for the b-bosons and hence their number. We have calculated ρ ± and have found that ρ + and ρ − are almost the same, so that the correlations are very small.
We now discuss the phase transitions along the lines C-F of Fig. 6(a) in more detail, at a temperature T = 0.15. There are three different phase-transitions: 1. a transition between the insulating state and the V-SF(B) phase and along the line C in Fig. 6  where ρ b becomes non-zero while ρ a remains zero (see Fig. 9). We estimate the critical value zt b = 3.44 ± 0.02 using the scaling analysis.
2. a transition between the insulating state and the VI-SF phase appears along the lines D and E in Fig. 6(a). For fixed zt a = 0.7, the superfluid densities for both a− and b-bosons are zero at small values of zt b , and become finite above a critical value, which is the same for the two types of bosons, as we can see from 3. a transition is between two supersolid phases along line F in Fig. 6(a). Along this line zt b = 3.6 and ρ b is always finite. As zt a increases ρ a is zero up to a critical value zt a = 0.55 ± 0.01. (see Fig. 10).

V. CONCLUSIONS
Our quantum Monte Carlo simulations of a two-species bosonic Hubbard model of a supersolid show a phase diagram qualitatively consistent with previous mean field results. 2 The attractive interaction between a vacancy and interstitial may facilitate the superfluidity in a bosonic solid, even when single vacancies or interstitials are gapped. Quantum fluctuations which are ignored in the mean-field calculations stabilize the superfluid phase over a larger parameter regime. Unlike the modified spin wave calculations 4 which finds first order phase transitions at finite temperatures, the quantum Monte Carlo calculations show consistency with the expected scaling behavior at second order phase transitions in the U (1) universality class, for temperatures above T = 0.10.