Directed-Loop Quantum Monte Carlo Method for Retarded Interactions

TL;DR

This paper extends the directed-loop quantum Monte Carlo method to efficiently handle retarded interactions in fermion-boson and spin-boson models, enabling faster simulations of complex quantum phase transitions.

Contribution

The authors develop an exact quantum Monte Carlo algorithm that incorporates retarded interactions using path integrals, improving efficiency over previous methods.

Findings

Overcomes long autocorrelation times in simulations.

Achieves exponential speedup in sampling.

Enables study of large chains up to 1282 sites.

Abstract

The directed-loop quantum Monte Carlo method is generalized to the case of retarded interactions. Using the path integral, fermion-boson or spin-boson models are mapped to actions with retarded interactions by analytically integrating out the bosons. This yields an exact algorithm that combines the highly-efficient loop updates available in the stochastic series expansion representation with the advantages of avoiding a direct sampling of the bosons. The application to electron-phonon models reveals that the method overcomes the previously detrimental issues of long autocorrelation times and exponentially decreasing acceptance rates. For example, the resulting dramatic speedup allows us to investigate the Peierls quantum phase transition on chains of up to $1282$ sites.

Figures

Tables

Table 1. Table 1: Vertex weights W v 1 , v 2 subscript 𝑊 subscript 𝑣 1 subscript 𝑣 2 W_{v_{1},v_{2}} for the spinless Holstein model for all possible combinations of vertex types v 1 subscript 𝑣 1 v_{1} and v 2 subscript 𝑣 2 v_{2} .

	1	2	3	4	5	6
1	$\lambda t\left(C+\frac{1}{2}\right)$	$\lambda tC$	$\lambda tC$	$\lambda t\left(C-\frac{1}{2}\right)$	$t/2$	$t/2$
2	$\lambda tC$	$\lambda t\left(C+\frac{1}{2}\right)$	$\lambda t\left(C-\frac{1}{2}\right)$	$\lambda tC$	$t/2$	$t/2$
3	$\lambda tC$	$\lambda t\left(C-\frac{1}{2}\right)$	$\lambda t\left(C+\frac{1}{2}\right)$	$\lambda tC$	$t/2$	$t/2$
4	$\lambda t\left(C-\frac{1}{2}\right)$	$\lambda tC$	$\lambda tC$	$\lambda t\left(C+\frac{1}{2}\right)$	$t/2$	$t/2$
5	$t/2$	$t/2$	$t/2$	$t/2$	$0$	$0$
6	$t/2$	$t/2$	$t/2$	$t/2$	$0$	$0$

Equations

$$Z=\sum _{\alpha }\sum _{n=0}^{\infty }\frac{{\beta }^n}{n!}\sum _{S_n}⟨\alpha |\prod _{p=1}^n{\hat{H}}_{a_p,b_p}|\alpha ⟩.$$

$$\hat{H}=-t\sum _i{\hat{B}}_{i,i+1}+{\omega }_0\sum _i{\hat{a}}_i^{†}{\hat{a}}_i^{}+\gamma \sum _i{\hat{\rho }}_i^{}({\hat{a}}_i^{†}+{\hat{a}}_i^{})$$

$${\mathcal{S}}_{\text{ret}}=-2\lambda t\iint d{\tau }_{1}d{\tau }_2\sum _i{\rho }_i({\tau }_1)P({\tau }_1-{\tau }_2){\rho }_i({\tau }_2).$$

$${\mathcal{S}}_1=-\sum _{\nu }{w}_{\nu }{h}_{\nu }.$$

$$Z=\sum _{n=0}^{\infty }\sum _{C_n}\frac{Z_0}{n!}{w}_{{\nu }_1}\dots w_{{\nu }_n}{⟨h_{{\nu }_1}\dots h_{{\nu }_n}⟩}_0,$$

$$\sum _{S_n}{Z}_0{⟨h_{{\nu }_1}\dots h_{{\nu }_n}⟩}_0=\sum _{S_n}\sum _{\alpha }⟨\alpha |\prod _p{\hat{H}}_{a_p,b_p}|\alpha ⟩.$$

$${\mathcal{S}}_1=-\iint d{\tau }_{1}d{\tau }_{2}P({\tau }_1-{\tau }_2)\sum _{a_1,a_2,b}{h}_{a_1{a}_2,b}({\tau }_1,{\tau }_2).$$

$$h_{10,b}({\tau }_1,{\tau }_2)$$

$$h_{01,b}({\tau }_1,{\tau }_2)$$

$$h_{22,b}({\tau }_1,{\tau }_2)$$

$$b_1+a+b=W_{1,v_2},$$

$$a+b_2+c=W_{2,v_2},$$

$$b+c+b_3=W_{5,v_2}.$$

$$a$$

$$b$$

$$c$$

$$a$$

$$b$$

$$c$$

Full text

Introduction.—In the absence of general exact solutions for strongly correlated quantum systems, the development of efficient numerical methods is a central objective. For the one-dimensional (1D) case, the density-matrix renormalization group (DMRG) method [1; 2] has become the standard. However, it is much less efficient for higher dimensions, finite temperatures, or long-range interactions, so that quantum Monte Carlo (QMC) methods are often advantageous. The latter yield high-precision results for rather general 1D fermion and spin Hamiltonians. In particular, the cost for QMC simulations of path integrals in the stochastic series expansion (SSE) representation [3] scales linearly with system size $L$ and inverse temperature $\beta=1/k_{\text{B}}T$. Autocorrelation times are short due to the use of cluster updates (operator loops or directed loops) [4; 5; 6]. While usually restricted to 1D fermionic models by the sign problem, such QMC methods were successfully applied in higher dimensions to models of spins [7; 8; 9] and bosons [10; 11].

Retarded interactions (i.e., nonlocal in imaginary time) impose significant limitations regarding system size, temperature, and parameters. They typically arise from a coupling to bosonic modes (e.g., phonons [12; 13] or spin fluctuations [14]), or from dynamical screening [15]. Because such problems are generically nonintegrable, much of their understanding comes from numerical investigations. While impurity problems can be solved very efficiently [16], lattice problems remain a challenge. In the DMRG [17; 18; 19], the large bosonic Hilbert space becomes a limiting factor, especially for finite-temperature or dynamical properties. For SSE-based QMC methods [20; 21; 22; 23; 24; 25; 26], the absence of global updates for the bosons (except for a specific form of spin-phonon coupling [27]) makes simulations significantly less efficient than for fermions or spins. Autocorrelation times increase strongly near phase transitions [21]. Extremely long autocorrelation times also affect determinant QMC methods, where even the sampling of free bosons can be challenging [28]. In fact, for electron-phonon models, the continuous-time interaction-expansion (CT-INT) QMC method [29; 30; 31] with a $(\beta L)^{3}$ scaling and local updates produces better results for correlators at the same numerical cost [23]. On the other hand, directed-loop methods remain efficient for interactions that are long-ranged in space [32; 33].

In this Letter, we overcome these limitations by formulating the problem in imaginary time, integrating out the bosons analytically, and using directed-loop updates to efficiently sample the resulting problem with a retarded interaction. This novel approach combines the advantages of global updates available in the SSE representation and the analytical integration over the bosons possible in the action-based CT-INT method.

Method.—The SSE representation [3] corresponds to a high-temperature expansion of the partition function,

[TABLE]

The Hamiltonian $\hat{H}$ is written as a sum of local operators, $\hat{H}=-\sum_{a,b}\hat{H}_{a,b}$, where $a$ specifies an operator type and $b$ the bond between sites $i(b)$ and $j(b)$. The expansion (1) is sampled stochastically. A configuration with expansion order $n$ corresponds to a string of $n$ operators, specified by the index sequence $S_{n}=\{[a_{1},b_{1}],\dots,[a_{n},b_{n}]\}$, and a state $\ket{\alpha}$ from a complete basis in which $\hat{H}$ is nonbranching, i.e., $\hat{H}_{a,b}\ket{\alpha}\sim\ket{\alpha^{\prime}}$. If $\hat{H}$ contains both off-diagonal ($a=1$) and diagonal ($a=2$) operators, two types of updates are sufficient to achieve ergodicity. For the diagonal updates, it is convenient to fix the length of the operator string to $N$ by inserting $N-n$ unit operators ($a=0$). Then, $S_{N}$ can be traversed sequentially and updates $\hat{H}_{0,b}\leftrightarrow\hat{H}_{2,b}$ be proposed. The corresponding configuration weights are directly obtained from the propagated state $\ket{\alpha(p)}\sim\prod_{l=1}^{p}\hat{H}_{a_{l},b_{l}}\ket{\alpha}$. The global directed-loop updates interchange diagonal and off-diagonal operators on an extensive number of bonds [5].

For concreteness, we explain our method for the highly nontrivial 1D spinless Holstein Hamiltonian [12]

[TABLE]

with the electronic hopping $\hat{B}_{i,i+1}=(\hat{c}^{\dagger}_{i}\hat{c}^{\vphantom{\dagger}}_{i+1}+\mathrm{H.c.})$ and the density $\hat{\rho}^{\vphantom{\dagger}}_{i}=(\hat{c}^{\dagger}_{i}\hat{c}^{\vphantom{\dagger}}_{i}-1/2)$; $\hat{c}^{\dagger}_{i}$ and $\hat{a}^{\dagger}_{i}$ are the usual fermionic and bosonic creation operators acting at lattice site $i$; the chemical potential is zero. In the existing SSE approach to fermion-boson models [20; 21; 22; 23], the fermions are sampled as explained above. Because of the absence of off-diagonal terms, the bosons are updated by local moves with a cutoff for the local occupation number. Even using tempering, autocorrelation times increase significantly with the coupling $\gamma$ [21], and acceptance rates decrease exponentially for large $\omega_{0}$ [23], which severely restricts applications. Our new approach eliminates these problems. It is based on the coherent-state path integral [34], where the Gaussian integrations over the bosonic fields are carried out [35] to obtain a fermionic action with retarded interaction [we define $\lambda=\gamma^{2}/(2\omega_{0}t)$]

[TABLE]

$P(\tau)=\omega_{0}\cosh[\omega_{0}(\beta/2-\tau)]/[2\sinh(\omega_{0}\beta/2)]$ is the free boson propagator with $\tau\in[0,\beta)$ and $P(\tau+\beta)=P(\tau)$. Similar interactions arise for other fermion-boson models.

In an action-based formulation, the SSE representation corresponds to an expansion of $Z=\int\mathcal{D}(\bar{c},c)\,e^{-\mathcal{S}_{0}-\mathcal{S}_{1}}$ around $\mathcal{S}_{0}=\int d\tau\sum_{i}\bar{c}_{i}(\tau)\,\partial_{\tau}\,c_{i}(\tau)$. For a general action, we write $\mathcal{S}_{1}$ as a sum over vertices,

[TABLE]

A vertex is specified by a superindex $\nu$, a weight $w_{\nu}$, and the Grassmann representation $h_{\nu}$ of an operator. The partition function becomes

[TABLE]

where $C_{n}=\{\nu_{1},\dots,\nu_{n}\}$ encodes a configuration of order $n$, $\left\langle O\right\rangle_{0}=Z_{0}^{-1}\int\mathcal{D}(\bar{c},c)\,e^{-\mathcal{S}_{0}}O$ with $Z_{0}=\int\mathcal{D}(\bar{c},c)\,e^{-\mathcal{S}_{0}}$, and time-ordering is implicit. The expansion (5) converges for any finite $\beta$ and $L$ [29].

For problems without retardation, we have $\mathcal{S}_{1}=-\int d\tau\sum_{a,b}H_{a,b}(\tau)$, i.e., $\nu=\{a,b,\tau\}$, $w_{\nu}=d\tau$, and $h_{\nu}=H_{a,b}(\tau)$. The relation between Eq. (5) and the SSE representation is established by mapping the time-ordered expectation value to an operator string:

[TABLE]

On the right-hand side, the time labels are obsolete. Therefore, the $\tau$-integrations contained in $\sum_{C_{n}}$ can be carried out and give $\beta^{n}$, leading to Eq. (1). A mapping between the SSE and a time-ordered interaction expansion in the full Hamiltonian was introduced in Ref. [25].

For retarded interactions such as Eq. (3), $\mathcal{S}_{1}$ contains in particular the interaction $\mathcal{S}_{\text{ret}}$; $\left\langle h_{\nu_{1}}\dots h_{\nu_{n}}\right\rangle_{0}$ can still be mapped to an operator string to calculate the weight of a configuration. However, the fact that the weight $w_{\nu}$ depends on imaginary time demands an explicit sampling of the $\tau$-integrals, as well as time-ordering of the fields. Since $\mathcal{S}_{1}$ consists of bilinears $\overline{c}(\tau)c(\tau)$, this reordering does not change the sign of the configuration.

Algorithm.—We illustrate our algorithm for the spinless Holstein model; other models (see below) only require minimal modifications. As the standard directed-loop method is well documented [5], we focus on the differences. For technical details see also the SM.

(i) Configuration space.—A configuration consists of a state $\ket{\alpha}=\ket{n_{1},\dots,n_{L}}$ in the local occupation number basis, an expansion order $n$, and an ordered vertex list $C_{n}=\{\nu_{1},\dots,\nu_{n}\}$. For a coherent representation of both the bilinear hopping terms with one time argument and the biquadratic interaction terms with two time arguments, it is expedient to write each vertex $\nu_{k}$ as two subvertices $\nu_{k,1}$ and $\nu_{k,2}$ and add a unit “operator” $\mathbb{1}_{b}(\tau)$ with a dummy time variable to the hopping terms. Each subvertex then has local variables $\nu_{k,j}=\{a_{k,j},b_{k,j},\tau_{k,j}\}$ (see Fig. 1). For the Holstein model, $b_{k,1}=b_{k,2}=b_{k}$. To lighten the notation we drop the index $k$ from here on. We write Eq. (4) as

[TABLE]

The off-diagonal hopping vertices are given by

[TABLE]

whereas the diagonal interaction vertices read

[TABLE]

with $j(b)=i(b)+1$. To arrive at the form (7), we multiplied the off-diagonal terms in Eq. (8) with the bosonic propagator and exploited $\int_{0}^{\beta}d\tau_{2}\,P(\tau_{1}-\tau_{2})=1$ for the dummy time variables. This essentially promotes the hopping terms to retarded interactions and yields a vertex weight $\mathcal{W}_{\nu}=w(\tau_{1},\tau_{2})\,W[h_{a_{1}a_{2},b}(\tau_{1},\tau_{2})]$. Here, $w(\tau_{1},\tau_{2})=P(\tau_{1}-\tau_{2})\,d\tau_{1}d\tau_{2}$ irrespective of the operator types $a_{1}$, $a_{2}$. As a result, $P(\tau_{1}-\tau_{2})$ only plays a role for the diagonal updates but drops out of the directed-loop equations, allowing for a simple and efficient implementation. In contrast, $W[h_{a_{1}a_{2},b}(\tau_{1},\tau_{2})]$ depends on time only implicitly via the world line configuration and its values are given in the SM. Finally, the constant $C=1/2+\delta$ ($\delta\geq 0$) in Eq. (9) ensures positive weights.

(ii) Diagonal updates.—For retarded interactions, the operator string cannot be traversed sequentially because each vertex update requires knowledge of the propagated state at two distinct positions in the string. However, the occupation number at $\{i,\tau\}$ is completely determined by the initial state $\ket{\alpha}$ and the number of off-diagonal operators that act between [math] and $\tau$ and involve site $i$. During the diagonal updates, we construct an ordered list containing the time arguments of the operators $B_{b(i)}(\tau)$ for each $i$. Sorting this list takes $\mathcal{O}(L\beta\log\beta)$ operations, after which any propagated state can be quickly calculated. Diagonal updates involve adding or removing a single vertex $h_{22,b}(\tau_{1},\tau_{2})$ using the Metropolis-Hastings algorithm [36; 37] with acceptance rates $A_{C\to C^{\prime}}=\min\left(R_{C\to C^{\prime}},1\right)$. For the addition of a new vertex $R_{C_{n}\to C_{n+1}}={L\beta W[h_{22,b}(\tau_{1},\tau_{2})]}/{(n_{\text{diag}}+1)}$, whereas for the removal $R_{C_{n}\to C_{n-1}}=1/R_{C_{n-1}\to C_{n}}$. Here, $n_{\text{diag}}$ is the number of diagonal vertices in $C_{n}$. Sampling $\tau_{1}$, $\tau_{2}$ according to $P(\tau_{1}-\tau_{2})$ by inverse transform sampling ensures high acceptance rates for any $\omega_{0}$.

(iii) Directed-loop updates.–Directed-loop updates are very similar for retarded and instantaneous interactions [5]. In the latter case, $\ket{\alpha(p)}$ is updated along a closed path connecting a subset of the vertices. Starting at a leg $l_{\text{i}}$ of a randomly chosen vertex, the choice of the exit leg $\l_{\text{e}}$ determines how the vertex changes as a result of the flipping of the occupation numbers $n_{l_{\text{i}}},n_{l_{\text{e}}}$ to $1-n_{l_{\text{i}}},1-n_{l_{\text{e}}}$. Thereby, the operator type of the vertex can change from $a=1$ to $a=2$ or vice versa. From $\l_{\text{e}}$ the loop continues to the next vertex until it closes. The probabilities for choosing $l_{\text{e}}$ are determined by the directed-loop equations for a general vertex [5], which can be derived from the requirement of local detailed balance.

Our generalization to retarded interactions exploits (i) the subvertex structure introduced above, (ii) the fact that the update of a subvertex only changes the world line configuration locally into another allowed configuration, and (iii) our choice of the weight $w(\tau_{1},\tau_{2})$ that removes any time dependence from the directed-loop equations. Because of (i) and (ii) each subvertex becomes an independent entry to the usual linked vertex list [5] that also includes the unit operators. While (ii) allows us to update subvertices individually, the retarded interaction (3) leads to an update probability that also depends on the other subvertex connected via $P(\tau)$. These conditions hold for the Holstein model (see SM) but also for, e.g., Fröhlich, Su-Schrieffer-Heeger, and spin-phonon models [38]. Finally, for the spinless Holstein model, the directed-loop equations can be solved exactly and backtracking is absent for $\lambda\leq 1$ (see SM).

(iv) Observables.—Electronic observables are calculated exactly as in the SSE representation [39]. Bosonic estimators are obtained using generating functionals [40]. Dynamic correlation functions are also accessible.

Application.—To demonstrate the potential of our new method, we first discuss its efficiency. In standard SSE simulations of the Holstein-Hubbard model—the spinful analog of Eq. (2)—the integrated autocorrelation time $\tau_{\text{int}}$ essentially diverges with $\lambda$ [21]. Although reduced by parallel tempering, $\tau_{\text{int}}$ exceeds $100$ at intermediate coupling already for moderately difficult parameters ($\omega_{0}=t$, $L=16$, $\beta t=2L$) [21]. Figure 2 shows $\tau_{\text{int}}$ for our method for $L=18$ and $\beta t=2L$ 111A sweep consisted of two blocks of diagonal and directed-loop updates. For each block of diagonal updates, we attempted approximately $2\left\langle n_{\text{diag}}\right\rangle$ updates. The number of loop updates was fixed by touching approximately $2\left\langle n\right\rangle$ subvertices of type $a=1,2$., covering the entire range of phonon frequencies from adiabatic to antiadiabatic and the entire range of couplings from weak to strong. Periodic boundary conditions were used. Remarkably, $\tau_{\text{int}}$ is of order 1 both for the spinful and the spinless Holstein model 222For the spinful Holstein model, each subvertex obtains an additional spin variable $\sigma_{j}$. Including this variable for the dummy unit “operators” in the off-diagonal vertices leads to an additional prefactor of $1/2$.. Autocorrelations in fact decrease with increasing $\lambda$, with no visible signature of the Peierls phase transition. The data shown are for the total energy, for other observables $\tau_{\text{int}}$ is even smaller. Similar autocorrelation times were observed for larger systems.

Having established its numerical efficiency, we used the method to obtain high-precision results for the half-filled spinless Holstein model (2). The latter provides a generic framework to study the Peierls transition of 1D electrons coupled to quantum phonons. From previous work [45; 46; 47; 48; 44], the model is known to exhibit a Berezinskii–Kosterlitz–Thouless quantum phase transition with dynamical exponent $z=1$ between a Luttinger liquid and a charge-density-wave (CDW) insulator with a $q=2k_{\text{F}}=\pi$ modulation of charge density and lattice deformations. Since $z=1$ we keep $\beta/L=\text{const}$.

Figure 3 shows the real-space density correlator $C_{\rho}(r)=\langle\hat{\rho}_{r}\hat{\rho}_{0}\rangle$ (using the conformal distance $\xi$, see caption) and the fidelity susceptibility $\chi_{\text{F}}$ [50; 40], a finite-temperature extension of the quantum fidelity and an unbiased diagnostic for quantum phase transitions [51; 52]. We simulated systems of up to $L=1282$ sites with $\beta t\geq 2L$. Real-space correlation functions were previously reported for $L\lesssim 50$ [23], DMRG results for other quantities were available up to $L=256$ [44].

Figure 3(a) reveals the theoretically predicted power-law decay of density correlations in a spinless, repulsive Luttinger liquid [53]. The dominant contribution to $C_{\rho}(r)$ is the oscillating term $\cos(2k_{\text{F}}r)\,r^{-2K}$ (we only plot even $r$). The nonuniversal exponent is determined by the Luttinger parameter $K$. As expected for the Mott transition of a spinless Luttinger liquid, $K$ decreases with increasing $\lambda$ until it reaches the critical value $K=1/2$ for $\lambda_{c}=0.68(1)$ [44]. This can be seen by comparing to the dashed line in Fig. 3(a) that shows a $1/r$ power law. The inset shows estimates for $K$ from power-law fits (see caption for details). For $\lambda>\lambda_{c}$, $K$ scales to zero and the system exhibits long-range CDW order.

In Fig. 3(b), we plot the density correlator at the largest distance $r=L/2$, whose thermodynamic limit serves as an order parameter for the quantum phase transition. We find a nonzero extrapolated order parameter for $\lambda\gtrsim 0.68$, in accordance with Fig. 3(a) and previous estimates [44]. The transition can also be detected from the fidelity susceptibility shown in Fig. 3(c). Because statistical errors are generally larger for $\chi_{\text{F}}$, the maximum system size was $L=162$. In contrast to previous work, the directed-loop algorithm permits us to reach sufficiently large values of $L$ and $\beta$ to observe the cusp at $\lambda_{c}$ predicted theoretically [54]. The latter sharpens and converges (slowly, similar to other 1D models [54]) to $\lambda_{c}$ with increasing $L$. More generally, Figs. 3(b) and 3(c) are important because they establish the usefulness of the order parameter and $\chi_{\text{F}}$ to detect the CDW transition without reference to bosonization results. They can therefore be used for spinful electron-phonon models, the analysis of which is complicated by the existence of a spin gap in the metallic phase [55; 56; 23] and the absence of a reliable theory for the Mott transition of a Luther-Emery liquid [55]. Moreover, our method can access the system sizes necessary to resolve the spin gap [23].

Conclusions and Outlook.—We have introduced a highly efficient directed-loop QMC method for systems with retarded interactions. For the electron-phonon models considered, our algorithm outperforms any other existing method, including the DMRG. Because of the global updates, autocorrelations are negligible and there is no need for tempering or machine learning. The method permits us to study fermion-boson models with the same accuracy as purely fermionic models. It can be used to solve a number of open problems in the field of electron-phonon physics. These include the phase diagrams of models with competing electron-phonon and electron-electron interactions, the specific heat of quantum Peierls chains, and the dimensional crossover as a function of temperature in Peierls materials. The method can also be extended to spin-boson models via a fermionic path-integral representation with a suitable constraint. Sign-free simulations of spins or hardcore bosons with retarded interactions can be carried out in any dimension and are important for the understanding of correlated quantum systems with dissipation (for recent numerical work see, e.g., Refs. [57; 58]). Finally, it will be interesting to explore if the method permits us to study quasi-1D materials beyond the low-energy regime by including higher bands via a frequency-dependent interaction.

I Supplemental Material

In the following, we provide details on our implementation of the directed-loop algorithm for the half-filled spinless Holstein model that may be helpful for readers wishing to implement the method themselves or modify an existing directed-loop code.

I.1 Vertex weights

In the SSE representation, the weight of a Monte Carlo configuration, $W(C_{n})=\frac{1}{n!}\prod_{p=1}^{n}\mathcal{W}_{\nu_{p}}$, factorizes into a product of individual vertex weights $\mathcal{W}_{\nu}=w(\tau_{1},\tau_{2})\,W[h_{a_{1}a_{2},b}(\tau_{1},\tau_{2})]$. The explicit time dependence of the vertex is in $w(\tau_{1},\tau_{2})=P(\tau_{1}-\tau_{2})\,d\tau_{1}d\tau_{2}$, which is independent of the operator type and therefore has to be considered only during the diagonal updates. The remainder $W[h_{a_{1}a_{2},b}(\tau_{1},\tau_{2})]=W_{v_{1},v_{2}}$ is fully determined by the vertex types $v_{1},v_{2}\in\{1,\dots,6\}$ that in turn specify the change of the world-line configuration at each subvertex. Figure 4 shows the possible subvertex types for the spinless Holstein model, where $v\in\{1,\dots,4\}$ corresponds to unit and diagonal operators ($a=0,2$) and $v\in\{5,6\}$ to off-diagonal ones ($a=1$). The corresponding weights $W_{v_{1},v_{2}}$ are given in Table 1.

I.2 Solution of the directed-loop equations

For the directed-loop updates, the configuration space of vertex types $v\in\{1,\dots,6\}$ shown in Fig. 4 is enlarged by assigning to a given vertex directed paths that connect an entrance leg $l_{\text{i}}\in\{1,\dots,4\}$ with an exit leg $l_{\text{e}}\in\{1,\dots,4\}$. The Monte Carlo weights of these assignments are determined by the directed-loop equations, which can be derived from the requirement of local detailed balance [5]. For retarded interactions, each vertex consists of two subvertices. While the directed-loop equations determine the weight of the total vertex, the loop is constructed locally by only assigning a directed path to one subvertex and leaving the other unchanged. As shown in Fig. 5(a), the occupation numbers on the sites included in the loop are then switched and the vertex type changes. For the vertices defined by Eqs. (8) and (9) of the main text, we have to distinguish two cases: the directed loop either hits a unit operator, or any other operator. For the former case, the path goes straight through the subvertex with probability $1$ and changes its vertex type. For the latter case, the directed-loop equations have to be solved explicitly.

As discussed in the main text, the world-line configuration is updated at each subvertex independently. However, while the vertex type of the other subvertex does not change, the same is not in general true for its dummy operator type $a_{2}$. For example, $a_{2}$ changes from [math] to $2$ when interchanging $h_{10,b}(\tau_{1},\tau_{2})\leftrightarrow h_{22,b}(\tau_{1},\tau_{2})$ in Fig. 5(a). This corresponds to the update from a hopping operator at $\tau_{1}$ and a unit operator at $\tau_{2}$ to a (diagonal) density-density interaction term at times $\tau_{1}$ and $\tau_{2}$. Because unit operators that change into diagonal operators are relevant for the weights $W_{v_{1},v_{2}}$ in later updates, it is important to keep track of such changes.

We illustrate the solution of the directed-loop equations for the assignment table given in Fig. 5(b) (for a detailed discussion of assignment tables see Ref. [5]). We only show the possible assignments for vertex type $v_{1}$, corresponding to the lower subvertex in Fig. 5(a). The second subvertex of type $v_{2}$ remains unaffected by this segment of the loop update but still enters the configuration weight. Each row in Fig. 5(b) shows the possible assignments for a fixed $l_{\text{i}}$ and the three possible exit legs $l_{\text{e}}$. The associated weights are symmetric around the diagonal because the corresponding assignments are related by inverting the direction of the path and flipping the occupation numbers on the sites touched by the loop. For the specific example of Fig. 5(b), we obtain for the corresponding weights

[TABLE]

The bounce weights $b_{i}$, $i\in\{1,2,3\}$, are related to the assignments on the diagonal, whereas $a$, $b$, and $c$ are the remaining weights. Our goal is to reduce the bounce weights and solve for $a$, $b$, and $c$. To this end, we write

[TABLE]

For concreteness, we choose $v_{2}=3$ and insert the weights given in Table 1. This leads to

[TABLE]

The bounce weights $b_{i}$ and the constant $C=1/2+\delta$ must be chosen such that $a$, $b$, and $c$ are positive. For $\lambda<1$, this is already fulfilled by $b_{1}=b_{2}=b_{3}=0$ and $\delta\geq(1-\lambda)/(4\lambda)$. In our simulations, we have chosen the lower bound. For $\lambda\geq 1$, the positivity of $c$ requires $b_{1}\geq(\lambda-1)t/2$, whereas the positivity of $b$ demands $b_{1}\leq(\lambda+1)t/2$. We have chosen the lower bound and $\delta=0$. This procedure has to be repeated for each type of background vertex $v_{2}$ and each possible assignment table for $v_{1}$. In the end, we find that the global constant $C$ has to be chosen as for the example given here.

Let us point out that exact solutions of the directed loop equations are neither common nor necessary for efficient simulations. Instead, the equations can be solved using linear programming techniques [6].