Sampling Electronic Fock States using Determinant Quantum Monte Carlo

For the Hubbard model, the determinant quantum Monte Carlo employs the Hubbard-Stratonovich decomposition to map the expectation of observables in an interacting system into a statistic average of measurements in an effective non-interacting system that couple to an auxiliary field [3, 4, 5, 6, 7]. This decomposition is expressed as

$\displaystyle Z=\mathrm{Tr}[\rho]=\sum_{x}\mathrm{Tr}[{\rho}_{x}],$

(1)

where $Z$ is the partition function, $x$ is the Hubbard-Stratonovich field, and ${\rho}_{x}=\mathcal{T}e^{-\int_{0}^{\beta}\sum_{ij}c^{\dagger}_{i}H_{ij}[x(\tau)]c_{j}d\tau}$ is the density matrix for the effective non-interacting system, associated with an imaginary-time-dependent auxiliary field $x(\tau)$. Here, $H[x(\tau)]$ represents the Hamiltonian for the effective non-interacting system.

Unlike the traditional DQMC, we can further project ${\rho}_{x}$ to the Fock-state basis $\ket{\eta}$ before evaluating the expectation values of observables. That is,

$\displaystyle Z=\sum_{x,\ket{\eta}}Z_{x,\ket{\eta}}=\sum_{x,\ket{\eta}}\matrixelement{\eta}{\rho_{x}}{\eta},$

(2)

where $\ket{\eta}$ is a binary vector that specifies a fermionic Fock state in the real-space representation. For example, $\ket{1010}$ represents the state $c_{1}^{\dagger}c_{3}^{\dagger}|0\rangle$. This projection mimics the snapshot sampling in the quantum gas microscope of quantum simulations [15]. With these projections, the expectation value of an observable $O$ is calculated by

$\displaystyle\langle{O}\rangle=Z^{-1}\sum_{x,\ket{\eta}}Z_{x,\ket{\eta}}\langle{O}\rangle_{x,\ket{\eta}},$

(3)

where $\langle{O}\rangle_{x,\ket{\eta}}=\matrixelement{\eta}{O\rho_{x}}{\eta}/\matrixelement{\eta}{\rho_{x}}{\eta}$. Here, $Z_{x,\ket{\eta}}$ is not positive semi-definite since the projection onto a random Fock state breaks the particle-hole symmetry [46]. Therefore sign re-weighting is needed by adjusting Eq. (3) into

$\displaystyle\langle{O}\rangle=\frac{\sum_{x,\ket{\eta}}\absolutevalue{Z_{x,\ket{\eta}}}{\rm sgn}\quantity(Z_{x,\ket{\eta}})\langle{O}\rangle_{x,\ket{\eta}}}{\sum_{x,\ket{\eta}}\absolutevalue{Z_{x,\ket{\eta}}}{\rm sgn}\quantity(Z_{x,\ket{\eta}})}.$

(4)

FDQMC utilizes $\absolutevalue{Z_{x,\ket{\eta}}}$ as the joint statistical weight for $x$ and $\eta$ to generate a significant number of $(x,\eta)$ pairs. The statistical average over these pairs allows an unbiased evaluation for the expectation value $\langle{O}\rangle$. We adopt the Metropolis-Hasting update for Markov-chain importance sampling. Specifically, for a proposed flip $(x,\eta)\rightarrow(x^{\prime},\eta^{\prime})$, the acceptance probability is calculated as $P_{\rm acc}=\min\{\absolutevalue{R_{\rm acc}},1\}$, where $R_{\rm acc}=Z_{x^{\prime},\ket{\eta^{\prime}}}/Z_{x,\ket{\eta}}$. As shown in Fig. 1, FDQMC updates introduce an additional dimension compared to traditional DQMC. It achieves importance sampling of Fock states with a signed uniform weight. Ideally, for large sample size, the physical distribution of Fock-state snapshots in the original system can be reproduced by partially canceling positive-weight samples with those carrying negative weights,

$\displaystyle\matrixelement{\eta}{\rho}{\eta}=\sum_{x}\absolutevalue{Z_{x,\ket{\eta}}}{\rm sgn}\quantity(Z_{x,\ket{\eta}}).$

(5)

In practice, insufficient sample size may result in a non-positive distribution. However, if the sign problem is not severe, these samples can accurately reproduce all high-order correlations in an unbiased manner through Eq. (4). A systematic investigation of the sign problem is presented in Supplementary Note 2. The update and measurement strategies are detailed in Methods.

Unlike traditional DQMC, direct pre-selection for any ensembles and physical constraints is easily achieved in FDQMC, through the Fock-state restrictive sampling (FRS). This capability is essential for accurately evaluating observables, particularly doping-sensitive observables, in finite systems where the grand-canonical and canonical ensembles differ significantly, and particle number fluctuations can introduce substantial noise. FRS scheme is analogous to the post-selection method used in analyzing quantum gas microscope experiments [33, 34, 35, 36, 37, 38, 39, 40], but is performed before the Monte Carlo update without generating unused samples. The grand-canonical ensemble is simulated by default if no restrictions are applied. In this work, we examine three different ensembles. The canonical ensemble with a fixed particle number can be simulated by randomly swapping a site-$i$ particle and a site-$j$ hole for each $\eta$-update (see Fig. 2a). Building upon this, the spin-selected ensemble further mandates a consistent total spin (in the z-direction) and is realized by restricting particle-hole swaps within each spin sector (see Fig. 2b). Finally, the non-doublon ensemble excludes double occupation on the top of the spin-selected ensemble, often used in quantum simulations. In the context of the Hubbard model, this ensemble serves as an extension of the $t$-$J$ model, including all higher-order spin-exchange processes. The non-doublon ensemble is achieved by swapping up-spins, down-spins, and holes individually (see Fig. 2c). As will be elaborated later, the non-doublon ensemble is more effective in signaling high-order correlations. Additionally, constraints on fixing on-site fermion number enable FDQMC to simulate spin-fermion interactions, such as Kondo coupling, through slave fermions (see Fig. 2d).

We apply FDQMC to the single-band Fermi-Hubbard model in a 2D square lattice, whose Hamiltonian is

$\displaystyle\mathcal{H}_{\rm H}=-t\sum_{\langle{ij}\rangle,s}(c^{\dagger}_{is}c_{js}+{\rm h.c.})+U\sum_{i}n_{i\uparrow}n_{i\downarrow}\,.$

(6)

Here, $t$ denotes the hopping between nearest neighbors, and $U>0$ represents the on-site repulsive interaction. In this section, we examine static two-point and higher-order correlations across various ensembles.

Leveraging the capability to select specific ensembles, the simulated results of FDQMC can be benchmarked against the quantum gas microscope experiments in Ref. 47. Figure 3a shows the distributions of staggered magnetization at half-filling, defined as

$\displaystyle M^{\rm(st)}_{z}=\frac{2}{N}\sum_{{\bf\it i}}(-1)^{i_{x}+i_{y}}S^{\rm z}_{{\bf\it i}}.$

(7)

The simulated distributions are statistically derived from $\sim 4\times 10^{6}$ Fock-state samples by FDQMC, while the experimental results are obtained using $\sim$ 250 snapshots via quantum gas microscope [47]. Here, $S_{{\bf\it i}}^{\rm z}=({n}_{{\bf\it i}\uparrow}-{n}_{{\bf\it i}\downarrow})/2$ represents the $z$-component of the spin at site ${\bf\it i}$, $N$ denotes the total number of sites, and the Hubbard interaction is set to $U=7.2\,t$ in align with experiments. We select an $8\times 10$ periodic square lattice and employ the canonical ensemble to compare with the experimental results measured within a circular central region containing approximately 80 sites. At high temperatures, the simulated $M^{\rm(st)}_{z}$ closely aligns with the distribution obtained from experiments, appearing as a Gaussian envelope center at zero with spin symmetry. As the temperature falls below the spin-exchange energy $J\sim 0.55\,t$, the distribution noticeably broadens due to quantum fluctuations driven by antiferromagnetic (AFM) correlations [48]. These fluctuations are reflected numerically by more distinct Fock-state configurations within the same sample volume. Therefore, the experimental distribution starts to deviate from a smooth and symmetric distribution, constrained by its 250 snapshots. In contrast, the distribution obtained by FDQMC can reach exact solutions with minimal statistical error, benefiting from the extensive sample volume ($\sim 4\times 10^{6}$). For a detailed comparisons, please refer to Supplementary Note 1.

When a single hole is introduced into the AFM background, it can influence magnetization [49, 50, 51, 52]. Here, we examine a single-hole-doped Hubbard model with $U=8t$ on a $6\times 6$ periodic square lattice. Given the odd number of electrons with the presence of a single hole, we opt for the spin-selected and non-doublon ensembles and set the total spin to be $1/2$. With this selected orientation, the magnon dressing of the spin polaron can be visualized through the (connected) spin-hole correlation [53, 36]

$\displaystyle C^{\rm(con)}_{\rm h,s}({\bf\it d})=2\langle{{n}_{{\bf\it 0}}^{\rm h}{S}_{{\bf\it d}}^{\rm z}}\rangle/\langle{{n}_{{\bf\it 0}}^{\rm h}}\rangle-2\sum_{{\bf\it r}}\langle{S_{{\bf\it r}}^{\rm z}}\rangle/N,$

(8)

where ${n}_{{\bf\it 0}}^{h}=(1-{n}_{{\bf\it 0}\uparrow})(1-{n}_{{\bf\it 0}\downarrow})$ is the hole operator at the origin. As shown in Fig. 3b, the high-temperature system is disordered, with no magnetization except that the additional spin moment accumulates near the nearest neighbor of the hole. With the decrease of temperature, the AFM order starts to develop, and the motion of the hole is dressed by the disturbance of the AFM correlations, resulting in the checkerboard distribution of $C^{\rm(con)}_{\rm h,s}({\bf\it d})$ throughout the system.

The magnetization is further influenced by the mobility of the hole. We delve into this effect by adding a pinning potential $V$ at the origin site to tune the mobility of the dopant. This leads to the modified Hamiltonian

$\displaystyle H_{V}=H_{\rm H}+V(n_{{\bf\it 0}\uparrow}+n_{{\bf\it 0}\downarrow})\,.$

(9)

The pining potential can be realized experimentally using an optical tweezer in an optical lattice [54, 55, 36]. While similar to the $V=0$ results at high temperatures, the staggered magnetization starts to develop at a higher temperature with a strong pinning potential. In the lowest temperature ($T=0.2t$), the magnetization resembles a Néel state with a missing down-spin at the origin. Given the reduced hole’s mobility with a pronounced $V$, it serves as a geometric defect at low temperatures.

A significant advance of quantum simulations is the analysis of multi-point high-order correlations, which extends beyond the capabilities of conventional spectroscopic measurements in solid-state materials. In the context of the Hubbard model, connected high-order correlations have been utilized to uncover the hidden orders, polaronic wavefunctions, and entanglement [33, 34, 35, 36, 37, 38, 39, 40, 45].

Following the formalism in Ref. [56], we examine the property of a single-hole-doped Hubbard model using FDQMC on an $6\times 6$ square lattice through the analysis of the third- and fifth-order correlations. The connected part of the third-order correlation is defined as [36, 56]

$\displaystyle B^{\rm(con)}({\bf\it r},{\bf\it r}^{\prime},{\bf\it r}^{\prime\prime})=4\langle{n_{{\bf\it r}}^{h}S_{{\bf\it r}^{\prime}}^{\rm z}S_{{\bf\it r}^{\prime\prime}}^{\rm z}}\rangle/\langle{n_{{\bf\it r}}^{h}}\rangle-4\langle{S_{{\bf\it r}^{\prime}}^{\rm z}S_{{\bf\it r}^{\prime\prime}}^{\rm z}}\rangle,$

(10)

highlighting the impact of a hole (at site ${\bf\it r}$) on the spin-spin correlation $\langle{S_{{\bf\it r}^{\prime}}^{\rm z}S_{{\bf\it r}^{\prime\prime}}^{\rm z}}\rangle$ (see the inset of Fig. 4a). The temperature dependence of $B^{\rm(con)}({\bf\it r},{\bf\it r}^{\prime},{\bf\it r}^{\prime\prime})$, for the diagonal spin correlations $({\bf\it r}^{\prime},{\bf\it r}^{\prime\prime})=(\hat{{\bf\it x}},\hat{{\bf\it y}})$ with respect to the hole at ${\bf\it r}=0$, is shown in Fig. 4a. Its disconnected part, i.e. $4\langle{S_{{\bf\it r}^{\prime}}^{\rm z}S_{{\bf\it r}^{\prime\prime}}^{\rm z}}\rangle$, is positive in the AFM background, and $B^{\rm(con)}({\bf\it r},{\bf\it r}^{\prime},{\bf\it r}^{\prime\prime})$ indicates the enhancement or diminishment of this correlation near a hole. Without the pinning potential, $B^{\rm(con)}$ consistently exhibits a negative value, reflecting the disturbance on the AFM spin correlations by the hole’s motion. This connected correlation serves as a fingerprint for a spin polaron. Nonetheless, with the introduction of a strong pinning potential $V=5t$, the sign of $B^{\rm(con)}$ transitions to positive at sufficiently low temperature ($T<0.25\,t$). Such a flip signals an “anti-screening” effect of the hole at low $T$, strengthening the spin correlations near the hole. This effect is attributed to the reduction of spin fluctuations with fewer neighboring sites, when the hole is immobile and becomes effectively a geometric defect [56]. Such a transition from polaronic screening to anti-screening only occurs when $V$ is adequately large to surpass the kinetic energy of the hole (see Fig. 4b). It is important to note that, although the results from different ensembles vary quantitatively, the critical temperature and pinning potential for the transition remain unaffected by the choice of ensemble.

Another high-order correlation depicting the single-hole dynamics is the fifth-order hole-spin-ring correlation, defined as [57, 56]

$\displaystyle C_{\rm\diamondsuit}=2^{4}\langle{n_{{\bf\it 0}}^{\rm h}{S}_{{\bf\it r}+\hat{{\bf\it x}}}^{\rm z}{S}_{{\bf\it r}+\hat{{\bf\it y}}}^{\rm z}{S}_{{\bf\it r}-\hat{{\bf\it x}}}^{\rm z}{S}_{{\bf\it r}-\hat{{\bf\it y}}}^{\rm z}}\rangle/\langle{n^{\rm h}_{{\bf\it 0}}}\rangle\,,$

(11)

where the hole at the origin is encircled by a spin ring. Similar to the third-order correlation, $C_{\rm\diamondsuit}$ remains negative for a mobile hole across all temperatures (see Fig. 4c), reflecting the string excitation caused by the formation of spin polaron [57]. This negativity stems from spin correlations of the AFM background, hence intensifying at lower temperatures. With a pinning potential in place, this fifth-order correlation aligns with the $V=0$ scenario at high temperatures due to the lack of AFM order. Yet, as the temperature drops significantly below $J$, the emergence of AFM correlation and the suppression of quantum fluctuations by the geometric defect help the development of a pronounced spin-ring correlation surrounding the immobile hole. This effect is evidenced by the substantial positive values of $C_{\rm\diamondsuit}$ at low temperatures ($T<0.3\,t$). The transition from negative to positive also occurs around $V_{c}\sim 2t$ (see Fig. 4d), indicating that the coincident underlying anti-screening physics as observed in the third-order correlation $B^{\rm(con)}$.

Expanding our analysis to higher dopings, the fermion-sign problem becomes more pronounced at low temperatures, restricting our simulations to relatively high temperatures. Here, we choose $U=8t$ and $T=0.5t$ on the periodic $8\times 8$ lattice, matching the conditions of cold-atom experiments in Ref. 40, where $U=7.4(8)t$ and $T=0.52(5)t$. An even number of doped holes are used for the FDQMC simulation to ensure that the total spin of the spin-selected and non-doublon ensembles is zero. This setting reflects the experimental reality and simplifies the simulation, as all low-order correlations involving an odd number of spin operators are nullified.

Extending the simulation of the connected third-order correlation $B^{\rm(con)}({\bf\it r},{\bf\it r}+{\bf\it r}^{\prime},{\bf\it r}+{\bf\it r}^{\prime\prime})$ into finite doping, we focus on two key distances indicative of spin polaron wavefunction: the nearest-neighbor-spin $B^{\rm(con)}(\mathrm{NN\!-\!spin})=B^{\rm(con)}({\bf\it r},{\bf\it r}+\hat{{\bf\it y}},{\bf\it r}+\hat{{\bf\it x}}+\hat{{\bf\it y}})$ and the nearest-diagonal-spin $B^{\rm(con)}(\mathrm{ND\!-\!spin})=B^{\rm(con)}({\bf\it r},{\bf\it r}+\hat{{\bf\it x}},{\bf\it r}+\hat{{\bf\it x}}+\hat{{\bf\it y}})$ (illustrated in the insets of Figs.  5a and b). Without the pinning potential, the system is translational symmetric, and the choice of ${\bf\it r}$ is irrelevant. As shown in Fig. 5a, $B^{\rm(con)}(\mathrm{NN\!-\!spin})$ obtained from all ensembles exhibits a rapidly decrease from a significantly positive value, transitioning to negative at $\sim 18\%$ doping — a change potentially linked to the temperature-independent quasi-particle interruption observed in ARPES studies of cuprates [58]. This sign flip reflects the breakdown of spin polaron with increasing doping. When comparing results from the canonical ensemble with cold atom experiments in Ref.[40], a consistent agreement is observed throughout all dopings. This consistency further validates the efficacy of FDQMC samples in mirroring quantum simulation snapshots. A similar agreement is observed for $B^{\rm(con)}(\mathrm{ND\!-\!spin})$, as presented in Fig. 5b. Results from the non-doublon ensemble, however, deviate from the canonical ensemble and experimental results in the low doping regime, attributed to the exclusion of doublon-hole fluctuations.

The non-monotonic doping dependence of $B^{\rm(con)}(\mathrm{NN\!-\!spin})$ indicates that the doped Hubbard model maintains a strongly correlated state beyond the breakdown of the spin polaron, which has also been suggested by the persistent spin fluctuations observed in cuprates [59, 60, 61]. Particularly, potential interactions between holes may be mediated by the overlap of two spin-polarons [62]. Therefore, we examine the fourth-order correlations involving two holes [40]

		$\displaystyle\quad D^{\rm(con)}({\bf\it r}_{1},{\bf\it r}_{2},{\bf\it r}_{3},{\bf\it r}_{4})$		(12)
		$\displaystyle=\frac{1}{\langle{n^{\rm h}_{{\bf\it r}_{1}}\mkern-2.0mun^{\rm h}_{{\bf\it r}_{2}}}\rangle}\mkern-2.0mu\Big{[}\langle{n^{\rm h}_{{\bf\it r}_{1}}\mkern-2.0mun^{\rm h}_{{\bf\it r}_{2}}\mkern-2.0muS^{\rm z}_{{\bf\it r}_{3}}\mkern-2.0muS^{\rm z}_{{\bf\it r}_{4}}}\rangle-\langle{n^{\rm h}_{{\bf\it r}_{1}}}\rangle\langle{n^{\rm h}_{{\bf\it r}_{2}}S^{\rm z}_{{\bf\it r}_{3}}S^{\rm z}_{{\bf\it r}_{4}}}\rangle-\langle{n^{\rm h}_{{\bf\it r}_{2}}}\rangle$
		$\displaystyle\cdot\mkern-2.0mu\langle{n^{\rm h}_{{\bf\it r}_{1}}\mkern-2.0muS^{\rm z}_{{\bf\it r}_{3}}\mkern-2.0muS^{\rm z}_{{\bf\it r}_{4}}}\rangle\mkern-2.0mu-\mkern-2.0mu\langle{\mkern-1.0mun^{\rm h}_{{\bf\it r}_{1}}\mkern-2.0mun^{\rm h}_{{\bf\it r}_{2}}\mkern-1.0mu}\rangle\langle{S^{\rm z}_{{\bf\it r}_{3}}\mkern-2.0muS^{\rm z}_{{\bf\it r}_{4}}}\rangle\mkern-2.0mu+\mkern-2.0mu2\langle{\mkern-1.0mun^{\rm h}_{{\bf\it r}_{1}}\mkern-1.0mu}\rangle\langle{\mkern-1.0mun^{\rm h}_{{\bf\it r}_{2}}\mkern-1.0mu}\rangle\mkern-2.0mu\langle{S^{\rm z}_{{\bf\it r}_{3}}\mkern-2.0muS^{\rm z}_{{\bf\it r}_{4}}}\rangle\mkern-2.0mu\Big{]}\,,$

This connected correlation quantifies the net effect of a pair of holes on adjacent spin correlations, compared to separated ones. When the examined two spins and two holes form a plaquette, i.e. NN-spin and ND-spin as shown in Figs. 5c and 5d, $D^{\rm(con)}$ manifests significant values across a wide range of doping (up to $\sim 60$%). In particular, $D^{\rm(con)}(\mathrm{NN\!-\!spin})$ is consistently negative, reflecting a tendency for spin polarons to share spin defects. This correlation monotonically decreases as the system is doped away from the AFM phase. At the same time, $D^{\rm(con)}(\mathrm{ND\!-\!spin})$ becomes significantly negative only near 20% doping, where the spin polarons break down. Both the NN- and the ND- $D^{\rm(con)}$ suggest the preference of spin-singlet around the closest proximity of the hole pair, consistent with experimental findings [40].

When compared against various theories, it has been found that analytical wavefunctions ansatzes fail to capture the doping evolution observed in experiments [40]. Numerical simulations using finite-temperature ED have successfully reproduced third-order correlations $B^{\rm(con)}$ and qualitatively traced the trend of the fourth-order correlations $D^{\rm(con)}$. However, discrepancies of a factor of 2 to 2.5 are present in the ED simulations, largely due to the finite-size effects. Using the FDQMC method at the same size and temperature as the experiments, we manage to closely match experimental results for $D^{\rm(con)}(\mathrm{NN\!-\!spin})$ and significantly reduce the discrepancy of $D^{\rm(con)}(\mathrm{ND\!-\!spin})$ to around $50$%. Since the DQMC simulations are unbiased at this temperature and system size, the remaining mismatch likely stems from the inhomogeneity and the uncertainty of model parameters in experiments or the distinction in boundary conditions. Upon comparing across different ensembles, we find that high-order correlations evaluated in the canonical and the spin-selected ensemble are similar, whereas the non-doublon ensemble leads to more pronounced correlations below 40% doping. This indicates that the non-doublon ensemble is more suitable for elucidating genuine hole-spin correlations and uncovering their entanglement, due to the exclusion of irrelevant Fock states that involve double occupation.

While simulating dynamical correlations presents challenges with analog quantum simulators [63], the Fock state sampling can be extended to the analysis of unequal-time correlations and allows for the emulation of spectroscopies similar to traditional DQMC. For a specific configuration of $(x,\eta)$, the unequal-time Green’s function is calculated as (assuming $\tau_{1}\geq\tau_{2}$)

	$\displaystyle\langle{c_{i}(\tau_{1})c_{j}^{\dagger}(\tau_{2})}\rangle_{x,\ket{\eta}}$	$\displaystyle=$	$\displaystyle\quantity[B_{x}(\tau_{1},\tau_{2})G_{x,\ket{\eta}}(\tau_{2},\tau_{2})]_{ij},$		(13)
	$\displaystyle\langle{c_{j}^{\dagger}(\tau_{1})c_{i}(\tau_{2})}\rangle_{x,\ket{\eta}}$	$\displaystyle=$	$\displaystyle\quantity[\big{(}I-G_{x,\ket{\eta}}(\tau_{2},\tau_{2})\big{)}B^{-1}_{x}(\tau_{1},\tau_{2})]_{ij}$

using the numerically stable method introduced in Ref. 64. Other dynamical observables are derived from Green’s functions using Wick’s theorem. Here, we discuss two representative examples: The single-particle spectrum $A({\bf\it k},\omega)$ measures the evolution of an individual electron, and the dynamical spin structure factor $S({\bf\it q},\omega)$ measures the propagation of a spin excitation (see definitions in Methods). Both $A({\bf\it k},\omega)$ and $S({\bf\it q},\omega)$ are analytically continued using the maximum entropy method [65]. For the sake of visualization, we normalize $S({\bf\it q},\omega)$ for each momentum, denoted as $\tilde{S}({\bf\it q},\omega)$, while the original data are shown in Supplementary Note 3.

Building on above analyses, we first study the half-filled Hubbard, utilizing an $8\times 8$ cluster within the canonical ensemble. As shown in Fig. 6a, a Mott gap of $\sim 4t$ is evident in $A({\bf\it k},\omega)$. The Hubbard interaction also leads to the formation of 2D AFM order, as manifested by the magnon dispersion in $\tilde{S}({\bf\it q},\omega)$ (see Fig. 6b). Different from small-cluster ED simulations, the $\tilde{S}({\bf\it q},\omega)$ evaluated from FDQMC correctly portrays the splitting between the nodal $(\pi/2,\pi/2)$ and the anti-nodal $(0,\pi)$ magnons, a discrepancy stemming from higher-order spin-exchange processes at large systems.

To demonstrate the advantage of the FRS scheme in FDQMC, we turn our attention to the Kondo lattice model, described by the Hamiltonian

$$\mathcal{H}_{\rm K}\mkern-2.0mu=\mkern-2.0mu-t\mkern-4.0mu\sum_{\langle{ij}\rangle,s}\mkern-2.0mu(c^{\dagger}_{is}c_{js}+{\rm h.c.})+\frac{J_{\rm K}}{2}\mkern-4.0mu\sum_{i,s,s^{\prime}}\mkern-2.0muc^{\dagger}_{is}{\bf\it\sigma}_{ss^{\prime}}c_{is^{\prime}}\mkern-2.0mu\dotproduct\mkern-2.0mu{\bf\it S}^{f}_{i}\,.$$

(14)

Here, the spin of an itinerant electron couples to a localized spin-1/2 moment ${\bf\it S}^{f}_{i}$ with strength $J_{\rm K}$. By decomposing ${\bf\it S}^{f}_{i}$ using the slave-fermion representation $(1/2)\sum_{s,s^{\prime}}f^{\dagger}_{is}{\bf\it\sigma}_{ss^{\prime}}f_{is^{\prime}}$, the Kondo coupling in Eq. (14) is equivalent to an effective local, fermionic interaction

$\displaystyle-\frac{J_{\rm K}}{4}\quantity(\sum_{s}c^{\dagger}_{is}f_{is}+{\rm h.c.})^{2}+{\rm const.}\,,$

(15)

which is feasible for simulation using DQMC algorithms [66]. However, the slave-fermion decomposition requires the constraint $n^{f}_{i\uparrow}+n^{f}_{i\downarrow}=1$, a condition not inherently met by traditional DQMC at finite temperatures. Common solutions to this dilemma include switching to the Anderson model with $f$-site repulsion $U_{f}\propto 1/{J_{\rm K}}$, or introducing a strong $U_{f}$ to suppress both double and zero occupancies [66]. These approaches, however, induce approximations and potential bias. The strong artificial repulsion also affects the numerical stability, resulting in severe sign problems when doping. In contrast, the FRS scheme of FDQMC strictly adheres to the slave-fermion constraint, thereby serving as an exact finite-temperature solver for Kondo-type spin-fermion interactions.

Figures 6c-f present FDQMC-simulated spectral results of a half-filled Kondo-lattice model in a $8\times 8$ cluster. We choose the Kondo coupling ${J_{\rm K}}=2t$ and consider itinerant electrons in the canonical ensemble. The Kondo resonance occurs at a low temperature $T=0.2t$, resulting in the hybridization gap and the heavy-fermion dispersion in $A({\bf\it k},\omega)$ (see Fig. 6c). At the same time, $\tilde{S}({\bf\it q},\omega)$ exhibits pronounced magnon dispersion at $T=0.2\,t$ (see Fig. 6d). The excitation energy at ${\bf\it q}=(\pi,\pi)$ remains finite, revealing the spin-gapped phase for $J>J_{c}\approx 1.45t$ [66]. At high temperatures, e.g. $T=t$ presented in Figs. 6e and f, the hybridization gap closes, resulting in diminished magnon excitations. The decoupling between itinerant electrons and local spins reflects the asymptotic freedom of the Kondo lattice ($J_{\rm K}>0$) in the high-temperature limit.

For each update epoch, we alternately perform updates for the auxiliary field $x$ and the Fock state $\ket{\eta}$. The update algorithms leverage two different formulations of $Z_{x,\ket{\eta}}$. On one side, $Z_{x,\ket{\eta}}$ is a determinant

$\displaystyle Z_{x,\ket{\eta}}=\det\matrixquantity[P_{\ket{\eta}}^{\mathrm{T}}B_{x}P_{\ket{\eta}}]\,.$

(16)

Here, $\quantity[P_{\ket{\eta}}]_{ij}=\delta_{i,p_{j}}$ is the projection matrix of a Fock state $\ket{\eta}$, with $p_{j}$ denoting the site of the $j$-th particle. $B_{x}$ represents the single-particle propagator $\mathcal{T}e^{-\int_{0}^{\beta}H[x(\tau)]\differential{\tau}}$, associated with $x$. Eq. (16) shares the same mathematical structure used in the zero-temperature projective quantum Monte Carlo (PQMC) [3, 5, 6, 7]. Hence, the $x$-direction update follows the same strategy as PQMC.

At the same time, $Z_{x,\ket{\eta}}$ is proportional to the multi-point correlation under the auxiliary field configuration $x$ in the grand-canonical ensemble,

$\displaystyle Z_{x,\ket{\eta}}=\tr[\rho_{x}]\langle{\prod_{i=1}^{N_{m}}n_{i}(\eta)}\rangle_{x}\,.$

(17)

In this formula, $n_{i}(\eta)$ denotes the electron (hole) density at site $i$, if the site is (is not) occupied in the Fock state $\ket{\eta}$ and $N_{m}$ represents the system size. To facilitate the rank-1 update for $\ket{\eta}$, we construct an auxiliary matrix

where $\eta_{i}$ denotes electron density of $\ket{\eta}$ at site $i$ and $G^{(\rm gr)}_{x}=(I+B_{x})^{-1}$ is the equal-time Green’s function under auxiliary field $x$ in a grand-canonical ensemble. The numerically stable inversion of Eq. (Update strategy) is presented in Supplementary Note 5. The acceptance ratio for a proposed a $\eta$-flip at site $i$, namely $\eta^{\prime}_{j}=\eta_{j}+(-1)^{\eta_{i}}\delta_{ij}$, is

$\displaystyle R_{\rm acc}=-1-(-1)^{\eta_{i}}\quantity[M_{\ket{\eta}}]_{ii}.$

(18)

Upon acceptance of the flip, then $M_{\eta}$ is updated as

$\displaystyle\quantity[M_{\ket{\eta^{\prime}}}]_{jk}=\quantity[M_{\ket{\eta}}]_{jk}+\frac{(-1)^{\eta_{i}}}{R_{\rm acc}}\quantity[M_{\ket{\eta}}]_{ji}\quantity[M_{\ket{\eta}}]_{ik}.$

(19)

Each epoch of the Monte Carlo updates consists of first updating $x$ across all spatial and temporal sites, followed by proposing $\eta$-flips for each spatial site. The derivation of Eq. (Update strategy), (18), and (19) is in Supplementary Note 6.

The computational complexities for constructing the initial $M_{\ket{\eta}}$ and for performing the iterations of $\eta$-flips both scale as $\order{N_{m}^{3}}$. With the Sherman-Morrison fast update [6], the complexity of updating the auxiliary field $x$ is $\order{\beta N_{m}^{3}}$. Due to the imaginary time dimension, the cost for the latter overwhelms that of the $\eta$-upates. Therefore, the computational cost for FDQMC is only marginally more than that of the standard DQMC.

In the canonical ensemble, where the particle number is fixed, the Fock-state is updated by randomly swapping a particle at site $i$ with a hole at site $j$, as shown in Fig. 2a. The acceptance ratio for such a swap is given as

$\displaystyle\begin{aligned} &R^{\rm swap}_{\rm acc}=1+(-1)^{\eta_{i}}\quantity[M_{\ket{\eta}}]_{ii}+(-1)^{\eta_{j}}\quantity[M_{\ket{\eta}}]_{jj}\\ &+\mkern-2.0mu(-1)^{\eta_{i}+\eta_{j}}\quantity(\quantity[M_{\ket{\eta}}]_{ii}\mkern-1.0mu\quantity[M_{\ket{\eta}}]_{jj}\mkern-3.0mu-\mkern-3.0mu\quantity[M_{\ket{\eta}}]_{ij}\mkern-1.0mu\quantity[M_{\ket{\eta}}]_{ji}).\end{aligned}$

(20)

Upon acceptance of the swap, $\ket{\eta}$ and $M_{\ket{\eta}}$ are updated by successively imposing single-site $\eta$-flips [i.e. Eq. (19)] at the sites of $i$ and $j$.

The complexity of widely used canonical DQMC algorithms [67, 68, 69], which derive canonical ensemble properties by projecting from the grand-canonical ensemble using a Fourier projector, scales as $\order{\beta N_{m}^{4}}$ for observables beyond two-point correlations. In contrast, FDQMC avoids this additional overhead by directly sampling the canonical ensemble, keeping the complexity at $\order{\beta N_{m}^{3}}$. Meanwhile, truncation algorithms are applicable to FDQMC to further reduce the complexity of simulating dilute fermionic systems [70].

The expectation value of general observable $O$ is evaluated by decomposing into the Green’s function using Wick’s theorem. The equal-time Green’s function $[G_{x,\ket{\eta}}]_{ij}(\tau,\tau)=\langle{c_{i}(\tau)c_{j}^{\dagger}(\tau)}\rangle_{x,\ket{\eta}}$ for a specific $(x,\ket{\eta})$ is

$$G_{x,\ket{\eta}}(\tau,\mkern-2.0mu\tau)\mkern-2.0mu=\mkern-2.0muI-B_{x}\mkern-2.0mu(\tau,\mkern-2.0mu0)P_{\ket{\eta}}\mkern-4.0mu\quantity[\mkern-2.0muP_{\ket{\eta}}^{\mathrm{T}}B_{x}(\beta,\mkern-2.0mu0)P_{\ket{\eta}}\mkern-2.0mu]^{-1}\mkern-6.0muP_{\ket{\eta}}^{\mathrm{T}}B_{x}(\beta,\mkern-2.0mu\tau),$$

(21)

where the single-particle propagator from $\tau_{1}$ to $\tau_{2}$ is $B_{x}(\tau_{2},\tau_{1})=\mathcal{T}e^{-\int_{\tau_{1}}^{\tau_{2}}H[x(\tau)]\differential{\tau}}$. Note that $G_{x,\ket{\eta}}$ differs from the grand-canonical $G^{\rm(gr)}_{x}$ with unspecified $\ket{\eta}$. For observables $O=\sum_{\ket{\eta}}w_{\ket{\eta}}\outerproduct{\eta}{\eta}$ that are diagonal in real space, e.g., the charge and spin correlations, $\langle{O}\rangle_{x,\ket{\eta}}$ simplifies to $w_{\ket{\eta}}$. This simplification significantly reduces the computational cost for higher-order correlations by directly evaluating the sampled Fock states.

The single-particle spectrum is computed via $A({\bf\it k},\omega)=A_{+}({\bf\it k},\omega)+A_{-}({\bf\it k},\omega)$, where the particle-addition spectrum $A_{+}({\bf\it k},\omega)$ and the particle-removal spectrum $A_{-}({\bf\it k},\omega)$ are extracted from

	$\displaystyle\langle{c_{{\bf\it k}s}^{\dagger}(\tau)c_{{\bf\it k}s}(0)}\rangle$	$\displaystyle=\int A_{+}({\bf\it k},\omega)\mathrm{e}^{-\tau\omega}\differential{\omega},$		(22)
	$\displaystyle\langle{c_{{\bf\it k}s}(\tau)c^{\dagger}_{{\bf\it k}s}(0)}\rangle$	$\displaystyle=\int A_{-}({\bf\it k},\omega)\mathrm{e}^{\tau\omega}\differential{\omega},$		(23)

separately, with $c_{{\bf\it k}s}^{\dagger}=(1/\sqrt{N})\sum_{{\bf\it r}}c_{{\bf\it r}s}^{\dagger}\mathrm{e}^{\mathrm{i}{\bf\it k}\dotproduct{\bf\it r}}$. The dynamic spin structure factor $S({\bf\it q},\omega)$ is obtained through

$\displaystyle\langle{S_{{\bf\it q}}^{\rm z}(\tau)S_{{\bf\it q}}^{\rm z}(0)}\rangle=\int S({\bf\it q},\omega)\mathrm{e}^{-\tau\omega}\differential{\omega}$

(24)

with $S_{{\bf\it q}}^{\rm z}=(1/\sqrt{N})\sum_{{\bf\it r}}S_{{\bf\it r}}^{\rm z}\mathrm{e}^{\mathrm{i}{\bf\it q}\dotproduct{\bf\it r}}$. Both spectra in real frequency are calculated using the analytic continuation through the maximum entropy method [65].