Inchworm Quantum Monte Carlo

• Inchworm Quantum Monte Carlo is a diagrammatic Monte Carlo method that incrementally extends propagators to reorganize series expansions and mitigate the dynamical sign problem.
• It leverages causal structure and recursive formulation to transform the computational scaling from exponential to quadratic in real-time as well as linear or weakly superlinear in inverse temperature for imaginary-time simulations.
• Variants like SBCI, DCCI, and steady-state implementations extend iQMC to diverse systems including multi-orbital impurity models, open quantum systems, and DMFT, ensuring robust and scalable quantum simulations.

Inchworm Quantum Monte Carlo (Inchworm QMC, often iQMC) denotes a family of diagrammatic Monte Carlo methods in which a perturbative expansion is reorganized so that propagators are extended incrementally in short intervals while previously computed restricted propagators are reused. The defining step is an “inchworm” recursion: diagrams wholly contained in earlier intervals are absorbed into dressed propagators, and only the remaining “inchworm proper” contributions that connect to the newly added interval are sampled explicitly. In the original real-time formulation this largely overcomes the dynamical sign problem, changing the scaling from exponential to quadratic, and subsequent developments extended the same principle to imaginary time, steady-state formulations, multi-orbital impurity models, open quantum systems, DMFT, and quasi-Monte Carlo integration (Cohen et al., 2015, Chen et al., 2016, Erpenbeck et al., 2022, Erpenbeck et al., 2024).

1. Origin and conceptual basis

The method first appeared as an answer to the central obstruction in real-time diagrammatic Monte Carlo: the dynamical sign problem. In conventional real-time expansions, the perturbation order needed to reach a target time grows with that time, while oscillatory factors on the Keldysh contour cause the average sign to decay exponentially. The original Inchworm algorithm exploits the causal structure of contour propagators and reuses information obtained at earlier times to extend the propagation to longer times, thereby changing the computational scaling from exponential to quadratic in the propagation time (Cohen et al., 2015).

The early methodological development quickly bifurcated into several related lines. For the spin-boson model, two expansions were formulated: a system-bath coupling inchworm (SBCI) and a diabatic-coupling cumulant inchworm (DCCI), the latter motivating a cumulant formulation with improved scaling (Chen et al., 2016). A companion benchmark study established that the two formulations are complementary in parameter space and that at least one of them converges rapidly across a broad range of couplings, bath cutoffs, and temperatures (Chen et al., 2016). In parallel, the method was generalized to the full Keldysh contour with forward, backward, and equilibrium branches, allowing the computation of Green’s functions, spectral functions, and currents for quantum impurities out of equilibrium (Antipov et al., 2016).

A second line of development clarified the mathematical structure of the method. For open quantum systems coupled to Gaussian baths, inchworm resummation was established rigorously on the basis of Dyson-series rearrangement, and an integro-differential formulation was introduced that makes the method resemble a memory-kernel evolution equation with a linked-diagram kernel (Cai et al., 2018). That perspective also enabled deterministic time-integration schemes, such as Heun’s method, to be incorporated into inchworm propagation.

A third line broadened the domain of applicability. Real-time inchworm became a DMFT impurity solver on the L-shaped Keldysh contour for the infinite-coordination Bethe lattice (Dong et al., 2017). In imaginary time, the same incremental strategy was adapted to multiorbital impurity models with general interactions and hybridizations, where it suppresses the severe temperature-dependent sign problem of conventional CT-HYB for off-diagonal hybridizations (Eidelstein et al., 2019). Later work introduced an interaction-expansion inchworm solver for impurity and lattice models (Li et al., 2022), a direct steady-state formulation that searches for a time-translation-invariant fixed point instead of propagating through the transient (Erpenbeck et al., 2022), and a steady-state multi-orbital construction that combines the steady-state and multiorbital inchworm ideas in a numerically exact framework (Erpenbeck et al., 2024).

2. Diagrammatic and recursive structure

The formal starting point depends on the chosen expansion, but the basic structure is common. In hybridization-expansion impurity problems one writes

$H = H_S + H_B + H_{SB},$

with $H_S$ acting on a finite impurity Fock space and $H_B$ a noninteracting bath. In the interaction picture with $H_0=H_S+H_B$, the restricted propagator can be expanded as

$$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$

and expressed as a time-ordered series in insertions of $H_{SB}(\tau)$; Wick contraction of bath operators produces determinants or products of hybridization functions such as $\Delta^<$ and $\Delta^>$ (Erpenbeck et al., 2024).

The inchworm step begins by assuming that the propagator is known up to a cutoff $T$. One then writes the propagator at a later time $t>T$ as a sum of classes: bare propagation from $H_S$0 to $H_S$1, terms with one hybridization bridge between $H_S$2 and $H_S$3, and higher-order terms containing all remaining inchworm-proper diagrams. Diagrammatically, only contractions with at least one vertex before and after $H_S$4 are included; contractions entirely before or after the cutoff are omitted because they are already contained in the previously computed dressed propagators. In the formulation of the multi-orbital steady-state solver this yields a recursion of the form

$H_S$5

while the two-branch propagator satisfies an analogous map involving both $H_S$6 and $H_S$7 at earlier times (Erpenbeck et al., 2024).

The original real-time formulation can be summarized in the same language. If a contour propagator $H_S$8 is known up to an intermediate time $H_S$9, then $H_B$0 for $H_B$1 is decomposed into a piece with no new vertices beyond $H_B$2, a no-crossing contribution in the new interval, and a crossing contribution with at least one hybridization line spanning $H_B$3. Only the latter two need to be sampled stochastically, because the earlier-time subsegments are absorbed into the known propagators (Cohen et al., 2015).

In the open-system formulation, taking the continuous-step limit leads to an integro-differential equation. The reduced propagator $H_B$4 obeys an equation of the form

$H_B$5

where the kernel $H_B$6 collects irreducible linked diagrams that straddle the interval $H_B$7 (Cai et al., 2018). This form makes explicit that inchworm is not only a Monte Carlo sampling rule but also a reorganization of the dynamics into a memory equation with dressed short-time propagators.

3. Main algorithmic families

The best-known family is the real-time hybridization-expansion impurity solver on the Keldysh contour. Here one computes restricted propagators for impurity many-body states, marches forward on a discretized contour, and later measures occupations, contour-ordered Green’s functions, currents, and spectral functions. At finite truncation in the number of crossing hybridization lines, the method reproduces the $H_B$8-crossing approximations, with $H_B$9 corresponding to NCA and $H_0=H_S+H_B$0 to OCA; the full inchworm CT-HYB is recovered as $H_0=H_S+H_B$1 with sampling up to a maximum total order $H_0=H_S+H_B$2 (Antipov et al., 2016).

A distinct but related family is built for non-adiabatic dynamics in bosonic environments. In the SBCI formulation the perturbation is the system-bath coupling, whereas in DCCI the perturbation is the diabatic coupling $H_0=H_S+H_B$3 and the bath is treated through cumulant resummation. DCCI yields a memory equation in terms of COP-cumulants $H_0=H_S+H_B$4, and the inchworm version samples only those cumulant partitions that straddle the new interval. The benchmark study shows that SBCI converges best at small $H_0=H_S+H_B$5 or large $H_0=H_S+H_B$6, while DCCI converges more easily for large $H_0=H_S+H_B$7 (Chen et al., 2016, Chen et al., 2016).

Imaginary-time inchworm for multiorbital impurity models adapts the same incremental logic to the interval $H_0=H_S+H_B$8. Instead of sampling the full CT-HYB expansion over the entire inverse-temperature interval, the algorithm computes partial propagators on a fine imaginary-time grid and reuses them at later steps. The aim is not a real-time dynamical sign problem but the severe sign deterioration with $H_0=H_S+H_B$9 caused by off-diagonal hybridizations and general interactions. The formal local Hilbert-space size remains exponential in the number of orbitals, but the temperature scaling is transformed from exponential sign degradation to behavior that is linear or weakly superlinear in $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$0 (Eidelstein et al., 2019).

The interaction-expansion inchworm formulation uses the expansion in the two-body interaction $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$1 rather than in the hybridization. It introduces auxiliary objects $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$2 and $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$3, where all interaction vertices are restricted to $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$4, and constructs $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$5 from $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$6 by sampling only those connected diagrams containing at least one vertex in $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$7 while excluding two-particle reducible Type-1 subgraphs. One important distinction from bold-line skeleton methods is explicit in this literature: the problem of convergence to unphysical fixed points, which hampers so-called bold-line methods, is not encountered in inchworm Monte Carlo (Li et al., 2022).

Steady-state inchworm departs more radically from transient propagation. In the nonequilibrium steady state one assumes time-translation invariance and, in the steady-state impurity formulations, $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$8-independence of the two-branch restricted propagator. The algorithm then iterates a self-consistent fixed-point map for $$\phi^\beta_\alpha(t)=\mathrm{Tr}_B\{\rho_B \langle \beta|e^{-iHt}|\alpha\rangle\}$$9 directly on a grid of time differences up to a cutoff, sometimes called the coherence time, instead of stepping through long transients (Erpenbeck et al., 2022, Erpenbeck et al., 2024). The 2024 multi-orbital implementation applies this construction in an atomic-state basis of size $H_{SB}(\tau)$0 and is designed specifically to mitigate both the dynamical and multi-orbital sign problems (Erpenbeck et al., 2024).

4. Observables and application areas

The observable content of inchworm QMC is broader than propagator evaluation alone. In the full-Keldysh impurity formulation, single-time observables such as occupations are obtained directly from restricted propagators, while contour-ordered two-time Green’s functions are measured by sampling diagrams with additional operator insertions. The current from lead $H_{SB}(\tau)$1 can then be written as

$H_{SB}(\tau)$2

and the spectral function can be extracted either from real-time Green’s functions or, in one formulation, from currents through auxiliary leads attached at frequency $H_{SB}(\tau)$3 (Antipov et al., 2016).

For transport and full counting statistics, arbitrary lead geometries enter only through the hybridization functions. The counting-field extension modifies the lead coupling on the Keldysh contour and yields a moment-generating function $H_{SB}(\tau)$4, from which the current and noise are obtained as long-time slopes of the first and second cumulants. This framework was used to study charge transport and its fluctuations through molecular junctions with finite-bandwidth leads (Ridley et al., 2019).

In DMFT, inchworm serves as an impurity solver inside the self-consistency loop. On the infinite-coordination Bethe lattice the equilibrium real-time self-consistency relation takes the particularly simple form

$H_{SB}(\tau)$5

Once the impurity Green’s function is obtained on the real-time contour, the retarded self-energy can be extracted from a discretized contour-Dyson equation, and direct Fourier transformation of $H_{SB}(\tau)$6 bypasses analytic continuation (Dong et al., 2017). The steady-state formulation was later used for correlated materials driven by a bias voltage, where the response of a correlated material was shown to differ qualitatively from the splitting of the Kondo resonance in bias-driven quantum dots (Erpenbeck et al., 2022).

Multi-orbital impurity physics is a major application domain because standard CT-HYB is often restricted by sign problems for off-diagonal hybridizations. The equilibrium multiorbital solver addresses general interactions and hybridizations without truncating the impurity bath to a discrete set of levels (Eidelstein et al., 2019). The steady-state multiorbital method extends this capability to equilibrium and nonequilibrium steady states with finite bias or thermal gradients, and computes Green’s functions by an inchworm-proper expansion analogous to that for the steady-state propagator (Erpenbeck et al., 2024).

The method has also been generalized beyond fermionic impurity problems. For bosonic baths, inchworm was formulated for open quantum systems and spin chains, where the Dyson expansion is reorganized into an integro-differential evolution for reduced propagators or single-spin half-propagators (Cai et al., 2018, Wang et al., 2023). In a later extension to spin chains with off-diagonal coupling, tensor-train representations and the transfer tensor method were used to control operator growth and long-time memory cost (Sun et al., 2024).

5. Sign problems, scaling, and accelerations

The central computational claim of inchworm QMC is not that sign problems disappear, but that the dominant cancellations are pushed into short intervals and absorbed into dressed propagators. In the original formulation this transforms the cost of reaching time $H_{SB}(\tau)$7 from exponential in $H_{SB}(\tau)$8 to approximately quadratic in $H_{SB}(\tau)$9 (Cohen et al., 2015). In the full-Keldysh impurity implementation, error bars typically grow sub-exponentially, nearly linearly in time, rather than exhibiting the catastrophic exponential blow-up of bare real-time expansions (Antipov et al., 2016).

In multiorbital settings the same principle attacks a different obstruction. Off-diagonal hybridizations in CT-HYB generate a severe sign problem even in equilibrium. The equilibrium multiorbital inchworm solver shows that the error per unit $\Delta^<$0, $\Delta^<$1, remains essentially constant as $\Delta^<$2 increases, implying linear or low-order polynomial scaling in $\Delta^<$3 rather than exponential behavior (Eidelstein et al., 2019). In the 2024 steady-state multiorbital construction, the inchworm expansion resums an infinite set of bare diagrams at each order, so the cancellations causing both the dynamical and multi-orbital sign problems are strongly suppressed; off-diagonal $\Delta^<$4 no longer induce an exponentially hard sign problem because sign-canceling clusters are absorbed into the dressed propagators rather than sampled explicitly (Erpenbeck et al., 2024).

The residual computational bottlenecks depend on the formulation. In the steady-state multiorbital solver, single-branch propagation scales linearly in the number of time steps $\Delta^<$5, and the two-branch steady-state iteration also scales linearly in $\Delta^<$6 per iteration, with roughly $\Delta^<$7–$\Delta^<$8 iterations for convergence. The required inchworm order grows slowly with coherence length, while inclusion-exclusion enumeration of proper diagrams costs $\Delta^<$9 with $\Delta^>$0. Propagator operations scale as $\Delta^>$1 with $\Delta^>$2 (Erpenbeck et al., 2024).

A substantial part of later work focuses on accelerating the stochastic integration itself. In imaginary time, quasi-Monte Carlo was introduced by mapping Sobol low-discrepancy sequences from the hypercube to the ordered simplex with the Root transform. For sufficiently smooth integrands this changes the asymptotic error scaling from the Monte Carlo $\Delta^>$3 law to observed $\Delta^>$4 behavior, and in the single-orbital benchmark produced two orders of magnitude higher accuracy at $\Delta^>$5 (Strand et al., 2023). For bosonic baths, inclusion-exclusion methods reduce the computational complexity of summing pairings from double factorial to exponential, and for rounded-box linked diagrams the optimized cost is $\Delta^>$6 with $\Delta^>$7 (Yang et al., 2021). More recently, the bath influence functional has been approximated in tensor-train form so that deterministic quadrature can replace Monte Carlo sampling of high-dimensional ordered integrals; in that construction the resulting complexity scales linearly with the number of dimensions and couples naturally to the transfer tensor method for long-time propagation (Wang et al., 14 Jun 2025).

6. Benchmarks, misconceptions, and limitations

Benchmark studies consistently show that inchworm reorganizations extend accessible times, temperatures, or orbital complexity far beyond the corresponding bare expansions. For the Anderson impurity model in the Kondo regime, bare real-time expansion remained accurate only to about $\Delta^>$8, whereas inchworm QMC was stable to $\Delta^>$9 with roughly constant error bars (Cohen et al., 2015). In the voltage-quench study on the full Keldysh contour, inchworm with $T$0–$T$1 yielded converged populations and currents with $T$2–$T$3 error bars and captured the build-up and splitting of the Kondo resonance at late times (Antipov et al., 2016). In the imaginary-time quasi-Monte Carlo study of a two-orbital Kanamori model with off-diagonal hybridization, the average sign was found to be approximately $T$4 within Monte Carlo noise, whereas CTHYB can have $T$5; at $T$6, inchworm QMC at order $T$7 already lay on top of the exact diagonalization reference for $T$8 (Strand et al., 2023).

The multi-orbital steady-state solver provides a particularly direct validation of the modern framework. For two decoupled spin-full orbitals, one noninteracting and one with $T$9, the method reproduces both the analytic noninteracting spectral function and the single-orbital steady-state inchworm result, with convergence around $t>T$0. For a fully interacting two-orbital model at a high-symmetry point that diagonalizes into two single-orbital problems, diagonal and off-diagonal spectral functions converge to the exact solutions by $t>T$1. Benchmarks on two orbitals report accurate spectral functions in roughly $t>T$2–$t>T$3 Monte Carlo samples (Erpenbeck et al., 2024). In DMFT for the bilayer Hubbard model, relative errors of CT-HYB and inchworm were compared at fixed wall time: at $t>T$4, CT-HYB was reported as completely overwhelmed, with errors larger than $t>T$5, whereas inchworm still delivered $t>T$6 (Goldberger et al., 2023).

Several common misconceptions are corrected by the literature itself. First, inchworm does not remove every sign problem. The original impurity work states explicitly that it tames the dynamical sign problem, not the intrinsic fermionic sign problem of generic lattice fermions away from half filling (Cohen et al., 2015). Second, inchworm is not merely a skeleton or bold-line fixed-point iteration under another name. In the interaction-expansion formulation, the partially dressed $t>T$7 is advanced by a forward expansion in $t>T$8, and the convergence to unphysical fixed points that can afflict bold-line methods is not encountered (Li et al., 2022). Third, numerically exact does not mean parameter-free. The full-Keldysh implementation requires systematic convergence checks in $t>T$9, $H_S$00, and run-to-run fluctuations, with multiple independent full-contour calculations recommended in order to capture error propagation correctly; in practice $H_S$01–$H_S$02 independent runs are used (Antipov et al., 2016).

The main limitations are therefore structural rather than conceptual. The local Hilbert-space cost remains exponential in the number of orbitals, since propagators are matrices in an atomic-state basis of size $H_S$03 (Erpenbeck et al., 2024). In real-time impurity formulations, the cost still grows exponentially in inverse temperature and in local-Hilbert-space size (Antipov et al., 2016). In imaginary-time quasi-Monte Carlo, factorial growth of diagram topologies restricts the practical expansion order to approximately $H_S$04–$H_S$05 unless additional inclusion-exclusion resummations are used (Strand et al., 2023). For spin-boson problems, the choice between SBCI and DCCI remains parameter dependent: SBCI is efficient at weak system-bath coupling, whereas DCCI excels at strong coupling and large $H_S$06 (Chen et al., 2016, Chen et al., 2016). These constraints explain why inchworm QMC is best understood not as a single algorithm, but as a general resummation strategy whose effectiveness depends on the chosen expansion, observable, and physical regime.