NRG and CT-QMC Impurity Solvers: Validity and Error

Numerical renormalization group (NRG) and continuous-time quantum Monte Carlo (CT-QMC) answer complementary questions about an impurity model. NRG follows the flow across exponentially separated energy scales and supplies real-frequency transitions, but it approximates the bath, truncates states, and usually broadens a discrete spectrum. CT-QMC samples a continuous-time diagrammatic expansion without Trotter error, but it returns noisy imaginary-time or Matsubara data and may suffer from autocorrelation and a sign problem. Analytic continuation adds a separate inverse-problem uncertainty. A credible comparison therefore freezes one action and one observable, exposes each solver’s controls, and meets first on quantities that require neither NRG broadening nor Monte Carlo continuation.

Required background. The Anderson model supplies a shared benchmark. Convergence certification and covariance and resampling supply numerical error discipline.

Helpful background. Spectral moments supplies exact high-frequency checks.

Freeze the benchmark before choosing a solver

Use the spin-degenerate Anderson model and specify the bath through the same hybridization function $\Delta(z)$ in both codes. A concrete interacting benchmark is

$$D=1,\qquad \Gamma=0.1D,\qquad U=0.4D,\qquad \varepsilon_d=-\frac{U}{2},\qquad \beta D=100,$$

with a flat, particle–hole-symmetric bath of half-bandwidth $D$, zero field, and energies measured from the chemical potential. The bath-edge convention and the definition of $\Gamma=-\operatorname{Im}\Delta^R(0)$ are part of the input, not implementation details.

This benchmark has exact checks before either solver runs:

• particle–hole and spin symmetry give $n_d=1$ and $n_{d\sigma}=1/2$;
• fermionic anticommutation gives $G_\sigma(0^-)-G_\sigma(0^+)=1$;
• the leading high-frequency terms are

$$G_\sigma(i\omega_n) =\frac{1}{i\omega_n} +\frac{\varepsilon_d+U\langle n_{d\bar\sigma}\rangle}{(i\omega_n)^2} +O(\omega_n^{-3}) =\frac{1}{i\omega_n}+O(\omega_n^{-3});$$

• at $U=0$, each solver must reproduce the exact resonant-level result for the bath representation it actually uses,

$$G_{0\sigma}(i\omega_n) =\frac{1}{i\omega_n-\varepsilon_d-\Delta(i\omega_n)}$$

where $\Delta$ is that representation’s hybridization function. CT-QMC checks the specified continuous-bath target directly. For NRG, each $(\Lambda,z)$ calculation first checks the corresponding discretized $\Delta_{\Lambda z}$; only the controlled $\Lambda\to1^+$ and $z$-averaged limit is compared with the continuous-bath result.

The interacting comparison below uses the local susceptibility

$$\chi_{\mathrm{loc}} =(g\mu_{\mathrm B})^2\int_0^\beta\mathrm d\tau\, \langle S_z(\tau)S_z(0)\rangle,$$

at zero field. This is not the thermodynamic impurity susceptibility $\chi_{\mathrm{full}}-\chi_{\mathrm{bath}}$. Comparing one definition from NRG with the other from CT-QMC would create a solver discrepancy before any numerics begin.

How NRG reaches exponentially small scales

For a discretization parameter $\Lambda>1$ and shift $0<z\le1$, NRG groups the positive and negative bath energies into intervals whose edges scale as $D\Lambda^{-(n+z)}$. A unitary transformation maps one representative mode from each interval to a Wilson chain with asymptotic hopping

$$t_n\propto D\Lambda^{-n/2},$$

so adding site $N$ resolves a characteristic scale

$$\omega_N\sim D\Lambda^{-(N-1)/2}.$$

The Hamiltonian is diagonalized iteratively. After each site is added, high rescaled-energy multiplets are discarded so the next Hilbert space remains finite. Those discarded states are not simply thrown away when observables are evaluated. Each discarded shell state is tensored with the untouched states of the later Wilson-chain “environment”; the union over shells forms the complete approximate many-body basis. The full-density-matrix construction then assigns the appropriate thermal and environment weight to each shell. Anders and Schiller 2005, pp. 196801-1–196801-4 supplies the complete-basis construction; Weichselbaum and von Delft 2007, Eqs. (5)–(13) develops the full density matrix.

An NRG convergence grid changes one control at a time:

1. vary $\Lambda\to1^+$ and increase the number of interleaved $z$ shifts;
2. depending on the implementation, increase the number of kept multiplets, raise the rescaled energy cutoff, or tighten a separately defined discarded-weight tolerance;
3. exploit exact symmetries without accidentally removing perturbations under study;
4. move the temperature relative to the shell scale and verify known fixed-point spectra; and
5. compare exact sum rules and moments before applying a spectral kernel.

The discarded-weight diagnostic estimates whether important density-matrix support lies near the truncation boundary; it is a convergence indicator, not a rigorous error bound. Weichselbaum 2011, §§II–III defines this diagnostic.

For a fermionic operator $d_\sigma$, complete-basis transitions give discrete Lehmann weights. Evaluate $G(i\omega_n)$ directly from those unbroadened poles,

$$G_\sigma(i\omega_n) =\sum_{mn}\frac{(P_m+P_n) |\langle m|d_\sigma|n\rangle|^2} {i\omega_n+E_m-E_n},$$

Here $P_m$ denotes the appropriate full-density-matrix probability, including its shell and environment weight. Broadening is a later visualization or reconstruction step. A log-Gaussian width, its crossover to a linear kernel near $\omega=0$, and $z$ averaging can change a peak height or width without changing the underlying poles. Integrated weights and high-frequency moments should converge before pointwise line shapes are trusted. Bulla, Costi, and Pruschke 2008, §§II–IV reviews the discretization, truncation, thermodynamics, and spectral constructions.

What CT-HYB and CT-INT actually sample

CT-HYB expands the partition function in the impurity–bath hybridization. A configuration contains an ordered set of impurity creation and annihilation times; its weight is a local many-body trace times a determinant of hybridization functions connecting those times. CT-INT instead expands in interaction vertices, and its configuration weight contains a determinant of noninteracting propagators. The two algorithms are continuous in imaginary time, so there is no time-slice extrapolation, but their expansion orders, matrix conditioning, sign behavior, and efficient parameter regimes differ. Werner et al. 2006, Eqs. (2)–(5) introduces CT-HYB; Rubtsov, Savkin, and Lichtenstein 2005, §§II–III develops CT-INT.

Successive Monte Carlo configurations are correlated. For a scalar estimator $O$, define

$$\tau_{\mathrm{int}}=\frac12+\sum_{t=1}^{\infty}\rho_O(t), \qquad N_{\mathrm{eff}}\simeq\frac{N}{2\tau_{\mathrm{int}}}.$$

Blocking lengths must exceed the measured autocorrelation scale, and independent replicas should agree. A vector such as $G(\tau_i)$ needs its full covariance matrix: nearby time bins are not independent data points.

If weights have signs or phases, sampling $|w|$ gives

$$\langle O\rangle =\frac{\langle sO\rangle_{|w|}}{\langle s\rangle_{|w|}}, \qquad s=\frac{w}{|w|}.$$

The variance grows rapidly as the average sign decreases; more samples do not change that exponential complexity. The report therefore includes warmup, update acceptance and ergodicity tests, expansion-order distribution, matrix-stabilization checks, $\tau_{\mathrm{int}}$, binning, independent replicas, covariance, and average sign or phase. Gull et al. 2011, §§II, V, and IX compares the algorithm families and their error mechanisms.

Imaginary-axis information and continuation

A measured imaginary-time Green function is related to the local spectrum by

$$G(\tau)=-\int_{-\infty}^{\infty}\mathrm d\omega\, \frac{e^{-\tau\omega}}{1+e^{-\beta\omega}}\mathcal A(\omega).$$

Here $\mathcal A=-\pi^{-1}\operatorname{Im}G^R$ is the unit-normalized local density of states used throughout this impurity chapter; equivalently, the volume-wide spectral convention is $A=2\pi\mathcal A$. The kernel exponentially suppresses fine real-frequency structure, so many distinct spectra fit noisy $G(\tau)$ within its covariance. No increase in sample count turns analytic continuation into a well-conditioned inversion.

A continuation result inherits its default model or prior, regularizer, frequency grid, covariance treatment, and hyperparameter choice. Validate it by (i) forward transforming the proposed spectrum into $G(\tau)$ or $G(i\omega_n)$, including covariance; (ii) predicting held-out imaginary-axis points; (iii) testing synthetic spectra with comparable noise and resolution; and (iv) checking moments not imposed as priors. Normalizing $\int\mathcal A\,\mathrm d\omega=1$ during continuation is not an independent validation of that normalization. Jarrell and Gubernatis 1996, §§2–4 explains the inverse problem and maximum-entropy regularization.

The next figure separates raw solver outputs from reconstruction-dependent spectra.

Cross-solver validation begins before either real-frequency reconstruction. NRG supplies unbroadened complete-basis poles and CT-QMC supplies imaginary-axis estimates with their covariance; they first meet at definition-matched thermodynamics, $G(i\omega_n)$, and exact moments. Broadened NRG and analytically continued CT-QMC spectra are compared only after their separate reconstruction controls pass. Original workflow schematic; not to scale.

The same relationships in a text-first form are:

Stage	NRG lane	CT-HYB lane	Comparison or error rule
Frozen input	represent the declared $\Delta(z)$ by $(\Lambda,z)$ intervals and a Wilson chain	expand the same continuous-bath action in hybridization events	hold $U,\varepsilon_d,\mu,T,D$, spin/orbital normalization, and observable definitions fixed
Numerical controls	vary $\Lambda$, interleaved $z$ shifts, chain length, kept multiplets or energy cutoff, and symmetries	test warmup, acceptance, ergodicity, matrix stabilization, $\tau_{\mathrm{int}}$, binning, replicas, covariance, and average sign	keep discretization/truncation envelopes distinct from blocked sampling errors and sign diagnostics
Raw output	full-density-matrix thermodynamics and unbroadened complete-basis poles	$G(\tau)$ or $G(i\omega_n)$ with covariance and local estimators with blocked errors	retain the native representation and its correlations
First merge	direct pole sums or transforms	direct imaginary-axis and equal-time estimators	compare $n_d$, double occupancy, definition-matched $\chi_{\mathrm{loc}}$, $G(i\omega_n)$, identities, moments, and sum rules with no broadening or continuation
Optional real-frequency reconstruction	vary log-Gaussian/linear kernels separately from $\Lambda,z$, and truncation	vary prior or default model, covariance treatment, grid, regularizer, and hyperparameter	forward-test each reconstruction, then retain separate broadening and continuation envelopes
Real-frequency comparison	broadened NRG $\mathcal A(\omega)$	continued CT-HYB $\mathcal A(\omega)$	compare only after the separate controls pass; never collapse both reconstruction envelopes into one unexplained error bar

A definition-matched cross-solver benchmark

For the frozen particle–hole-symmetric benchmark, compare quantities that require no continuation or broadening first:

For the unit-normalized local spectrum, define $M_j=\int_{-\infty}^{\infty}\mathrm d\omega\,\omega^j\mathcal A(\omega)$. Then $M_0=1$ per spin, and particle–hole symmetry makes $M_1=0$ for this benchmark.

Observable	NRG check	CT-QMC check	Shared convention
Occupancy $n_d$	full-density-matrix trace, stable under $\Lambda,z$, and truncation	direct/equal-time estimator with covariance	spin sum; exact value $1$ here
Double occupancy $\langle n_\uparrow n_\downarrow\rangle$	unbroadened operator expectation	local estimator with blocked error	same $U,\beta$, and action sign
$\chi_{\mathrm{loc}}$	the same imaginary-time local-spin correlator or its Lehmann integral	integrated local-spin correlator with covariance	field couples only to $S_z$; not $\chi_{\mathrm{imp}}$
$G(i\omega_n)$	direct transform of unbroadened complete-basis poles	direct Matsubara estimator	same Fourier sign and frequency grid
Spectral weight and first moment	pole sums and large-$i\omega_n$ coefficients	equal-time discontinuity and large-$i\omega_n$ fit	$M_0=1$ per spin and $M_1=0$ here

For an NRG observable $X$, report separate envelopes from $(\Lambda,z)$, truncation, and—where relevant—thermal-shell placement. For a CT-QMC scalar, report the blocked statistical error, replica spread, sign, and any estimator or tail-fit variation. For $G(i\omega_n)$ use the CT-QMC covariance rather than treating frequencies as independent. If a combined uncertainty is quoted, state whether independent components were added in quadrature or conservative envelopes were added linearly; a truncation envelope is not automatically a Gaussian standard deviation.

A useful discrepancy diagnostic for a scalar is

$$R_X=\frac{|X_{\mathrm{NRG}}-X_{\mathrm{QMC}}|} {u_{\mathrm{NRG}}+u_{\mathrm{QMC}}}.$$

$R_X\le1$ means the declared conservative intervals overlap; it does not certify that both methods are unbiased. For a Matsubara vector, inspect the residual with the CT-QMC covariance and repeat the NRG control grid. Only after $n_d$, double occupancy, $\chi_{\mathrm{loc}}$, $G(i\omega_n)$, and the moments agree should one compare a real-frequency Kondo width. That last comparison must retain an NRG broadening envelope and a CT-QMC continuation ensemble separately.

The chapter structure diagram places both solvers after the model and its candidate fixed-point classification but before observable inference: numerical flow must test the candidate branch rather than assume it.

NRG or CT-QMC must verify the candidate infrared branch and quantify solver-specific errors before observable inference. The one-orbital Anderson route supplies one active hybridization eigenchannel; multichannel branches require separate microscopic input. Original workflow schematic, not to scale.

The impurity claim test matrix states the solver-specific stopping rules.

Common pitfalls

Lowering an NRG energy cutoff to “keep more.” In the common rescaled-energy convention, increasing the cutoff keeps more states. Only a separately defined discarded-weight tolerance is reduced.

Comparing different susceptibilities. $\chi_{\mathrm{loc}}$ and $\chi_{\mathrm{imp}}=\chi_{\mathrm{full}}-\chi_{\mathrm{bath}}$ answer different questions and can differ substantially.

Transforming a broadened NRG plot to Matsubara frequency. Use the unbroadened complete-basis poles. Broadening is not part of the Hamiltonian and should not contaminate a continuation-free comparison.

Using a continuation constraint as its own check. A unit-normalized spectrum proves nothing about $M_0$ if unit normalization was imposed. Check the equal-time discontinuity or the leading Matsubara coefficient independently.

Exercises

Estimate the required Wilson-chain length. With $D=1\,\mathrm{eV}$ and $\Lambda=2$, approximately which iteration reaches $\omega_N=10^{-6}\,\mathrm{eV}$?

Solution

Solve $2^{-(N-1)/2}=10^{-6}$. Since $\log_2(10^6)\simeq19.93$,

$$N\simeq1+2(19.93)\simeq41.$$

This is only a scale estimate. The usable resolution also depends on $z$ averaging, truncation, temperature, and the observable.

Convert autocorrelation into effective samples. A CT-HYB run records $N=2.0\times10^6$ measurements with $\tau_{\mathrm{int}}=25$. Estimate $N_{\mathrm{eff}}$. If the average sign is $0.20$, what rough sign-reweighted scale warns about the variance?

Solution

The autocorrelation estimate is

$$N_{\mathrm{eff}}\simeq\frac{2.0\times10^6}{2(25)}=4.0\times10^4.$$

Sign reweighting often inflates variance roughly as $\langle s\rangle^{-2}$. The corresponding warning scale is $N_{\mathrm{eff}}\langle s\rangle^2\simeq1.6\times10^3$. This is not a universal error formula, but it shows why the raw measurement count is misleading.

Check the symmetric high-frequency moment. Show why the coefficient of $(i\omega_n)^{-2}$ vanishes for $\varepsilon_d=-U/2$ at spin symmetry.

Solution

Spin and particle–hole symmetry give $\langle n_{d\bar\sigma}\rangle=1/2$. Hence

$$\varepsilon_d+U\langle n_{d\bar\sigma}\rangle =-\frac U2+\frac U2=0.$$

Both solvers must therefore approach $G(i\omega_n)=1/(i\omega_n)+O(\omega_n^{-3})$. A nonzero fitted first moment points to insufficient high-frequency range, a bath/action mismatch, broken symmetry, or an estimator error.

Diagnose a false agreement. NRG and CT-QMC give the same occupancy, but their real-frequency peak widths disagree by 30%. List the next variations.

Solution

For NRG, vary $\Lambda$, $z$ averaging, kept-state threshold, and broadening kernel separately, checking moments after each change. For CT-QMC, inspect covariance and autocorrelation, vary statistics, continuation method, default model, frequency grid, and synthetic-recovery tests. Compare $G(i\omega_n)$ directly before continuing. Occupancy agreement validates one integral but cannot determine a narrow spectral width.

References

• Anders, F. B., and Schiller, A. (2005). “Real-time dynamics in quantum-impurity systems: A time-dependent numerical renormalization-group approach.” Physical Review Letters 95, 196801. doi:10.1103/PhysRevLett.95.196801.
• Bulla, R., Costi, T. A., and Pruschke, T. (2008). “Numerical renormalization group method for quantum impurity systems.” Reviews of Modern Physics 80, 395–450. doi:10.1103/RevModPhys.80.395.
• Gull, E., Millis, A. J., Lichtenstein, A. I., Rubtsov, A. N., Troyer, M., and Werner, P. (2011). “Continuous-time Monte Carlo methods for quantum impurity models.” Reviews of Modern Physics 83, 349–404. doi:10.1103/RevModPhys.83.349.
• Jarrell, M., and Gubernatis, J. E. (1996). “Bayesian inference and the analytic continuation of imaginary-time quantum Monte Carlo data.” Physics Reports 269, 133–195. doi:10.1016/0370-1573(95)00074-7.
• Rubtsov, A. N., Savkin, V. V., and Lichtenstein, A. I. (2005). “Continuous-time quantum Monte Carlo method for fermions.” Physical Review B 72, 035122. doi:10.1103/PhysRevB.72.035122.
• Weichselbaum, A. (2011). “Discarded weight and entanglement spectra in the numerical renormalization group.” Physical Review B 84, 125130. doi:10.1103/PhysRevB.84.125130.
• Weichselbaum, A., and von Delft, J. (2007). “Sum-rule conserving spectral functions from the numerical renormalization group.” Physical Review Letters 99, 076402. doi:10.1103/PhysRevLett.99.076402.
• Werner, P., Comanac, A., de’ Medici, L., Troyer, M., and Millis, A. J. (2006). “Continuous-time solver for quantum impurity models.” Physical Review Letters 97, 076405. doi:10.1103/PhysRevLett.97.076405.