Numerically exact quantum gas microscopy for interacting lattice fermions
NNumerically exact quantum gas microscopy for interacting lattice fermions
Stephan Humeniuk ∗ and Yuan Wan
1, 2, 3 Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences, Beijing 100049, China Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China
A numerical method is presented for reproducing fermionic quantum gas microscope experimentsin equilibrium. By employing nested componentwise direct sampling of fermion pseudo-densitymatrices, as they arise naturally in determinantal quantum Monte Carlo (QMC) simulations, astream of pseudo-snapshots of occupation numbers on large systems can be produced. There is asign problem even when the conventional determinantal QMC algorithm can be made sign-problemfree, and every pseudo-snapshot comes with a sign and a reweighting factor. Nonetheless, this“sampling sign problem” turns out to be weak and manageable in a large, relevant parameterregime. The method allows to compute arbitrary distribution functions in Fock space and, from apractical point of view, facilitates the computation of complicated conditional correlation functions.
INTRODUCTION
The Hubbard model is a highly simplified, yet paradig-matic model of materials with strong correlations whichhas found an accurate physical realization in cold atomicgases in optical lattices [1]. Its phase diagram is stillpoorly understood, which has led to an intense synergyof numerical approaches [2, 3].Remarkably, fermionic quantum gas microscopes [4–9][see also Refs. [10, 11] and references therein] with single-site and single-atom resolution give access to the full dis-tribution function of occupation number states. This hasallowed the direct measurement of two-point correlationfunctions [12–14] and of more unconventional quantitiessuch as the full counting statistics (FCS) of macroscopicoperators [15] or the non-local string order parametercharacterizing spin-charge separation [16–18] in 1D. Con-ditional correlation functions around dopants [19, 20] andthe analysis of patterns in the snapshots [21] have re-vealed polarons in the doped Hubbard model, in and outof equilibrium [22]. Furthermore, time-dependent mea-surements give access to transport properties [23–25]. Inthis context, comparison with numerical simulations isnot only important for calibrating e.g. the temperaturein cold atoms experiments, but quite generally for re-liable benchmarking to prepare quantum simulators forparameter regimes where classical simulations are impos-sible.Yet, for fermions in D ≥ t − J model can give unbiased results for a small number ofdopants [27, 28] in the spin basis, but for the Fermi-Hubbard model at finite interaction path integral MonteCarlo simulations [29, 30] are only possible in one dimen-sion due to the fermionic sign problem which is extensivein the system size. We extend the determinantal QMC(DQMC) algorithm [31–33], by an inner loop, where Fockconfigurations are sampled directly , i.e. without autocor- reweighting factor sign reweighting spin up spin down up, down combined d hh h FIG. 1. Componentwise direct sampling in a single sampleof HS fields { s } . (a) Lattice sites (red circles) are orderedboustrophedonically and the joint distribution of their occu-pation numbers is written as a chain of conditional proba-bilities. For each sampled component (i.e. lattice site), thepseudo-probability distribution is reweighted. (b) Pseudo-snapshot generated for a given HS field configuration. Squarelattice with parameters U/t = 10 , βt = 10 , (cid:104) n (cid:105) = 1 , L = 16.Doublon-hole fluctuations, first revealed experimentally bybunching of “anti-moments” [12], are clearly visible in thissnapshot, which comes with a positive sign and a modestreweighting factor of R ≈ . relation time, from a fully tractable quasiprobability dis-tribution. A common technique for obtaining a tractablejoint probability distribution, which can be sampled di-rectly, is to model it as a product of conditional distri-butions [34], i.e. as a directed graphical model [35]. Thisidea is at the heart of autoregressive neural networks[34, 36], the generation of natural images pixel by pixel[37] and recent algorithms for simulation and generativemodelling of quantum systems [38–44].Alternative DQMC approaches, summing all Fockstates implicitly, exist for computing the FCS ofquadratic operators [45] and all elements of the reduced a r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p density matrix on small probe areas [46]. The nestedcomponentwise direct sampling technique presented hereis more versatile in that pseudo-snapshots can be pro-duced for probe areas as large as in current experimentswith the proviso that a (mild) sign problem is manage-able in experimentally relevant regimes.We calculate (i) the joint FCS of the staggered spin andpseudo-spin magnetization of the Hubbard model at halffilling, (ii) the distribution of the total number of holesand doubly occupied sites as a function of doping at hightemperature, and (iii) the magnetization environment ofa polaron, where we find qualitative agreement with arecent quantum Monte Carlo simulation for a single holein the t − J model [27]. An apparent discrepancy betweenRef. [27] and the quantum gas microscope experiment ofRef. [19], which we can also reproduce qualitatively, canbe pinpointed to a difference in doping regimes. NESTED AND COMPONETWISE DIRECTSAMPLING
We are considering the single-band Hubbard model H = − t (cid:88) (cid:104) i,j (cid:105) ,σ = ↑ , ↓ (ˆ c † i,σ ˆ c j,σ + h.c. )+ U (cid:88) i ˆ n i, ↑ ˆ n i, ↓ − µ (cid:88) i,σ ˆ n i,σ , (1)with the usual notation, and treat it within the DQMCframework [31–33]: After a Trotter-Suzuki decompo-sition of the density operator ˆ ρ ∼ exp( − βH ) into N τ = β/ ∆ τ imaginary times slices and a Hubbard-Stratonovich (HS) transformation for decoupling the in-teractions by introducing a functional integral over HSfields { s } , the density operator reads [47]ˆ ρ = 1 Z (cid:88) { s } (cid:89) σ = ↑ , ↓ (cid:32) w σ { s } e − (cid:80) i,j X σi,j ( { s } )ˆ c † i,σ ˆ c j,σ w σ { s } (cid:33) (2a) ≡ Z (cid:88) { s } (cid:89) σ = ↑ , ↓ w σ { s } ˆ ρ σ { s } . (2b)Here, formally we have exp( − X σ ) ≡ (cid:81) N τ l =1 B σl where B σl = e − ∆ τV σl ( { s l } ) e − ∆ τK is the matrix representationof the single-particle propagators for spin σ of the po-tential and kinetic part after HS transformation [33].Note that we use the conventional HS transformation[48] in which the discrete auxiliary fields couple to the S z -component of the electron spin, S zi = ˆ n i, ↑ − ˆ n i, ↓ ,which, as will be discussed below, is crucial. After in-tegrating out the fermionic degrees of freedom, w σ { s } =det (cid:0) + e − X σ ( { s } ) (cid:1) ≡ Z ( { s } ) is the contribution of thespin component σ to the Monte Carlo weight of the HSfield configuration { s } , which can also be interpreted asthe partition sum of the non-interacting fermion systemˆ ρ σ { s } . Since the kinetic and potential matrices in the ma-trix product leading to Eq. (2) do not commute, the re- sulting matrix e − Xσ ( { s } ) Z ( { s } ) is not Hermitian and (except in1D) not all diagonal matrix elements of ˆ ρ { s } , are semi-positive-definite; hence, ˆ ρ { s } is termed a pseudo -densitymatrix, while the total ˆ ρ in Eq. (2) is a true density ma-trix.The structure of the density matrix in Eq. (2) suggestsa nested sampling approach, in which the HS fields ofthe pseudo-density matrices ˆ ρ { s } are sampled using theMarkov chain of the conventional determinantal QMCalgorithm, while the occupation numbers can be sampled directly , i.e. without autocorrelation time, from each free-fermion pseudo-density matrix ˆ ρ { s } given that for fixed { s } their distribution function and all its marginals canbe calculated efficiently.From the chain rule of basic probability theory everyprobability distribution can be decomposed into a chainof conditional probabilities p { s } ( n , n , . . . , n D ) = D (cid:89) k =1 p { s } ( n k | n k − , n k − , . . . , n )(3)where some ordering of the random variables n , n , . . . , n D is implied. A sample from the jointdistibution is then generated by traversing the chainas p { s } ( n ) n ∼ p { s } ( n ) −−−−−−−−→ p { s } ( n | n ) n ∼ p { s } ( n ) −−−−−−−−→ p { s } ( n | n , n ) → . . . , where n k ∼ p { s } ( n k ) denotessampling variable n k from p { s } ( n k | n k − , n k − , . . . ) andthe sampled value is ”inserted“ into the next conditionalprobability along the chain [see Fig. 1(a)].Below we discuss how to calculate the conditionalquasiprobabilities in Eq. (3) for a free fermion pseudo-density matrix ρ σ { s } . Note that for fixed HS field con-figuration, the pseudo-density matrices for spin up anddown are statistically independent. Per HS sample, spinup and down occupancies are sampled independently andthen combined into full pseudo-snapshot, see Fig. 1(b).Henceforth, we drop the subscripts { s } and σ for nota-tional convenience. Direct sampling in the grand canonical ensemble
In the atomic microscopy, the quantity of central in-terest is the quasiprobability to find a given snapshot ofthe fermion occupation: p ( n , n , · · · n D ) = Tr( ρ Π n Π n · · · Π n D ) , (4)where n i = 0 , i . The projectors Π n i =0 = ˆ c i ˆ c † i and Π n i =1 = ˆ c † i ˆ c i projectonto the Fock states with occupation number n i . We maythink of n , n · · · n D as an ensemble of D random bi-nary variables. Provided that we can easily compute themarginals p ( n ) = Tr( ρ Π n ), p ( n , n ) = Tr( ρ Π n Π n ),etc. and thus the conditional quasiprobabilities, we maysample n , n · · · n D by a componentwise direct sampling. . . . - - . . - - . . . . . - - FIG. 2. Joint distribution P ( M z stag , Q z stag ) at half filling. (a) Weak Hubbard interaction U/t = 1 , βt = 4 , L = 12 , L A = 12.Even-odd oscillations, which are visible in the joint distribution, are smeared out in the marginal distributions (left and bottomof each panel). (b) Close to the metal-insulator crossover: U/t = 8 , βt = 5 , L = 12 , L A = 12. In the heat map outliers havebeen set to zero. (c) Strong Hubbard interaction: U/t = 14 , βt = 5 , L = 16 , L A = 8. In the grand canonical ensemble all marginal quasiprobability distributions of the occupation numberscan be computed straightforwardly, namely p ( n , n , · · · n k ) = ( − ) n + n + ··· n k det G − n G · · · G k G G − n · · · G k ... ... . . . ... G k G k · · · G kk − n k , (5)where G ij ≡ G σij ( { s } ) = (cid:104) c i,σ c † j,σ (cid:105) { s } is the equal-timesingle-paricle Green’s function of spin species σ for agiven HS field configuration { s } at a randomly chosenimaginary time slice. Eq. (5) is proven in the supplemen-tal material [49].We now decompose the high-dimensional quasiprob-ability distribution Eq. (4) into a chain of conditionalquasiprobability distributions, which by definition can becomputed as p ( n k +1 | n , n · · · n k ) = p ( n , n · · · n k , n k +1 ) p ( n , n , · · · n k ) . (6)Inserting Eq. (5) and using the determinant formula forblock matrices we find: p (0 | n , n · · · n k ) = G k +1 ,k +1 − ∆ G k +1 , (7a) p (1 | n , n · · · n k ) = 1 − G k +1 ,k +1 + ∆ G k +1 . (7b)with the “correction”∆ G k +1 = k (cid:88) i =1 G k +1 ,i ( G K,K − N K,K ) − G i,k +1 . (8)If the correction term ∆ G k +1 were zero, then the con-ditional probability p ( n k +1 | n , n , . . . ) would be simply given by the diagonal element G k +1 ,k +1 of the Green’sfunction and be independent of the other occupationnumbers[50]. Therefore the correction term is crucial forinter-site correlations. While traversing the chain of con-ditional probabilties, the block structure of the matrixwhose inverse is required in Eq. (8) can be exploited re-cursively such that no calculation of a determinant ormatrix inverse from scratch is necessary (see supplemen-tal material [49]). Reweighting
As said earlier, the pseudo-density operator ˆ ρ { s } is notHermitian and not all conditional quasiprobabilities inEq. (7) are non-negative. Therefore, we rewrite them as p ( n k ) = sign( p ( n k )) | p ( n k ) |N k × N k . (9)with the shorthand notation p ( n k ) ≡ p ( n k | n , n , . . . , n k − ). The sampling forcomponent k is then carried out using the valid prob-ability distribution q ( n k ) ≡ | p ( n k ) |N k with normalization N k = | p ( n k = 0) | + | p ( n k = 1) | [see Fig. 1(a)]. -0.22-0.2-0.18-0.16-0.14-0.12-0.1-0.08 1 1.5 2 2.5 3 3.5 4 4.5 -0.02-0.0100.010.020.030.040.050.060.5 1 1.5 2 2.5 3 3.5 4 -0.1-0.08-0.06-0.04-0.0200.020.04 0 0.5 1 1.5 2 2.5 3 3.5 4 β t = 2.5 β t = 2.0Blomquist and Carlstrom (2019) FIG. 3. Three-point spin-charge correlations C | d | ( | r | ), as depicted in the insets, around an isolated hole (see main text). U/t = 14, µ/t = − βt ∈ { , . } , system size L × L with L = 10. A total number of 10 pseudo-snapshots has been generatedwith 224 independent Markov chains. For µ/t = − L = 10, the average number of excess holes is (cid:104) N h (cid:105) − (cid:104) N d (cid:105) ≈ .
35 (seeFCS of N h and N d in [49]). The data agrees qualitatively with data for comparable parameters of temperature and interactions( βt = 2 . βJ = 0 .
66) from Ref. [27] where a single hole in the t − J model was simulated. Having sampled the entire chain of conditional proba-bilities for both spin components, the generated snap-shot is associated with a (signed) reweighting factor R = R ↑ R ↓ where R σ = sign( w σ { s } ) D (cid:89) k =1 (cid:0) sign( p { s } ( n σk )) N σk (cid:1) , (10)and all quantities O ( n ) = O ( n , n , . . . , n D ) that areevaluated on M generated snapshots { n i } Mi =1 need to bereweighted as (cid:104)O(cid:105) = (cid:80) Mi =1 R i O ( n i ) (cid:80) Mi =1 R i . (11)The joint pseudo-probabilitiy p ( n , n , . . . ) can be fac-tored in arbitrary order into components. However, thereweighted distribution q ( n , n , . . . ) and thus the mag-nitude and sign of the reweighting factor depends on thechosen factor ordering, leaving room for optimization ina given HS sample.The severity of the sign problem is basis dependent:We find that it strongly depends on the single-particlebasis for sampling and the chosen HS transformation.Sampling in the S x - or S y -basis of the electron spin (or inthe S z -basis in momentum space) leads to a very strongsign (phase) problem. Chosing the HS transformationthat couples to the electron charge density [33] ratherthan the electron spin also gives a very severe sign prob-lem.Even at half filling, where the DQMC algorithm can bemade sign-problem free, there is a sampling sign problem,with a non-uniform dependence on U/t . The average signdiminishes as
U/t increases from
U/t = 0 and reachesa minimum at an intermediate value (
U/t ) MIC ≈ −
7, where a metal-to-insulator crossover occurs [51]. For
U > ( U/t ) MIC the sign problem again gradually becomesmuch less severe. The dependence of the sampling sign problem on temperature, interaction strength and dopingis presented in [49].In all simulations, we use a Trotter discretization of∆ τ t = 0 .
02 and generate around 20 pseudo-shapshotsper HS sample on equidistant imaginary time slices. Thesimulation code has been verified by comparing with ex-act diagonalization results [49].
APPLICATIONSJoint full counting statistics (FCS)
At half filling, the Hubbard model has an enlarged( SU (2) × SU (2)) /Z = SO (4) symmetry [52], whichis the combination of spin-rotational and particle-holesymmetry and is generated by two commuting sets ofangular momentum operators describing the total spinand total pseudo-spin of the system, respectively. Asthe temperature is lowered, domains form with an or-der parameter of the same symmetry, which, apartfrom fluctuating in length, can rotate [15, 45] on an SO (4) sphere between antiferromagnetic, s -wave pair-ing and charge-density wave correlations, as one goesfrom one domain to a neighbouring domain. A jointhistogram of the operators representing different com-ponents of the order parameter should reflect that theyare projections of the same vector along different direc-tions in order parameter space. Fig. 2 shows the jointdistribution of the projections of the staggered mag-netization M z stag = (cid:80) N A i =1 ( − i (ˆ n i, ↑ − ˆ n i, ↓ ) and the stag-gered pseudo-spin Q z stag = (cid:80) N A i =1 ( − i (ˆ n i, ↑ + ˆ n i, ↓ −
1) ona square probe area of size N A = L A , as obtained fromre-weighted pseudo-snapshots. As U/t increases [Fig. 2(a-c)], the suppression of charge fluctuations manifestsitself in the narrowing of the pseudo-spin distribution.The even-odd effect visible in the joint distributions isdue to the fact that an even number of sites can only ac-comodate an even magnetization of spin- (pseudo-spin)moments.Note that P ( M z stag , Q z stag ) could also be obtained us-ing the generating function approach of Ref. [45]. Thisis not true for the FCS of the number of doublons N d = (cid:80) N sites ,A i =1 ˆ n i, ↑ ˆ n i, ↓ and the total number of holes N h = (cid:80) N sites ,A i =1 (1 − ˆ n i, ↑ ) (1 − ˆ n i, ↓ ) since these operatorsare non-quadratic in fermionic operators. We find thate.g. for U/t = 12 and βt = 4, the FCS of N h and N d asa function of doping are accurately modelled by (shifted)binomial distributions [49]. Three-point spin-charge correlator
A three-point spin-charge correlator in the referenceframe of the hole [19] was calculated for a single hole inthe t − J model at zero temperature with DMRG and trialwave functions [53] and at finite (high and low) temper-ature using worm algorithm QMC [27]. Presumably thesame correlator was measured experimentally in Ref. [19]for the Hubbard model at large interactions. For the pur-pose of meaningful comparison between the Hubbard and t − J model we calculate the slightly modified correlator: C | d | ( r ) = (cid:104)P h r S z r + r + d S z r + r − d (cid:105)(cid:104)P h r (cid:105) (12)where S z r = ˆ n r , ↑ − ˆ n r , ↓ , and the projector P h r =ˆ n h r ˆ n s r +ˆ e x ˆ n s r − ˆ e x ˆ n s r +ˆ e y ˆ n s r − ˆ e y × ˆ n s r +ˆ e x +ˆ e y ˆ n s r +ˆ e x − ˆ e y ˆ n s r − ˆ e x − ˆ e y ˆ n s r − ˆ e x +ˆ e y (13)with ˆ n h r = (1 − ˆ n r , ↑ )(1 − ˆ n r , ↓ ) and ˆ n s r = (1 − ˆ n r , ↑ )ˆ n r , ↓ +(1 − ˆ n r , ↓ )ˆ n r , ↑ selects configurations with an empty siteat position r surrounded from eight sides by spin-onlystates, which serves to exclude nearest- and next-nearestneighbour doublon-hole pairs from the statistics. Theconditional correlation functions in the reference frameof the hole, Eq. (12), are implemented straightforwardlyby applying a filter to the pseudo-snapshots, whereas im-plementing their evaluation using Wick’s theorem wouldbe a little bit more cumbersome (although this wouldgive better statistics since all Fock states are summedimplicitly).Fig. 3 shows overall qualitative agreement of C | d | ( r )for our data for the Hubbard model and that of Ref. [27]for the t − J model, with some notable differences in C √ and C in the immediate vicinity of the hole. Acareful comparison of the t − J model with the Hubbardmodel would require a renormalization of all correlatorsin the former by a polynomial in t/U [54] (although cer-tain qualitative model differences may not be capturedperturbatively [55]). There are qualitative differences to the experimentaldata of Ref. [19], which were already noted in Ref. [27].The data of Ref. [19] can be reproduced by nested compo-nentwise sampling, pointing, however, to a very differentconclusion: namely, that magnetic polarons have disap-peared for the relatively high doping level of Ref. [19] (see[49]). CONCLUSION
In conclusion, the presented method for generatingpseudo-snapshots allows both theorist to take part in theexploration of fermionic quantum microscopy and exper-imentalists to use numerical simulations in a more ver-satile way [56]. The nested componentwise direct sam-pling method is generic to all fermionic Monte Carlomethods that are based on the “free fermion decomposi-tion” [47] and is easily adapted to projector DQMC [33]for accessing zero-temperature properties where Slaterdeterminants rather than thermal free-fermion pseudo-density matrices are sampled directly in the inner loopof the Markov chain. For Hubbard models around in-termediate interaction strength
U/t ≈ ( U/t ) MIC furtherimprovements are required to reduce the sampling signproblem. There, the interest lies in potentially observingnon-Gaussian fluctuations [57] or characterizing attrac-tion and spin-correlations between dopants [28] in thecrossover from a polaronic metal to a Fermi liquid [20].While a general solution of the sampling sign problem isunlikely, it remains to be investigated how more generalHS decouplings and representations of the electron oper-ator [58], which were successful at eliminating the signproblem of the Monte Carlo weights, affect the samplingsign problem and whether snapshots in another single-particle basis can be generated.
ACKNOWLEDGMENTS
We thank Xiaopeng Li and Lei Wang for motivatingdiscussions and acknowledge Emil Blomquist and JohanCarlstr¨om for providing the raw data for Fig. 3 as well ascomments on the manuscript. This work is supportedby the International Young Scientist Fellowship fromthe Institute of Physics, Chinese Academy of Sciences,Grant No. 2018004 (Humeniuk) and by the NationalScience Foundation of China, Grant No. 11974396, andthe Strategic Priority Research Program of the ChineseAcademy of Sciences, Grant No. XDB33020300 (Wan).The simulations were carried out on TianHe-1A at theNational Supercomputer Center in Tianjin, China. ∗ [email protected][1] C. Gross and I. Bloch, Science , 995 (2017).[2] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik,G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M.Henderson, C. A. Jim´enez-Hoyos, E. Kozik, X.-W. Liu,A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria,H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn,S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull(Simons Collaboration on the Many-Electron Problem),Phys. Rev. X , 041041 (2015).[3] T. Sch¨afer, N. Wentzell, I. ˇSimkovic, Fedor, Y.-Y.He, C. Hille, M. Klett, C. J. Eckhardt, B. Arzhang,V. Harkov, F.-M. Le R´egent, A. Kirsch, Y. Wang, A. J.Kim, E. Kozik, E. A. Stepanov, A. Kauch, S. Ander-gassen, P. Hansmann, D. Rohe, Y. M. Vilk, J. P. F.LeBlanc, S. Zhang, A. M. S. Tremblay, M. Ferrero,O. Parcollet, and A. Georges, arXiv:2006.10769 [cond-mat.str-el].[4] L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V.Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein,Phys. Rev. Lett. , 193001 (2015).[5] E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaude-cerf, G. D. Bruce, and S. Kuhr, Nature Physics , 738(2015).[6] M. F. Parsons, F. Huber, A. Mazurenko, C. S.Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, andM. Greiner, Phys. Rev. Lett. , 213002 (2015).[7] A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Sa-lomon, I. Bloch, and C. Gross, Phys. Rev. Lett. ,263001 (2015).[8] G. J. A. Edge, R. Anderson, D. Jervis, D. C. McKay,R. Day, S. Trotzky, and J. H. Thywissen, Phys. Rev. A , 063406 (2015).[9] P. T. Brown, D. Mitra, E. Guardado-Sanchez, P. Schauß,S. S. Kondov, E. Khatami, T. Paiva, N. Trivedi, D. A.Huse, and W. S. Bakr, Science , 1385 (2017).[10] T. Hartke, B. Oreg, N. Jia, and M. Zwierlein, Phys. Rev.Lett. , 113601 (2020).[11] J. Koepsell, S. Hirthe, D. Bourgund, P. Sompet, J. Vi-jayan, G. Salomon, C. Gross, and I. Bloch, Phys. Rev.Lett. , 010403 (2020).[12] L. W. Cheuk, M. A. Nichols, K. R. Lawrence, M. Okan,H. Zhang, E. Khatami, N. Trivedi, T. Paiva, M. Rigol,and M. W. Zwierlein, Science , 1260 (2016),arXiv:1606.04089 [cond-mat.quant-gas].[13] M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo,L. Pollet, I. Bloch, and C. Gross, Science , 1257(2016), arXiv:1605.05661 [cond-mat.quant-gas].[14] M. F. Parsons, A. Mazurenko, C. S. Chiu, G. Ji,D. Greif, and M. Greiner, Science , 1253 (2016),arXiv:1605.02704 [cond-mat.quant-gas].[15] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons,M. Kan´asz-Nagy, R. Schmidt, F. Grusdt, E. Demler,D. Greif, and M. Greiner, Nature (London) , 462(2017).[16] T. A. Hilker, G. Salomon, F. Grusdt, A. Omran, M. Boll,E. Demler, I. Bloch, and C. Gross, Science , 484(2017), arXiv:1702.00642 [cond-mat.quant-gas].[17] G. Salomon, J. Koepsell, J. Vijayan, T. A. Hilker, J. Ne-spolo, L. Pollet, I. Bloch, and C. Gross, Nature ,5660 (2018). [18] J. Vijayan, P. Sompet, G. Salomon, J. Koepsell,S. Hirthe, A. Bohrdt, F. Grusdt, I. Bloch, andC. Gross, Science , 186 (2020), arXiv:1905.13638[cond-mat.quant-gas].[19] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A.Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross,Nature , 358362 (2019).[20] J. Koepsell, D. Bourgund, P. Sompet, S. Hirthe,A. Bohrdt, Y. Wang, F. Grusdt, E. Demler, G. Sa-lomon, C. Gross, et al. , (2020), arXiv:2009.04440 [cond-mat.quant-gas].[21] C. S. Chiu, G. Ji, A. Bohrdt, M. Xu, M. Knap, E. Dem-ler, F. Grusdt, M. Greiner, and D. Greif, Science ,251256 (2019).[22] G. Ji, M. Xu, L. H. Kendrick, C. S. Chiu, J. C.Br¨uggenj¨urgen, D. Greif, A. Bohrdt, F. Grusdt, E. Dem-ler, M. Lebrat, et al. , (2020), arXiv:2006.06672 [cond-mat.quant-gas].[23] M. A. Nichols, L. W. Cheuk, M. Okan, T. R. Hartke,E. Mendez, T. Senthil, E. Khatami, H. Zhang, andM. W. Zwierlein, Science , 383 (2019).[24] P. T. Brown, D. Mitra, E. Guardado-Sanchez,R. Nourafkan, A. Reymbaut, C.-D. H´ebert, S. Bergeron,A. M. S. Tremblay, J. Kokalj, D. A. Huse, P. Schauß, andW. S. Bakr, Science , 379 (2019), arXiv:1802.09456[cond-mat.quant-gas].[25] R. Anderson, F. Wang, P. Xu, V. Venu, S. Trotzky,F. Chevy, and J. H. Thywissen, Phys. Rev. Lett. ,153602 (2019).[26] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn,Journal of Experimental and Theoretical Physics , 310(1998).[27] E. Blomquist and J. Carlstr¨om, (2019),arXiv:1912.08825 [cond-mat.str-el].[28] E. Blomquist and J. Carlstr¨om, (2020),arXiv:2007.15011 [cond-mat.quant-gas].[29] J. E. Hirsch, R. L. Sugar, D. J. Scalapino, andR. Blankenbecler, Phys. Rev. B , 5033 (1982).[30] J. E. Hirsch, Phys. Rev. B , 3216 (1986).[31] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar,Phys. Rev. D , 2278 (1981).[32] E. Loh Jr. and J. Gubernatis, in Electronic Phase Tran-sitions , Modern Problems in Condensed Matter Sciences,Vol. 32, edited by W. Hanke and Y. Kopaev (North-Holland, Amsterdam, 1992) Chap. 4, pp. 177–235.[33] F. F. Assaad, in
Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms, PublicationSeries of the John von Neumann Institute for Computa-tion (NIC) (Edited by J. Grotendorst, D. Marx, and A.Muramatsu (NIC, Jlich, 2002)).[34] H. Larochelle and I. Murray, in
Proceedings of the Four-teenth International Conference on Artificial Intelligenceand Statistics (2011) pp. 29–37.[35] J. Pearl,
Probabilistic reasoning in intelligent systems:networks of plausible inference (Elsevier, 2014).[36] B. Uria, M.-A. Cˆot´e, K. Gregor, I. Murray, andH. Larochelle, The Journal of Machine Learning Research , 7184 (2016), arXiv:1605.02226 [cs.LG].[37] A. van den Oord, N. Kalchbrenner, and K. Kavukcuoglu,arXiv:1601.06759 [cs.CV].[38] O. Sharir, Y. Levine, N. Wies, G. Carleo, andA. Shashua, Phys. Rev. Lett. , 020503 (2020).[39] D. Wu, L. Wang, and P. Zhang, Phys. Rev. Lett. ,080602 (2019). [40] Z. Wang and E. J. Davis, (2020), arXiv:2003.01358[quant-ph].[41] A. J. Ferris and G. Vidal, Phys. Rev. B , 165146(2012).[42] Z.-Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang,Phys. Rev. X (2018).[43] P. Clifford and R. Clifford, in Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algo-rithms (SIAM, 2018) pp. 146–155.[44] X. Li, G. Zhu, M. Han, and X. Wang, Physical ReviewA (2019).[45] S. Humeniuk and H. P. B¨uchler, Phys. Rev. Lett. ,236401 (2017), arXiv:1706.08951 [cond-mat.str-el].[46] S. Humeniuk, Phys. Rev. B , 115121 (2019),arXiv:1907.12434 [cond-mat.str-el].[47] T. Grover, Phys. Rev. Lett. , 130402 (2013).[48] J. E. Hirsch, Phys. Rev. B , 4059 (1983).[49] Supplemental material.[50] An approximation which considers only the diagonal el-ements [59] of the Green’s function may give a particlenumber distribution with the correct average and vari-ance, but all correlations between sites in a given HSsample will be lost.[51] A. J. Kim, F. ˇSimkovic, and E. Kozik, Phys. Rev. Lett. , 117602 (2020).[52] C. N. Yang and S. C. Zhang, Mod. Phys. Lett. B , 759(1990).[53] F. Grusdt, A. Bohrdt, and E. Demler, Phys. Rev. B ,224422 (2019).[54] J.-Y. P. Delannoy, M. J. P. Gingras, P. C. W. Holdsworth,and A.-M. S. Tremblay, Phys. Rev. B , 115114 (2005).[55] T.-P. Choy and P. Phillips, Phys. Rev. Lett. (2005).[56] An implementatoin of the algorithm will be made avail-able online.[57] M. Moreno-Cardoner, J. F. Sherson, and G. DeChiara, New Journal of Physics , 103015 (2016),arXiv:1510.05959 [cond-mat.stat-mech].[58] Z.-X. Li and H. Yao, Annual Review of Condensed Mat-ter Physics , 337356 (2019).[59] E. Khatami, E. Guardado-Sanchez, B. M. Spar, J. F.Carrasquilla, W. S. Bakr, and R. T. Scalettar, (2020),arXiv:2002.12310 [cond-mat.str-el]. SUPPLEMENTAL MATERIALNumerically exact quantum gas microscopy for interacting lattice fermions
Stephan Humeniuk , Yuan Wan Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences, Beijing 100049, China Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China
Here we provide a proof of a key equation in the main text, details for an efficient implementation of the algorithm,a parameter scan of the sampling sign problem, and additional data on the polaron problem. The sections of thesupplemental material are organized as follows.1. Inductive proof of the expression for the marginals of the joint quasiprobabilities Eq. (5) in the main text.2. Exploiting the block structure in Eq. (8) of the main text.3. Parameter dependence of the sampling sign problem at and away from half filling.4. Code verification using exact diagonalization on a small system and by symmetry inspection of spin-only mi-crostates around an isolated hole on a large system.5. Additional data for the polaron problem and comparison with the data of the quantum gas microscope experi-ment of Ref. [19].
INDUCTIVE PROOF OF EQ. (5) IN THE MAIN TEXT
Let G (0) i α ,j β = (cid:104) ˆ c i α ˆ c † j β (cid:105) (S1)be the single-particle Green’s function of a free fermion system, and α, β ∈ { , . . . , D } . First we prove the well-knownfact that Wick’s theorem for higher-order correlation functions can be expressed in the compact determinant form (cid:68)(cid:16) ˆ c i ˆ c † j (cid:17) (cid:16) ˆ c i ˆ c † j (cid:17) · · · (cid:16) ˆ c i n ˆ c † j n (cid:17)(cid:69) = det (cid:16) G (0) I = { i ,i ,...,i n } ; J = { j ,j ,...,j n } (cid:17) ≡ (cid:104) G (0) (cid:105) I ; J . (S2)We use the convention that [ A ] I ; J refers to the determinant of the submatrix of A whose row index (column index)runs in the set I ( J ), which is the “inclusive” definition of the minor.The proof goes by induction. The case n = 1 is true by virtue of the definition (S1). In the induction step, we useWick’s theorem, writing all non-vanishing contractions for a product of n + 1 pairs of fermionic operators as (cid:68)(cid:16) ˆ c i ˆ c † j (cid:17) · · · (cid:16) ˆ c i n +1 ˆ c † j n +1 (cid:17)(cid:69) = − G (0) i ,j n +1 (cid:68)(cid:16) ˆ c i n +1 ˆ c † j (cid:17) (cid:16) ˆ c i ˆ c † j (cid:17) · · · (cid:16) ˆ c i n ˆ c † j n (cid:17)(cid:69) − n − (cid:88) k =2 G (0) i k ,j n +1 (cid:68)(cid:16) ˆ c i ˆ c † j (cid:17) (cid:16) ˆ c i ˆ c † j (cid:17) · · · (cid:16) ˆ c i n +1 ˆ c † j k (cid:17) (cid:16) ˆ c i k +1 ˆ c † j k +1 (cid:17) · · · (cid:16) ˆ c i n ˆ c † j n (cid:17)(cid:69) − G (0) i n ,j n +1 (cid:68)(cid:16) ˆ c i ˆ c † j (cid:17) · · · (cid:16) ˆ c i n − ˆ c † j n − (cid:17) (cid:16) ˆ c i n +1 ˆ c † j n +1 (cid:17)(cid:69) + G (0) i n +1 ,j n +1 (cid:68)(cid:16) ˆ c i ˆ c † j (cid:17) · · · (cid:16) ˆ c i n ˆ c † j n (cid:17)(cid:69) . (S3)The minus-sign, e.g. in the second line of Eq. (S3), comes from the permutation of ˆ c i n +1 with ˆ c † j k . Using the inductionhypothesis Eq. (S2), valid for n pairs of fermionic operators, the remaining correlators in Eq. (S3) can be expressedas determinants (cid:68)(cid:16) ˆ c i ˆ c † j (cid:17) · · · (cid:16) ˆ c i n +1 ˆ c † j n +1 (cid:17)(cid:69) = − G (0) i ,j n +1 · det G (0) i n +1 ,j G (0) i n +1 ,j · · · G (0) i n +1 ,j n G (0) i ,j G (0) i ,j · · · G (0) i ,j n ... ... . . . ... G (0) i n ,j G (0) i n ,j · · · G (0) i n ,j n − G (0) i ,j n +1 · det G (0) i ,j G (0) i ,j · · · G (0) i ,j n G (0) i n +1 ,j G (0) i n +1 ,j · · · G (0) i n +1 ,j n ... ... . . . ... G (0) i n ,j G (0) i n ,j · · · G (0) i n ,j n − · · ·− G (0) i n ,j n +1 · det G (0) i ,j G (0) i ,j · · · G (0) i ,j n ... ... . . . ... G (0) i n − ,j G (0) i n − ,j · · · G (0) i n − ,j n G (0) i n +1 ,j G (0) i n +1 ,j · · · G (0) i n +1 ,j n + G (0) i n +1 ,j n +1 · det G (0) i ,j G (0) i ,j · · · G (0) i ,j n G (0) i ,j G (0) i ,j · · · G (0) i ,j n ... ... . . . ... G (0) i n ,j G (0) i n ,j · · · G (0) i n ,j n . (S4)By rearranging rows in (S4), it can be concluded that (S4) is the expansion of the ( n + 1) × ( n + 1) determinantdet G (0) i ,j G (0) i ,j · · · G (0) i ,j n +1 G (0) i ,j G (0) i ,j · · · G (0) i ,j n +1 ... ... . . . ... G (0) i n +1 ,j G (0) i n +1 ,j · · · G (0) i n +1 ,j n +1 (S5)along the last column according to Laplace’s formula. This completes the inductive proof. Except for the lastdeterminant in Eq. (S4), row exchanges are necessary to obtain the correct submatrix structure. In the determinantaccompanying the single-particle Green’s function G (0) i k ,j n +1 in Eq. (S4) we need to perform ( n − k ) row exchangeswhich results in a factor ( − n − k such that the total sign is ( − n − k +1 , which is identical to the alternating factor( − n +1+ k coming from Laplace’s formula.Eq. (5) differs from Eq. (S2) in that instead of pairings (ˆ c i k ˆ c † j k ) there are projectors of the formΠ k = n k ˆ c † j k ˆ c j k + (1 − n k )ˆ c j k ˆ c † j k . (S6)Depending on whether n k = 0 or n k = 1, only one of either terms in (S6) applies. Let us replace for the momentonly the k -th pairing (ˆ c i k ˆ c † j k ) in Eq. (S2) by Π k . If n k = 0, the resulting expression is covered by Eq. (S2). The case n k = 1 is different since ˆ c † j k is to the right of ˆ c j k . Using ˆ c † j k ˆ c j k = 1 − ˆ c j k ˆ c † j k , one obtains (cid:104) (ˆ c i ˆ c † j · · · Π k · · · (ˆ c i n ˆ c † j n ) (cid:105) = ( − n k n k (cid:104) det (cid:16) G (0) I \ j k ,J \ j k (cid:17) − det( G (0) I,J ) (cid:105) (S7)To make the connection with Eq. (5) we consider the matrix (cid:16) ˜ G (0) (cid:17) ij = (cid:16) G (0) (cid:17) ij − n k δ i,j k δ j,j k (S8)and develop its determinant with respect to the k -th column:det( ˜ G ) = n (cid:88) l =1 ( − l + k (cid:16) ˜ G (0) (cid:17) i l ,j k (cid:104) ˜ G (0) (cid:105) I \ i l ,J \ j k (S9)= n (cid:88) l =1 ( − l + k (cid:16) G (0) (cid:17) i l ,j k (cid:104) G (0) (cid:105) I \ i l ,J \ j k − n k ( − k (cid:104) G (0) (cid:105) I \ j k ,J \ j k (S10)Here, angular braces ( · ) denote matrix elements and square brackets [ · ] denote the inclusive definition of the minor.Note that the minors (cid:104) ˜ G (0) (cid:105) I \ i l ,J \ j k = (cid:2) G (0) (cid:3) I \ i l ,J \ j k since ˜ G (0) and G (0) only differ in the element ( j k , j k ) which isexcluded from the minors. Now, one can recognize in the first sum Eq. (S10) the Laplace expansion of G (0) withrespect to the k -th column and undo it again to recover det( G (0) ):det (cid:16) ˜ G (0) (cid:17) = det (cid:16) G (0) (cid:17) − n k (cid:104) G (0) (cid:105) I \ j k ,J \ j k . (S11)Eq. (S11) is identical to Eq. (S7), which proves that (cid:104) (ˆ c i ˆ c † j ) · · · Π k · · · (ˆ c i n ˆ c † j n ) (cid:105) = det G (0) i ,j G (0) i ,j · · · · · · G (0) i ,j n G (0) i ,j G (0) i ,j · · · · · · G (0) i ,j n ... ... . . . · · · ...... ... G (0) j k ,j k − n k · · · ...... ... ... . . . ... G (0) i n ,j G (0) i n ,j · · · · · · G (0) i n ,j n . (S12)Repeatedly replacing each pairing (ˆ c i l ˆ c † j l ) in Eq. (S2) by a projector Π l of the form (S6) and repeating the derivationfrom (S8) to (S12), with a Laplace expansion carried out with respect to the l -th column, completes the proof ofEq. (5). EXPLOITATION OF BLOCK MATRIX STRUCTURE
In the following, it is understood that G ≡ G (0) denotes a free-fermion Green’s function, and we drop the superscript.Eq. (5) implies that the expressions for the joint quasiprobability distributions of successive numbers of componentsare related by a block matrix structure. Using the formula for the determinant of a block matrix and noticing that G ( y, y ) is just a number: p ( x , x , . . . , x k − ; y ) = 1 Z det X k − G ( x , y ) G ( x , y )... G ( x k − , y ) G ( y, x ) G ( y, x ) · · · G ( y, x k − ) G ( y, y ) (S13)= 1 Z det( X k − ) G ( y, y ) − k − (cid:88) i,j =1 G ( y, x i ) (cid:2) X − k − (cid:3) i,j G ( x j , y ) . (S14)Given that X k − is itself a block matrix X k − = X k − G ( x , x k − ) G ( x , x k − )... G ( x k − , x k − ) G ( x k − , x ) G ( x k − , x ) · · · G ( x k − , x k − ) G ( x k − , x k − ) (S15)with G ( x k − , x k − ) just a number and assuming that the inverse X − k − is already known, one can make use of theformula for the inversion of a block matrix to compute the inverse of X k − in an economical way. We define g ≡ G ( x k − , x k − ) − k − (cid:88) i,j =1 G ( x k − , x i ) (cid:2) X − k − (cid:3) i,j G ( x j , x k − ) (S16)and recognize that g = p ( x , x , . . . , x k − ) p ( x , x , . . . , x k − ) = p ( x k − | x k − , . . . , x , x ) , (S17)which means that we have computed g already previously when sampling the ( k − X − k − = (cid:32) X − k − + g − (cid:126)u ⊗ (cid:126)v T − g − (cid:126)u − g − (cid:126)v T g − (cid:33) , (S18)where [ (cid:126)u ] i = k − (cid:88) j =1 (cid:2) X − k − (cid:3) ij G ( x j , x k − ) , (S19)[ (cid:126)v T ] j = k − (cid:88) i =1 (cid:2) X − k − (cid:3) ij G ( x k − , x i ) , (S20)and (cid:2) (cid:126)u ⊗ (cid:126)v T (cid:3) ij = [ (cid:126)u ] i [ (cid:126)v T ] j . (S21)Thus, the update X − k − → X − k − requires the computation of (cid:126)u , (cid:126)v T , and the exterior product (cid:126)u ⊗ (cid:126)v T , which is oforder O (( k − ). It is easy to see that the sampling of N p particle positions requires O ( (cid:80) N p i =1 i ) = O ( N p ) floatingpoint operations. It is not necessary to compute any inverse or determinant from scratch. SIGN PROBLEM FOR NESTED COMPONENTWISE DIRECT SAMPLING AT AND AWAY FROMHALF FILLING
The severity of the conventional sign problem in the determinantal QMC algorithm is measured by the average signof the Monte Carlo weight (cid:104) s (cid:105) = 1 N HS samples (cid:88) { s } sign( w ↑{ s } w ↓{ s } ) , (S22)where the sum is over auxiliary field configurations of the Hubbard-Stratonovich samples. To quantify the “samplingsign problem” in the nested componentwise sampling algorithm we introduce the average sign of a snapshot, (cid:104) ˜ s (cid:105) ,which is given by (cid:104) ˜ s (cid:105) = 1 M (cid:88) { s } sign( w ↑{ s } w ↓{ s } ) N snapshots per HS (cid:88) i =1 (cid:89) σ = ↑ , ↓ D (cid:89) k =1 sign( p ( i ) { s } ( n σk )) , (S23)where M = N HS samples × N snapshots per HS is the total number of generated snapshots. Here, it is understood thatsnapshots drawn from the same HS sample, the number of which is N snapshots per HS , are multiplied by the sign of thecorresponding Monte Carlo weight w ↑{ s } w ↓{ s } , in case that there is already a conventional sign problem at the levelof the determinantal QMC algorithm. Note that within one HS sample, { s } , snapshots for spin- ↑ and spin- ↓ can bepaired up arbitrarily into a full snapshot due to the statistical independence of the two spin species. β t=4, L A =12 β t=10, L A =12 β t=5, L A =16 β t=10, L A =16101001000100001000001e+061e+07 0 2 4 6 8 10 12 14av.av. max. FIG. S1. Sampling sign problem at half filling. The probe area is the total L × L square system ( L = L A ). (a) Average signof a pseudo-snapshot. (b) Average (unsigned) reweighting factor (full lines) and average maximum reweighting factor (dashedlines). The red dashed-dotted line indicates a rough estimate of the threshold of maximum reweighting factors below whichnumerically exact results can be obtained with modest computational effort. Furthermore, one can define the average (unsigned) reweighting factor R av. = 1 M (cid:88) { s } N snapshots per HS (cid:88) i =1 (cid:89) σ = ↑ , ↓ D (cid:89) k =1 N σ, ( i ) k, { s } (S24)and the average maximum reweighting factor, averaged over independent Markov chains, R av. max. = 1 N Markov chains N Markov chains (cid:88) m =1 max snapshots i from m-th Markov chain ( | R i | ) . (S25)The average maximum reweighting R av. max. is the most relevant indicator of the severity of the sign problem as itquantifies the ”inflation“ of value of an individual sample. Because the reweighting factor can fluctuate over severalorders of magnitude, the average sign alone is not sufficient for a characterization of the ”sampling sign problem“.The snapshots generated in quantum gas microscope experiments originate from independent experimental runswhereas the pseudo-snapshots of the nested componentwise sampling are affected by the autocorrelation time inherentin the sampling of Hubbard-Stratonovich field configurations via the standard determinantal QMC algorithm. Forlarge U/t , pseudo-snapshots taken from the same HS sample (and at the same imaginary time slice), differ mostly inthe positions of holes and doubly occupied sites, while the spin background stays largely fixed.Additionally, pseudo-snapshots come with a sign and reweighting factor due to the non-Hermiticity of the pseudo-density operator within each HS sample. Therefore, experimental snapshots and pseudo-snapshots are not directlycomparable. An important question is how many pseudo-snapshots M are required to obtain a comparable precisionof measurement quantities as from M (exp.) experimental snapshots. Taking into account the average maximumreweighting factor, the effective number of pseudo-snapshots M (cid:48) that is equivalent to M ( exp ) can be roughly estimated -0.200.20.40.60.811.2 0.2 0.4 0.6 0.8 1 U/t=4, β =4U/t=5, β =4U/t=6, β =4U/t=8, β =4 U/t=8, β =2U/t=10, β =2U/t=12, β =2U/t=14, β =2 FIG. S2. Conventional sign problem and sampling sign problem away from half filling: dependence on
U/t . (a,c) Average sign( (cid:104) s (cid:105) , full lines) and average sampling sign ( (cid:104) ˜ s (cid:105) , dashed lines); (b,d) average reweighting factor (full lines) and average maximumreweighting factor (dotted lines). The inverse temperature β is in units of 1 /t . The probe area is equal to the system size with L A = L = 8. -0.200.20.40.60.811.2 0.2 0.4 0.6 0.8 1 U/t=8, β =2U/t=8, β =3U/t=8, β =4 U/t=14, β =1U/t=14, β =2U/t=14, β =3U/t=14, β =4 FIG. S3. Conventional sign problem and sampling sign problem away from half filling: dependence on inverse temperature βt .System size L A = L = 8. as M (cid:48) = M/ ( τ DQMCAC × R av.max. ) ↔ M (exp.) . (S26)Here, τ DQMCAC is some measure of the autocorrelation of the Markov chain generated by the standard determinantalQMC algorithm.The ”sampling sign“ deteriorates exponentially with the probe area N sites ,A = L A . Yet, an extent of the probe area L A ∼ ξ ( T /t, U/t, L ), where ξ is the correlation length, is often sufficient for meaningful simulations. CODE VERIFICATION
For benchmarking purposes an irregular model instance of the Hubbard model on five sites (see inset in Fig. S4(b))is chosen H = − (cid:88) (cid:104) i,j, (cid:105) ,σ = ↑ , ↓ t σij (cid:16) c † i,σ c j,σ + h.c. (cid:17) + U (cid:88) i =1 n i, ↑ n i, ↓ − (cid:88) σ = ↑ , ↓ µ σ (cid:88) i =1 n i,σ , (S27)which breaks translational, point group and spin rotational symmetry. The hopping matrix is identical for both spinspecies and reads [ t σ ] i,j = t . . . . .
05 0 . . . .
05 0 1 . . . . . (S28). The onsite repulsion is U/t = 4 and the chemical potential for spin up and down is µ ↑ /t = 1 . µ ↓ /t = 1 . βt = 4 discretized into N τ = 256 Trotter time slices with ∆ τ t = βt/N τ =1 /
64. Fig. S4 compares the probabilities P ( s ) of all microstates s in the occupation number basis with resultsfrom exact diagonalization. An enlarged view of Fig. S4(a) is shown in Fig. S4(b). The occupation number state (cid:126)n = [ n , ↑ . . . n , ↑ ; n , ↓ . . . n , ↓ ] with n i,σ ≡ n α ∈ { , } is interpreted as a bitstring and is represented by the integer s = (cid:80) N sites − α =0 n α α . With the knowledge of the eigenenergies E i and eigenvectors | φ i (cid:105) of the Hamiltonian H theprobability of an occupation number state s is P ( s ) = dim( H ) (cid:88) i =1 e − βE i |(cid:104) s | ϕ i (cid:105)| . (S29)In order to demonstrate that the nested componentwise direct sampling method yields consistent microscopicstates also for much larger system sizes, Fig. S5 shows the probabilities of spin-only states on the eight sites aroundan isolated, mobile hole for a system of 10 ×
10 lattice sites at small doping. Only isolated holes, i.e. those surroundedexclusively by singly occupied sites, are considered (see Eq. (13) of the main text). With the labelling of sites arounda hole as shown in the right panel of Fig. S5, the spin environment is characterized by the state (cid:126)σ = [ σ σ . . . σ ]where σ p = 1 if there is a spin- ↑ at position p around the hole and σ p = 0 if it is a spin- ↓ . With this convention, thestate index in Fig. S5 is given as s = (cid:80) p =0 σ p p . Fig. S5 shows that there is a hierarchy of groups of states (rightpanel of Fig. S5). Furthermore, states related by symmetry appear with approximately the same probability, whichis not a built-in feature of the nested componentwise sampling algorithm and thus provides strong evidence that allrelevant occupation number states are sampled with their correct probabilities. FCS OF NUMBER OF DOUBLONS AND HOLES
Fig. S6 displays the doping dependence of the FCS of the number of doubly occupied sites N d = (cid:80) N sites ,A i =1 n i, ↑ n i, ↓ and holes N h = (cid:80) N sites ,A i =1 (1 − n i, ↑ ) (1 − n i, ↓ ) for U/t = 12 , βt = 4 on an L × L system with L = L A = 8. For smalldoping values, the FCS are well described by binomial distributions (dashed and dotted lines) which are generated by P ( s ) state s componentwise samplingED00.0020.0040.0060.0080.010.0120.0140.016 0 100 200 300 400 500 P ( s ) state s FIG. S4. Probabilities P ( s ) of all microstates s in the occupation number basis. Comparison with exact diagonalization (ED)for the parameters shown in the upper panel and for the hopping matrix of Eq. (S28). The total number of Fock samples is10 , generated by 56 independent Markov chains. P ( s ) spin-only state s µ /t=-2.0 etc.20 states (and their spin-inversion partners) etc.4 states (and spin-inversion partners) 4 states (and spin-inversion partners) etc.Néel environment FIG. S5. (a) Probabilities of spin-only states around an empty site as obtained from reweighted pseudo-snapshots in thereference frame of the hole. Parameters:
U/t = 12 , µ/t = − , βt = 2 . L × L square lattice with L = 8. For this systemsize, the average number of holes and doubles (cid:104) N h (cid:105) ≈ .
50 and (cid:104) N d (cid:105) ≈ .
31 so that the average number of excess holes is (cid:104) N h (cid:105) − (cid:104) N d (cid:105) ≈ .
19. With 56 independent Markov chains a total number of 70 × occupation number samples of theentire system was obtained. In the reference frame of the hole, those states around the hole position are filtered out which donot contain charge fluctuations (“spin-only states”). Per Hubbard-Stratonovich configuration, 200 direct sampling steps wereperformed during the Monte Carlo sweeps at different imaginary time slices. Trotter discretization ∆ τ t = 0 .
02. (b) Illustrationof the most important spin-only states around an isolated hole which are marked in panel (a). Dashed lines delineate differentantiferromagnetic domains. populating lattice sites independently with doublons (holes) with probability p d = (cid:104) n ↑ n ↓ (cid:105) ( p h = (cid:104) (1 − n ↑ )(1 − n ↓ (cid:105) ) (cid:105) ).according to the distribution B n,p ( x ) = (cid:32) nx (cid:33) p x (1 − p ) n − x with n = N sites ,A and p ∈ { p d , p h } . In the stronglydoped case (lower panel), the distribution of the number of holes can be modelled accurately by a “shifted binomial”distribution (dashed-dotted line), which is obtained by populating the lattice with holes arising from doublon-holefluctuations according to a binomial distribution with parameter p = p d (since the number of doublon-hole pairs is µ /t = -2.50
U/t=14L=6, β t = 2.5|r| (e) Koepsell et al. (2019)
FIG. S7. The three-point spin-charge correlation function C | d | =2 ( | r | ) of Eq. (12) computed on a 4 × × U/t = 14 , βt = 2 . , µ = − . , − . , . . . , − .
5. Theexperimental data of Ref. [19] are also shown for comparison. There, the system size is approximately 4 × × U/t ≈
14 and stated temperature T = 1 . J , whichcorresponds to T /t ≈ . βt = 2 .
5) for J = 4 t /U . In the quantum gas microscope experiment, there are on average1 . × µ/t = − .
0, on the other hand, the average number of holes (hole-doped) is (cid:104) N h (cid:105) = 0 .
72, which is small enoughfor measuring the spin-environment of an isolated dopant. (d) Randomly chosen snapshots for µ/t = − µ/t = − x , x=d,h P(N d )P(N h ) FIG. S8. FCS of the number of doubly occupied ( N d ) and empty sites ( N h ) for the same parameters as in Fig. 3 of the maintext, which are L = 10 , U/t = 14 , µ/t = − . , βt = 2 . ADDITIONAL DATA FOR THE POLARON PROBLEM
This section illustrates the significance of numerical simulations for exploring the parameter space. In agreementwith Ref. [27] for the t − J model and the results of Fig. 3 for the Hubbard model, one can identify the joint condition C | d | =2 ( | r | = 0) < C | d | =2 ( | r | = 1) > C | d | =2 ( | r | = 0) < C | d | =2 ( | r | ) = (cid:104) n ( h ) r S z r + r + d S z r + r − d (cid:105)(cid:104) n ( h ) r (cid:105) . (S30)as a function of chemical potential µ/t for L = 4 (c) and L = 6 (e), respectively, with the distribution of the number ofholes (b,g) and doubly occupied sites (a,f). Unlike in Eq. (12) of the main text, for the study of the hole environmentin Fig. S7 no projection operator restricting the sites around the hole to be singly occupied has been applied prior toevaluating the correlation function. The data points connected by a thick red line are taken from Fig. 4(c) of Ref. [19].From µ/t = − µ/t = − . C | d | =2 ( | r | ) changes markedly which leads to the conclusion that Fig. 4of Ref. [19] is consistent with the disappearance rather than the presence of magnetic polarons, which is due to therelatively high level of doping chosen in Ref. [19]. This picture is supported by visually comparing the two randomlyselected pseudo-snapshots for µ/t = − µ/t = −−