Diagrammatic Monte Carlo Method

• Diagrammatic Monte Carlo is a stochastic computational approach that evaluates perturbative Feynman diagram expansions to compute exact quantum many-body observables.
• It employs Markov chain sampling over complex diagram configurations to overcome factorial growth and avoid the truncation errors of deterministic methods.
• Variants like bare and bold expansions enable targeted simulations across diverse models, from lattice systems to open quantum and impurity models.

The diagrammatic Monte Carlo (DiagMC) method is a stochastic computational approach for evaluating perturbative series expansions expressed in terms of Feynman diagrams, targeting ground-state and dynamical properties of strongly correlated quantum many-body systems. By recasting sums over diagrams—potentially of high order and factorial complexity—into a Markov chain over configuration space, DiagMC enables numerically exact (i.e., unbiased and systematically improvable) calculations free from the truncation errors that plague deterministic approaches. Distinct DiagMC variants have been developed for quantum lattice models, quantum field theories, molecular electronic structure, polaron physics, quantum impurity and open systems, each tailored to the analytic structure, convergence behavior, and physical observable of interest.

1. Formal Structure of Diagrammatic Monte Carlo

DiagMC is grounded in the expansion of physical observables (e.g., single-particle Green’s functions, self-energies, free energies) in terms of Feynman diagrams built from bare or “bold” (dressed) propagators, interaction vertices, and symmetry factors. In general, for a Hamiltonian $H = H_0 + V$, the expectation value of an operator $O$ may be written as: $$\langle O \rangle = \sum_{n=0}^\infty \frac{(-1)^n}{n!} \int d\tau_1 \cdots d\tau_n \langle T\left[ V(\tau_1) \cdots V(\tau_n) O \right] \rangle_{0,c}$$ Each expansion term corresponds to a set of diagrams determined by the model’s Feynman rules: the types of propagator lines, interaction vertices, and symmetry constraints dictated by the underlying physics (e.g., fermionic or bosonic statistics, vertex connectivity, energy-momentum conservation).

The series may be truncated at a maximum order $n_\text{max}$, but the diagrammatic space typically grows factorially with $n_\text{max}$—rendering deterministic enumeration infeasible. DiagMC circumvents this by random-walking in the space of diagram topologies and internal parameters, collecting statistical estimates of observables. At each Monte Carlo step, a diagram $D$ is proposed and assigned a weight $W(D)$, commonly factored into products of propagators, vertex amplitudes, and overall signs.

2. Variants: Bare vs. Bold (Skeleton) Expansions

A central axis for classification is whether internal lines in sampled diagrams carry bare propagators ($G_0$, $D_0$) or fully self-consistent (“bold”) propagators, as well as whether the expansion targets specific irreducible subclasses.

Bare Expansion: Every internal line carries the non-interacting propagator. The self-energy $\Sigma$ or polarization $\Pi$ are computed by summing over all irreducible diagrams built from $G_0$ and vertex factors. The full Green's function is reconstructed via Dyson’s equation.

Skeleton (Bold) Expansion: Internal lines carry “bold” (dressed) propagators, updated self-consistently from the outcomes of the Monte Carlo walk. For example,

$$G^{-1}(k, i\omega_n) = G_0^{-1}(k, i\omega_n) - \Sigma(k, i\omega_n)$$

with $\Sigma$ defined non-perturbatively as the sum over two-particle irreducible (skeleton) diagrams. Each iteration of BDMC updates the bold propagators and self-energies until convergence, conserving physical constraints and partial resummations (Pollet et al., 2010, Kulagin et al., 2012, Davody, 2013, Mishchenko et al., 2014).

This approach enables efficient, sign-problem resistant computations in regimes characterized by strong local (on-site) correlations (cf. the Non-Crossing Approximation embedding for quantum impurity models (Gull et al., 2010)), and admits tight integration with dynamical mean-field theory (DMFT) for rapid convergence in correlated lattice systems (Pollet et al., 2010, Carlström, 2023).

3. Monte Carlo Sampling Algorithms and Configuration Space

For any diagram class (bare, skeleton, dual boson, etc.), the DiagMC algorithm is defined by its configuration space and update rules:

• Diagram configuration: DiagMC treats each Feynman diagram as a configuration specified by expansion order $n$, topology (vertex connectivities and line assignments), internal integration variables (interaction times $\{\tau_i\}$, momenta $\{k_i\}$), and external labels (e.g., $(k,\tau)$ for Green’s functions).
• Update rules: Standard moves include diagram insertion/removal (changing $n$), reconnection of propagator lines (altering topology), shifting internal variables, and worm-type updates to probe Green's function and vertex sectors. Proposal and acceptance ratios are determined to respect detailed balance: $$A(D \to D') = \min\left\{ 1, \frac{W(D')}{W(D)} \frac{P_\text{prop}(D' \to D)}{P_\text{prop}(D \to D')} \right\}$$ Each sampled diagram weight $W(D)$ combines products of propagator and interaction amplitudes, combinatorial symmetry factors, and sign factors from fermion loops or vertex permutations.
• Convergence and “sign blessing”: For fermionic models, the expansion coefficients oscillate in sign—often suppressing cancellations and promoting rapid (exponential) convergence up to a finite radius dependent on coupling strengths and temperature (Kulagin et al., 2012, Davody, 2013, Mishchenko et al., 2014). Statistical errors are monitored as functions of $n_\mathrm{max}$, and extrapolation is performed via fits of the form $S_n \sim A e^{-\alpha n}$ to estimate plateau values.
• Self-consistency: In BDMC-type schemes, outer loops self-consistently update bold propagators from accumulated estimators for self-energy or vertex functions, iterating until convergence. Data structures such as Young diagrams and hash-based tables may be used to index topologies efficiently.

4. Applications to Quantum Many-Body and Field Theories

DiagMC is versatile and has been implemented across diverse contexts:

• Quantum Magnets: In frustrated spin systems (e.g., the triangular-lattice Heisenberg model (Kulagin et al., 2012)), BDMC yields accurate thermodynamic observables, revealing quantum-to-classical correspondence in spatial correlators and temperature mapping.
• Quantum Field Theory: BDMC methods applied to $\phi^4$ theory in three dimensions (Davody, 2013) confirm Wilson-Fisher infrared fixed points at strong coupling, avoiding the sign problem by sampling bold line and vertex functions.
• Electron-Phonon Systems: Bold sampling in the Holstein model (Mishchenko et al., 2014) enables non-perturbative computation of polaron effective mass and residue, with vertex corrections obtained to high order (typically $N=4$–$6$).
• Molecular Electronic Structure: DiagMC for Møller–Plesset perturbation theory in molecules (Bighin et al., 2022) encodes Hugenholtz diagrams via adjacency matrices, yielding accurate correlation energies at quadratic scaling in the basis set size.
• Impurity Models and Open Quantum Systems: General impurity models with continuous hybridization functions are tackled via recursive determinant schemes at polynomial cost, even at low temperature (Li et al., 2020). Real-time impurity dynamics with local Markovian dissipation are sampled using single-contour, thermofield-based hybridization expansions, wherein dissipation reduces the sign problem and extends accessible simulation times (Vanhoecke et al., 2023).
• Disordered Semiconductors: Imaginary-time DQMC samples transport observables with dynamic and static disorder, applying a generalized Wick’s theorem for Gaussian variables and maintaining computational cost independent of system size (Wang et al., 2022).

5. Extensions: Embedding, Resummation, and Hybrid Diagrammatics

Hybrid Monte Carlo approaches incorporate analytic partial resummations and embedding techniques to maximize efficiency and reach strongly correlated regimes:

• DMFT Embedding: BDMC sampling may be performed exclusively for the non-local corrections to the self-energy, with local (mean-field) processes handled deterministically via DMFT (Pollet et al., 2010, Carlström, 2023). This integration reduces diagram count by orders of magnitude, expands the feasible convergence domain, and achieves unbiased calculation of self-energy and density of states in regimes inaccessible to conventional BDMC.
• Resummation Techniques: For rapidly diverging series, analytic continuation and resummation methods (Borel, Padé, conformal mapping) are used to push convergence boundaries. However, the sign structure and diagrammatic cancellation intrinsic to bold sampling often render further resummation unnecessary in practice.
• Variance Reduction and Histogram Techniques: Techniques such as flat histogram weighting (multicanonical, Wang–Landau) (Diamantis et al., 2013) and variance-reduction estimators (two-copy difference, binning, CDet-like resummations (Bighin et al., 2022)) are employed to suppress critical slowing down and enhance accuracy, especially in observables sensitive to long-time or high-order tails.

6. Computational Complexity and Scaling

DiagMC performance and scaling are determined by:

• Diagram order and topology: The factorial growth with expansion order $n$ is mitigated by stochastic sampling and, in bold or embedded schemes, by partially analytical resummations. In local or weakly correlated systems, memory and CPU usage may scale linearly or quadratically with system size.
• Basis Set and Bath Discretization: Continuous representation of hybridization functions obviates the need for bath discretization, enabling polynomial scaling in temperature and basis set size (Li et al., 2020, Bighin et al., 2022). For models with many orbitals, expansion order and variance may rise, but the approach remains systemically unbiased.
• Sign Problem: Oscillating diagram weights may suppress sign cancellations (“sign blessing”), promoting convergence, particularly in super-renormalizable theories and spin models (Kulagin et al., 2012, Davody, 2013). In open quantum systems, local dissipation further reduces the sign problem, enabling longer simulations in real time (Vanhoecke et al., 2023).

7. Physical Insights and Limitations

DiagMC methods have led to several physical insights unattainable by other techniques:

• Quantum-to-Classical Mapping: In frustrated quantum magnets, DiagMC reveals a direct mapping between quantum and classical spatial correlators, suggesting spin-liquid states absent magnetic order at low temperature (Kulagin et al., 2012).
• Critical Phenomena in Field Theories: Sampling bold Schwinger–Dyson equations for $\phi^4_3$ demonstrates direct convergence to Wilson-Fisher fixed points and dimensionless renormalized couplings, without recourse to analytic resummation (Davody, 2013).
• Transition Points and Emergent Phases: Order-by-order analysis of diagrammatic series enables estimation of phase boundaries, e.g., the charge density wave transition point in the extended Hubbard model via the dual boson DiagMC technique (Vandelli et al., 2020).

Limitations persist, particularly in regimes where the expansion parameter lies close to singularities, in strongly correlated metallic phases, or when high-order impurity vertices are computationally costly to obtain. Series divergence or severe sign problems may necessitate further theoretical advances in resummation, real-time contour methods, or integration with tensor network and quantum computing techniques. Nevertheless, the diagrammatic Monte Carlo method remains a foundational approach for unbiased, high-precision computation in many-body quantum physics and chemistry.