CT-HYB: Hybridization Expansion for DMFT

• Hybridization Expansion (CT-HYB) is a continuous-time quantum Monte Carlo method that stochastically expands the impurity partition function in hybridization, yielding numerically exact solutions for strongly correlated electron systems.
• The method employs pair and multi-operator (e.g., four-operator) updates to overcome ergodicity issues in broken symmetry phases, crucial for accurately capturing phenomena like d-wave superconductivity.
• Practical implementations incorporate advanced techniques such as sliding-window sampling, worm algorithms, and efficient determinant updates to enhance computational efficiency and measurement precision.

Hybridization Expansion (CT-HYB) is a continuous-time quantum Monte Carlo (CT-QMC) impurity solver pivotal in dynamical mean-field theory (DMFT) and its cluster extensions for strongly correlated electrons. CT-HYB stochastically samples the power series expansion of the impurity partition function in the impurity–bath hybridization, giving numerically exact access to single- and multi-orbital impurity problems with general local interactions. It is the method of choice for models with strong interactions and arbitrary interactions, and adapts to complex, multi-orbital, spinful, and cluster problems. Its performance and accuracy have led to its status as a workhorse solver in modern DMFT studies.

1. Diagrammatic Formulation and Partition Function Expansion

At the core of CT-HYB is the expansion of the partition function $Z$ of a general impurity Hamiltonian,

$$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$

where $d_i^\dagger$ and $a_\mu^\dagger$ denote creation operators on impurity and bath, respectively, $V_{\mu i}$ are hybridization amplitudes, and $\epsilon_\mu$ bath energies (Sémon et al., 2014).

Integrating out the noninteracting bath yields the hybridization function,

$$\Delta_{ij}(\tau-\tau') = \sum_\mu V^*_{\mu i} V_{\mu j} \frac{e^{-\epsilon_\mu(\tau-\tau')}}{1 + e^{-\beta \epsilon_\mu}},$$

and the action becomes a nonlocal-in-time impurity problem. Expanding in powers of the hybridization term generates a sum over all configurations of operator insertions and removals (bath excitations), each configuration $c$ comprising $k$ hybridization events specified by times and orbital indices:

$$Z = \sum_{k=0}^\infty \sum_{\{i_j, i'_j\}} \frac{1}{k!^2} \int d\tau_1 \ldots d\tau_k d\tau_1' \ldots d\tau_k' \text{Tr}_\text{imp}[ T_\tau e^{-\beta H_\mathrm{loc}} \mathcal{O} ] \cdot \det[ \Delta_{i_j' i_i}(\tau_j' - \tau_i) ].$$

The local trace and bath determinant together define the Markov chain weights for Monte Carlo sampling (Sémon et al., 2014).

2. Monte Carlo Sampling and Update Schemes

Monte Carlo sampling of the expansion is realized via Metropolis–Hastings moves in the configuration space. The standard moves are insertion/removal of a pair of impurity creation/annihilation operators corresponding to hybridization events:

• Pair insertion: Add $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$0 at times ($$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$1, $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$2). Acceptance ratio involves determinant and local trace ratios, weighted by a phase-space prefactor.
• Pair removal: Remove an existing hybridization line.

The acceptance probability for insertion, for instance, is

$$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$3

using efficient rank-1 Sherman–Morrison updates for determinants (Sémon et al., 2014). For systems with density–density interactions (segment representation), the updates are restricted to nonoverlapping interval insertions/removals, further simplifying trace computation (Nagai et al., 2016).

3. Ergodicity and Broken-Symmetry States

The standard pair-update scheme is nonergodic in quantum impurity problems with broken spatial or superconducting symmetries, such as anomalous baths. For instance, in cluster momentum space, anomalous hybridizations between distinct $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$4-points cannot be sampled via consecutive two-operator updates due to quantum number conservation constraints in $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$5. This deficiency manifests as the inability to reach certain configurations possessing anomalous pairings, e.g., configurations crucial for $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$6-wave superconductivity in cluster DMFT (Sémon et al., 2014).

Restoring ergodicity in such phases requires multi-operator updates. Specifically, the essential move is the four-operator insertion/removal: simultaneously inserting two creation and two annihilation operators in a configuration that collectively conserves all quantum numbers but cannot be constructed from pairwise updates. This extension enables direct sampling of anomalous configurations and is necessary, for instance, for access to d-wave order parameters in plaquette CDMFT (Sémon et al., 2014).

4. Practical Applications and Algorithmic Impact

In practical cluster DMFT studies, e.g. of the one-band Hubbard model, neglecting four-operator updates leads to a significant underestimate of the amplitude and spatial extent of superconducting order. As reported in (Sémon et al., 2014), inclusion of four-operator moves increases the $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$7-wave order parameter by approximately 20% and extends the superconducting dome to higher dopings, yielding results that closely match those obtained using exact diagonalization.

In terms of computational cost, multi-operator (rank-4) determinant updates for four-operator moves scale as $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$8 per attempt, introducing only a modest overhead. The Markov chain mixing time is unaffected, and the overall scaling per Monte Carlo step remains $$H_\mathrm{imp} = H_\mathrm{loc}(d_i^\dagger, d_i) + \sum_{i,\mu} (V_{\mu i} a_\mu^\dagger d_i + V_{\mu i}^* d_i^\dagger a_\mu) + \sum_\mu \epsilon_\mu a_\mu^\dagger a_\mu,$$9, where $d_i^\dagger$0 is the local Hilbert space dimension (Sémon et al., 2014).

5. Extensions and Algorithmic Implementations

Modern CT-HYB implementations include further algorithmic optimizations:

• Sliding-window sampling: Restricts local trace updates to a small interval in imaginary time, reducing effective per-update cost to $d_i^\dagger$1 except at window boundaries (Shinaoka et al., 2014).
• Segment and Krylov representations: Apply depending on the structure of $d_i^\dagger$2; segment for density–density interactions, Krylov (or Newton–Leja, MPS) for full rotationally invariant or off-diagonal interactions (Huang et al., 2012, Shinaoka et al., 2014).
• Worm sampling: Enables unbiased measurement of Green’s functions and higher-order correlators, circumvents missing diagram issues in non-ergodic samplings, and is essential for systems with discrete baths (Gunacker et al., 2015, Hausoel et al., 2022).
• Superstate/state sampling and skip-lists: Drastically improve computational efficiency for multi-orbital and large Hilbert space cases, achieving up to $d_i^\dagger$3–$d_i^\dagger$4 speedups (Kowalski et al., 2018, Sémon et al., 2014).

6. Measurement of Observables and Advanced Estimators

Green’s function measurement in CT-HYB exploits the line-removal estimator, but worm sampling is mandatory to obtain unbiased two-particle correlators and self-energies in the presence of discrete or sparsely coupled baths (Gunacker et al., 2015, Hausoel et al., 2022). Improved estimators for the self-energy and two-particle vertex use equations of motion, requiring measurement of higher-order correlators at negligible additional cost; this reduces noise and renders analytic continuations more stable (Hafermann et al., 2011). Legendre polynomial filtering is commonly employed for noise suppression in the high-frequency regime.

7. Summary Table: CT-HYB Update Types and Ergodicity

Update Type	Ergodic for standard problems	Ergodic for broken symmetries	Typical MC cost per move
Pair (2-operator)	Yes	No	$d_i^\dagger$5
Quadruple (4-operator)	Not necessary	Yes	$d_i^\dagger$6
Higher (6,8,...) operator	Rarely needed	Exotic cases	$d_i^\dagger$7

Incorporating multi-operator moves is essential for ergodicity and quantitative accuracy in phases breaking particle-number, spin, or spatial symmetries, enabling correct sampling of physically relevant configuration space at minimal added computational expense (Sémon et al., 2014).