Equilibrium phases of two-dimensional bosons in quasi-periodic lattices
EEquilibrium phases of two-dimensional bosons in quasi-periodic lattices
C. Zhang, A. Safavi-Naini,
1, 2 and B. Capogrosso-Sansone Homer L. Dodge Department of Physics and Astronomy,The University of Oklahoma, Norman, Oklahoma ,73019, USA JILA , NIST and Department of Physics, University of Colorado, 440 UCB, Boulder, CO 80309, USA
We report on results of Quantum Monte Carlo simulations for bosons in a two dimensional quasi-periodic optical lattice. We study the ground state phase diagram at unity filling and confirm theexistence of three phases: superfluid, Mott insulator, and Bose glass. At lower interaction strength,we find that sizable disorder strength is needed in order to destroy superfluidity in favor of the Boseglass. On the other hand, at large enough interaction superfluidity is completely destroyed in favorof the Mott insulator (at lower disorder strength) or the Bose glass (at larger disorder strength). Atintermediate interactions, the system undergoes an insulator to superfluid transition upon increasingthe disorder, while a further increase of disorder strength drives the superfluid to Bose glass phasetransition. While we are not able to discern between the Mott insulator and the Bose glass atintermediate interactions, we study the transition between these two phases at larger interactionstrength and, unlike what reported in [1] for random disorder, find no evidence of a Mott-glass-likebehavior.
Introduction:
Condensed matter systems, either man-ufactured or occurring in nature, posses, in general, acertain degree of disorder. Studying physical phenomenasuch as Anderson [2] localization, resulting from the pres-ence of disorder, is therefore of crucial importance. An-derson localization pertains to the case of non-interactingfermions. More realistic systems, though, consist of in-teracting particles. For interacting systems, the inter-play between disorder and interaction may result in novelphysical effects. For instance, when random disorder isadded to paradigmatic condensed matter models, suchas the Bose-Hubbard model or the BCS model for super-conductivity, it gives rise to disorder-driven phase transi-tions from a conducting to an insulating phase, resultingfrom the localization of bosons and cooper pairs, respec-tively [3–6]. While disorder driven phase transitions havebeen observed in a wide range of experimental systemssuch as films of adsorbed He on substrates [7, 8], bosonicmagnets [9–11], and thin superconducting films [12, 13],and in spite of a remarkable theoretical effort [14–21], athorough understanding of the effects of disorder in in-teracting quantum many body systems is lacking. Onthe one hand these systems are challenging to study the-oretically, on the other poor control over experimentalcondensed matter systems does not allow for thoroughexperimental investigation of these systems.Optical lattice systems of ultra-cold atoms andmolecules provide a unique possibility of engineeringmatter with an unprecedented level of control and flex-ibility over the parameters entering the hamiltonian[22–25]. Hence, optical lattice simulators have rapidly be-come an important tool in the study of disordered sys-tems where disordered or quasi-disordered optical lat-tice potentials are created using speckle patterns ormulti-chromatic incommensurate optical lattices respec-tively [26–29]. These techniques were employed in thefirst realizations of Anderson localization in one dimen- sion using non-interacting bosons in continuum [30], andin an optical lattice [31]. Subsequently, delocalizationinduced by weak repulsive interaction was observed inone dimension [32], while localization induced by stronginteraction was demonstrated in one and three dimen-sions [28, 33, 34].In the past decade the behavior of strongly interact-ing systems in the presence of random disorder has beenstudied extensively using a variety of theoretical meth-ods [4, 35–43]. For the most part these studies have con-sidered systems of bosonic particles trapped in one-, two-,or three-dimensional optical lattices. In the absence ofdisorder, these systems feature two phases: superfluid(SF) and Mott-Insulator (MI). In the presence of ran-dom disorder a third insulating but compressible phase,known as the Bose Glass (BG), is stabilized [3]. As aresult of finite disorder strength, no direct SF-MI phasetransition exists [38] and the BG always intervenes be-tween the MI and SF regions.Optical lattice systems can also be used to create quasi-periodic trapping potentials by employing bichromaticlattices which are formed by combining two optical lat-tices with incommensurate wavelengths [27]. In solidstate physics, quasi-periodic crystalline structures suchas photonics quasicrystals were found to have a nontrivialconnection to topological states of matter [44–46]. Non-intereacting bosons in quasi-periodic one-dimensional po-tentials are described by the analytically solvable Aubry-Andr´e model [47, 48], which features Anderson localiza-tion. This model was first realized experimentally in[31], where Anderson localization was confirmed by in-vestigating transport properties, as well as the spatialand momentum distributions. As in the case of randomdisorder, the introduction of interactions to the Aubry-Andr´e model increases the complexity of the system andgives rise to new physical phenomena. Most of the re-cent theoretical work has focused on one-dimensional a r X i v : . [ c ond - m a t . d i s - nn ] N ov systems of interacting bosons in quasi-periodic poten-tials, where DMRG methods can be successfully em-ployed [49–54]. These studies have identified a directSF-MI transition of Kosterlitz-Thouless type at weak dis-order [50, 53], in accordance to the predictions of theHarris-Luck criterion[55]. The criterion states that a per-turbative quasi-periodic disorder is irrelevant from therenormalization group point of view, leaving the natureof the phase transition unchanged when compared to thetransition in the clean system.In the following, we use Path Integral QuantumMonte Carlo by the Worm algorithm [56] to study two-dimensional (2D) lattice bosons in the presence of quasi-periodic disorder. We find that if the interaction strengthis smaller than the critical interaction strength corre-sponding to the 2D SF-MI transition in the clean system,sizable disorder is needed to destroy superfluidity. On theother hand, at any given disorder, one can find an inter-action strength above which superfluidity is completelydestroyed in favor of an insulating phase. At lower disor-der strength this insulating phase is a MI, while at largerdisorder strength it is a BG. Our numerical results forthe compressibility in the range of interaction strengthswhere SF has completely disappeared are consistent witha direct MI-BG phase transition and do not show anyevidence of a cross-over region characterized by Mott-glass-like behavior (or anomalous Bose glass), unlike thefindings of [1] for the case of random disorder. Finally,at intermediate interaction strengths, the system under-goes an insulator to superfluid transition upon increasingthe strength of the disorder. One can (re)enter the Boseglass phase by further increasing of disorder strength. Hamiltonian:
We study a system of bosons in a 2D lat-tice in the presence of quasi-periodic disorder, describedby the Hamiltonian H = − J (cid:88) (cid:104) i j (cid:105) ( a † i a j + h.c ) + U (cid:88) i n i ( n i − − µ (cid:88) i n i + (cid:88) i ∆ i n i . (1)The first term in the Hamiltonian is the kinetic energy,where a † i ( a i ) are the bosonic creation (annihilation) op-erators with the usual commutation relations, and J isthe hopping matrix element between sites i and j . Weuse (cid:104) . . . (cid:105) to denote nearest neighboring sites. Here U setsthe strength of the on-site repulsion and µ is the chemicalpotential, which in the absence of disorder, sets the num-ber of particles in the system. The quasi-periodic on-sitedisorder ∆ i is created by perturbing the primary opticallattice with a second incommensurate one. The net resultis a quasi-periodic external potential that couples to theon-site density n i . Hence, the on-site disorder takes theform ∆ i = ∆ cos(2 πβ d x i + φ x ) cos(2 πβ d y i + φ y ), where ∆is the strength of disorder, φ x,y is an arbitrary phase shift,and β d measures the degree of commensurability. Both ∆ and β d can be tuned experimentally, the former bytuning the relative heights of the primary and secondarylattices and the latter by varying the wave numbers ofthe two lattice potentials. The results presented belowcorrespond to the maximally incommensurate ratio givenby the choice β d = ( √ − / Results:
In the following we present a numerical studyof the Hamiltonian (1) at unit filling ( n = N/N sites = 1)by the means of quantum Monte Carlo simulations usingthe Worm Algorithm. In order to obtain accurate resultsin the thermodynamic limit we perform finite-size scalingon the simulations results. This process is challenging inthe presence of quasi-periodic disorder where the disorderis incommensurate with the lattice. Incommensurabilitymeans that one cannot produce comparable systems bysimply scaling the lattice size. To circumvent this prob-lem we have used system sizes L , with N site = L × L ,from the Fibonacci sequence [57]. Unlike the results re-ported for one-dimensional systems [57], we have foundthat for disorder strength ∆ (cid:38) J our results dependstrongly on the choice of ( φ x , φ y ). Hence, for each set ofparameters ( µ, ∆ , U, L ) we have run simulations with fiftydifferent choices for phases φ x,y ∈ [0 , π ). The resultspresented below are extracted from the fifty runs usingthe bootstrap method. We find that further averagingover ( φ x , φ y ) realizations simply reduces the statisticalerror.The ground state phase diagram of the system at unitfilling is shown in Fig. 1, where the horizontal and ver-tical axes correspond to UJ and ∆ J respectively. At lowerdisorder strength and for U/J (cid:46)
16 the system is inthe superfluid state associated with the presence of off-diagonal long-range order. The superfluid phase is char-acterized by finite compressibility and non-zero singleparticle condensate order parameter, (cid:104) ψ (cid:105) = (cid:104) a i (cid:105) (cid:54) = 0,associated with a finite superfluid stiffness ρ s . The su-perfluid stiffness is extracted from simulations using therelation ρ s = (cid:104) W (cid:105) /dL d − β , where W is the windingnumber in space, d is the spatial dimension ( d = 2 in ourcase)[58], and β = 1 /T is the inverse temperature. In allour simulations we have chosen β such that the systemis in its ground state, and have scaled β ∝ L z where z isthe dynamical critical exponent. The SF phase becomesunstable at stronger disorder strength and a transitionto the insulating BG phase occurs. The BG phase ischaracterized by vanishing superfluid stiffness and finitecompressibility κ .For 16 (cid:46) U/J (cid:46)
35 and at low disorder strength,the system is in an insulating phase and undergoes aphase transition in favor of the SF phase upon increas-ing the disorder strength. A similar phase transition ispresent if the trapping potential features random disor-der, where it has been shown that the presence of anintervening BG phase between the MI and SF is guaran-teed by the theorem of inclusions [38]. It should be notedthat this theorem does not apply to quasi-periodic disor-der and therefore the existence of the BG phase or thelack thereof should be confirmed by direct measurementof the compressibility. However, the parameter regimecorresponding to the range of interactions and disorderstrengths where the BG region may form is narrow. Asfor the case of random disorder, the compressibility ofthe BG in narrow regions would be too small to be de-tected numerically, making it impossible to distinguishbetween the MI and BG phases. We are therefore unableto discuss the onset of the BG phase, and generically re-fer to the dashed blue region separating the SF and thezero-disorder MI in Fig. 1 as an insulating phase. Furtherincreasing the disorder strength results in the destructionof the SF order in favor of the BG.
FIG. 1. (Color online) Ground state phase diagram of the sys-tem described by Eq. 1 at filling factor n = 1. The horizontaland vertical axes are the onsite interaction strength U/J anddisorder strength ∆ /J , respectively. Using these two param-eters as tuning knobs, the system can form a Mott-Insulator(MI), a superfluid (SF), and a Bose glass (BG). Simulationsresults for the SF-insulator phase boundary are shown usingsolid orange circles (the solid orange line is a guide to theeye), while solid purple squares (the dashed line is a guide tothe eye) correspond to the phase boundary between the MIand BG phases. At lower disorder and intermediate interac-tions we are unable to distinguish between the MI and theBG (dashed blue region). Figure 2 illustrates the finite size scaling procedureused to determine the SF-BG (or generic insulator) phaseboundary (solid orange circles in Fig. 1). Here we plotthe scaled superfluid stiffness ρ s L ( d + z − with z = 2, as afunction of ∆ /J at U/J = 22 and L = 21, 34, 55, and 89(red circles, blue squares, empty black squares, and blackdiamonds, respectively). In these simulations we haveused β = ( L/ z to scale the imaginary time dimension L τ . The dynamical critical exponent z was set to d = 2,following the prediction in Ref. [3] for random disorder,and the recent unambiguous confirmation using MonteCarlo techniques [59]. The drift in the position of theintersection point indicates that a correction to the finite FIG. 2. (Color online) Main plot: Scaled superfluid stiffness ρ s L − ( d + z − with z = 2, as a function of ∆ /J for U/J = 22and L = 21, 34, 55 and 89 using red circles, blue squares,empty black squares, and black diamonds, respectively. Inthese simulations we have used β = ( L/ to scale the imag-inary time dimension L τ . Inset: Data collapse using ν = 0 . a = − . ω = − .
9, and ¯∆ c = 10 .
21 corresponding to thecritical point extracted from the main plot. Here ¯∆ = ∆ /J .The symbols are the same as those used for the main plot. size scaling relation ρ s L ( d + z − = f ( L /ν ∆ J , βL − z ) where f ( x, const) is a universal function, must be included inorder to observe data collapse. After this correction istaken into account the scaling relation takes the form ρ s L ( d + z − = (1 + aL − ω ) f ( L /ν ∆ J , βL − z )[60]. The in-set of Fig. 2 shows L ρ s / (1 + aL − ω ) as a function of( ¯∆ − ¯∆ c ) L /ν , where ¯∆ = ∆ /J . From the best datacollapse we find ν = 0 . ± .
07 ( ν = 0 .
67 holds for theSF-insulator transition of a clean system), a = − . ± . ω = − . ± .
05, and ¯∆ c = 10 . ± .
05. This value of ν suggests that the quasi-periodic disorder is still irrelevantfor ∆ /J ∼
10. As noted above, the choice of dynamicalexponent z = d has only been predicted and confirmedfor random disorder. To ensure that our choice of thecritical exponent z does not affect the position of thetransition line, we have performed finite size scaling witha choice of z = 1 . z .Finally, we turn our attention to the region U/J (cid:38) z =0.75, 1, 1.25, and 1.5 at U/J = 45. The inset of Fig. 3 shows the scaled compress-ibility κL d − z with z = 0 .
75 (top panel) and z = 1 (lowerpanel) for L =21, 34, 55, and 89 using black squares, redcircles, blue triangles and green diamonds, respectively.While we have not performed an exhaustive scan over FIG. 3. (Color online) Main plot: κ vs. ∆ /J for U/J = 45, L = 21, and β = L/
2. The compressibility becomes finite at∆ /J ∼ . /J ∼ .
5. Inset: The topand bottom panels show κL d − z versus ∆ /J for z = 0 .
75 and z = 1 respectively, at U/J = 45. The scaled compressibilityis shown for L = 21 ,
34, 55, and 89 using black squares, redcircles, blue triangles and green diamonds, respectively. Ourdata indicates that 0 . (cid:46) z (cid:46) . different values of dynamical exponent z , we find thatthe best crossing corresponds to z = 0 .
75 and gives thecritical point at ∆ c /J = 23 . ± . κ , with κ ∼ exp( − b/T α ) + c ( α < c ∼ κ away from theSF lobe boundary at fixed β as a function of ∆ /J , and atfixed ∆ /J as a function of β . The main plot of figure 3shows κ vs. ∆ /J for U/J = 45, L = 21, and β = L/ /J ∼ . /J ∼ .
5. Fig. 4 shows κ as a functionof β at U/J = 45 and L = 21 for ∆ /J =22, 23.5, 23.6,23.7, 23.8, 23.9, and 26. Our data indicates that belowthe quantum critical point, ∆ c /J = 23 . ± .
05, the sys-tem is in the MI state and κ ∼ exp( − β ∆ G ), where ∆ G is the energy gap. Upon increasing the disorder strengththe system enters the BG phase as shown by a plateauedcompressibility at large enough β (see curves correspond-ing to ∆ /J =23.9 and 26). It should be noted that at∆ /J = 23 . U/J = 45, Figure 4 suggests that it would only extend within a range of disorder strength of width ∼ FIG. 4. (Color online) The plot shows log κ as a function of β ,at U/J = 45 and L = 21 for ∆ /J =22, 23.5, 23.6, 23.7, 23.8,23.9, and 26. Below the critical point ∆ /J = 23 . ± .
05 (seeinset of Fig. 3) the behavior is consistent with that of a MI.Above the transition, the behavior is consistent with that ofthe BG phase.
In conclusion, we have used Path Integral QuantumMonte Carlo by the Worm algorithm to study the phasediagram of bosons in a two-dimensional quasi-periodicoptical lattice. As in the case of random disorder, theground state phase diagram contains three phases: su-perfluid, Mott insulator, and Bose glass. At weaker in-teractions, the superfluid phase is favored and significantdisorder has to be introduced in order to destroy super-fluidity. At strong enough interactions, the superfluidphase has disappeared, and for weak enough disorder thesystem forms a Mott insulator. Upon increasing the dis-order strength the system undergoes a phase transitionfrom Mott insulator to Bose glass. We have used fi-nite temperature simulations to establish that there is noMott-glass-like Bose glass behavior separating the Mottinsulator from the Bose glass. Finally, at intermediateinteraction strengths and lower disorder, the compress-ibility of the Bose glass is too small to be measured nu-merically in finite systems. In this region we are unableto distinguish between a Mott insulator and a Bose glass.
Acknowledgments
We would like thank N. Prokof’evand S¸.G. S¨oyler for useful discussions. This work wassupported by the National Science Foundation (NSF) un-der the grant NSF-PIF-1415561. ASN was supported byAFOSR-MURI grant. The computing for this projectwas performed at the OU Supercomputing Center forEducation & Research (OSCER) at the University of Ok-lahoma (OU). [1] Y. Wang, W. Guo, and A. W. Sandvik, arXiv.org(2011).[2] P. W. Anderson, Phys. Rev. , 1492 (1958).[3] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S.Fisher, Phys. Rev. B , 546 (1989).[4] L. Pollet, Comptes Rendus Physique , 712 (2013).[5] A. Ghosal, M. Randeria, and N. Trivedi, Phys. Rev.Lett. , 3940 (1998).[6] Y. Dubi, Y. Meir, and Y. Avishai, Nature , 876(2007).[7] P. A. Crowell, F. W. Vankeuls, and J. D. Reppy, Phys.Rev. Lett. , 1106 (1995).[8] P. S. Ebey, P. T. Finley, and R. B. Hallock, J Low Temp.Phys. , 635 (1998).[9] T. Hong, A. Zheludev, H. Manaka, and Regnault, L. P.,Phys. Rev. B , 060410 (2010).[10] R. Yu, C. F. Miclea, F. Weickert, R. Movshovich, A. Pad-uan, V. S. Zapf, and T. Roscilde, Phys. Rev. B ,134421 (2012).[11] R. Yu, L. Yin, N. S. Sullivan, J. S. Xia, C. Huan, A. Pad-uan, N. F. Oliveira, S. Haas, A. Steppke, C. F. Miclea,F. Weickert, R. Movshovich, E. D. Mun, B. L. Scott,V. S. Zapf, and T. Roscilde, Nature , 379 (2012).[12] A. M. Goldman and N. Markovic, Phys. Today , 39(1998).[13] R. Schneider, A. G. Zaitsev, D. Fuchs, and H. vonL¨ohneysen, Phys. Rev. Lett. , 257003 (2012).[14] M. Fisher, Phys. Rev. Lett. , 923 (1990).[15] M. Fisher, G. Grinstein, and S. Girvin, Phys. Rev. Lett. , 587 (1990).[16] J. D. Reppy, J Low Temp. Phys. , 205 (1992).[17] K. Moon and S. Girvin, Phys. Rev. Lett. , 1328 (1995).[18] T. Giamarchi, C. Ruegg, and O. Tchernyshyov, Nat.Phys. , 198 (2008).[19] O. Nohadani, S. Wessel, and S. Haas, Phys. Rev. Lett. , 227201 (2005).[20] F. Yamada, H. Tanaka, T. Ono, and H. Nojiri, Phys.Rev. B , 020409 (2011).[21] E. Wulf, D. H¨uvonen, J. W. Kim, A. Paduan-Filho,E. Ressouche, S. Gvasaliya, V. Zapf, and A. Zheludev,Phys. Rev. B , 174418 (2013).[22] I. Bloch, Nat. Phys. , 23 (2005).[23] O. Morsch and M. Oberthaler, Rev. Mod. Phys. , 179(2006).[24] L. Sanchez-Palencia and M. Lewenstein, Nat. Phys. , 87(2010).[25] L. Fallani, C. Fort, and M. Inguscio, Adv. At. Mol. Opt.Phy. , 119 (2008).[26] J. E. Lye, L. Fallani, M. Modugno, D. S. Wiersma,C. Fort, and M. Inguscio, Phys. Rev. Lett. , 070401(2005).[27] L. Fallani, J. E. Lye, V. Guarrera, C. Fort, and M. In-guscio, Phys. Rev. Lett. , 130404 (2007).[28] M. White, M. Pasienski, D. McKay, S. Q. Zhou,D. Ceperley, and B. DeMarco, Phys. Rev. Lett. ,055301 (2009).[29] F. Jendrzejewski, A. Bernard, K. Muller, P. Cheinet,V. Josse, M. Piraud, L. Pezze, L. Sanchez-Palencia,A. Aspect, and P. Bouyer, Nat. Phys. , 398 (2012).[30] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Cl´ement, L. Sanchez-Palencia, P. Bouyer,and A. Aspect, Nature , 891 (2008).[31] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort,M. Zaccanti, G. Modugno, M. Modugno, and M. Ingus-cio, Nature , 895 (2008).[32] B. Deissler, M. Zaccanti, G. Roati, C. D’Errico, M. Fat-tori, M. Modugno, G. Modugno, and M. Inguscio, Nat.Phys. , 354 (2010).[33] B. Gadway, D. Pertot, J. Reeves, M. Vogt, andD. Schneble, Phys. Rev. Lett. , 145306 (2011).[34] M. Pasienski, D. McKay, M. White, and B. DeMarco,Nat. Phys. , 677 (2010).[35] W. Krauth, N. Trivedi, and D. Ceperley, Phys. Rev.Lett. , 2307 (1991).[36] S. Rapsch, U. Schollwock, and W. Zwerger, Europhy.Lett. , 559 (1999).[37] J. W. Lee, M. C. Cha, and D. Kim, Phys. Rev. Lett. ,247006 (2001).[38] L. Pollet, N. Prokof’ev, B. Svistunov, and M. Troyer,Phys. Rev. Lett. , 140402 (2009).[39] V. Gurarie, L. Pollet, N. V. Prokof’ev, B. V. Svistunov,and M. Troyer, Phys. Rev. B , 214519 (2009).[40] S. G. Soyler, M. Kiselev, N. V. Prokof’ev, and B. V.Svistunov, Phys. Rev. Lett. , 185301 (2011).[41] A. E. Niederle and H. Rieger, New J Phys , 075029(2013).[42] P. Schmitteckert, T. Schulze, C. Schuster, P. Schwab,and U. Eckern, Phys. Rev. Lett. , 560 (1998).[43] F. Lin, E. S. Sørensen, and D. M. Ceperley, Phys. Rev.B , 094507 (2011).[44] Y. E. Kraus, Y. Lahini, Z. Ringel, M. Verbin, and O. Zil-berberg, Phys. Rev. Lett. , 106402 (2012).[45] Y. E. Kraus and O. Zilberberg, Phys. Rev. Lett. ,116404 (2012).[46] M. Verbin, O. Zilberberg, Y. Kraus, Y. Lahini, andY. Silberberg, Phys. Rev. Lett. , 076403 (2013).[47] G. Andre and S. Aubry, Ann. Isr. Phys. Soc. , 133(1980).[48] M. Modugno, New J Phys , 033023 (2009).[49] P. J. Y. Louis and M. Tsubota, J Low Temp. Phys. ,351 (2007).[50] G. Roux, T. Barthel, I. P. McCulloch, C. Kollath,U. Schollwock, and T. Giamarchi, Phys. Rev. A ,023628 (2008).[51] X. Deng, R. Citro, E. Orignac, and A. Minguzzi, Eur.Phys. J B , 435 (2009).[52] M. Larcher, F. Dalfovo, and M. Modugno, Phys. Rev. A , 053606 (2009).[53] X. L. Deng, R. Citro, A. Minguzzi, and E. Orignac,Phys. Rev. A , 013625 (2008).[54] X. Cai, S. Chen, and Y. Wang, Phys. Rev. A , 043613(2011).[55] J. M. Luck, Europhy. Lett. , 359 (1993).[56] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn, J.Exp. Theor. Phys. , 310 (1998).[57] X. L. Deng, R. Citro, E. Orignac, A. Minguzzi, andL. Santos, Phys. Rev. B , 195101 (2013).[58] E. L. Pollock and D. M. Ceperley, Phys. Rev. B , 8343(1987).[59] Z. Yao, K. P. C. da Costa, M. Kiselev, and N. Prokof’ev,Phys. Rev. Lett. , 225301 (2014).[60] M. Nikolaou, M. Wallin, and H. Weber, Phys. Rev. Lett.97