Accelerating a Strong-Coupling Non-Equilibrium Steady-State Impurity Solver using (Quantics) Tensor Trains

Abstract: Including higher order diagrammatic corrections to the strong-coupling expansion is mainly limited by the evaluation of high-dimensional, time-ordered integrals. In this work we present and compare four different parametrizations of the integrands in order to obtain a low-rank (quantics) tensor-train representation using tensor cross interpolation. Particular emphasis is placed on a quantics time-difference formulation in which the required retarded convolutions are performed directly in quantics tensor-train form. Using controlled Gaussian benchmarks, we analyze the accuracy, bond dimensions, and computational scaling of the different approaches. We then validate the most promising formulations in self-consistent equilibrium and nonequilibrium DMFT calculations and demonstrate calculations up to the third order in the strong-coupling expansion. Finally, we extend the solver to impurity models with retarded density-density interactions and apply it within nonequilibrium extended DMFT. Our results show that tensor cross interpolation substantially reduces the cost of evaluating higher-order diagrams and provides a controlled, systematically improvable framework for nonequilibrium quantum impurity calculations.

• The paper develops four tensor-train parametrizations with tensor cross interpolation and introduces direct quantics convolutions, reducing high-dimensional Keldysh integral evaluation from cubic to nearly linear or logarithmic grid scaling.
• The paper validates a third-order nonequilibrium steady-state solver in equilibrium and photodoped DMFT, showing improved quasiparticle features over OCA and practical runtimes from seconds to minutes on a laptop.
• The paper extends the method to retarded interactions in EDMFT, where OCA and TOA provide controlled results beyond NCA and enable more reliable simulations of dynamical screening and nonequilibrium GW+EDMFT systems.

The paper "Accelerating a Strong-Coupling Non-Equilibrium Steady-State Impurity Solver using (Quantics) Tensor Trains" (2608.13146) addresses the principal computational bottleneck of high-order strong-coupling impurity solvers: the evaluation of high-dimensional, time-ordered integrals over the Keldysh contour. The authors present four tensor-train (TT) parametrizations of these integrands, constructed via tensor cross interpolation (TCI), and systematically compare their accuracy, bond dimensions, and runtime scaling. Building on the most efficient formulations, they construct a third-order nonequilibrium steady-state (NESS) impurity solver, validate it in equilibrium and photodoped DMFT calculations on the Bethe lattice, and extend it to impurity models with retarded density-density interactions for use in nonequilibrium extended DMFT (EDMFT).

Formalism and computational setting

The work builds on the self-consistent strong-coupling (hybridization) expansion around the atomic limit, formulated in terms of pseudo-particle propagators $\mathcal{G}$ and a pseudo-particle self-energy $\Sigma$ satisfying a Dyson equation on the Keldysh contour. An $n$th-order diagram requires integration over $2n-1$ cyclically ordered contour times, which is the dominant cost at third order and beyond. For nonequilibrium steady states, time-translational invariance allows elimination of the imaginary-time branch, so that all objects are reconstructed from greater and lesser components depending only on time differences.

A central technical choice is the Keldysh parametrization of Ref. Eckstein (2024): each contour time is represented by a physical time plus a Keldysh index $\sigma_j \in \{\pm 1\}$. This introduces an additional sum over $2^D$ Keldysh configurations but yields substantially lower TT bond dimensions than the cyclic parametrization, an advantage that more than compensates for the extra cost.

Four decomposition strategies

The authors compare four quasi-linear schemes for evaluating the time-ordered integrals:

1. Time decomposition: TCI sampling of the integrand on an extended rectangular domain of absolute times, followed by forward/backward recursions implementing the nested ordered sums. Runtime scales as $\mathcal{O}(D N_t \chi^2)$.
2. Quantics time decomposition: the same masked integrand sampled in fused quantics representation; crucially, the Heaviside time-ordering mask admits an exact fused QTT with bond dimension two, so the ordering constraint factorizes. Evaluation costs $\mathcal{O}(N_t \log_2 N_t\, \chi^3)$.
3. Time difference decomposition: the integrand is expressed in successive time differences, exploiting steady-state time-translational invariance; convolutions are evaluated via FFT as in Ref. Eckstein (2024).
4. Quantics time difference decomposition (the key methodological contribution): the time-difference integrand is decomposed in fused quantics form, and retarded convolutions are performed directly in QTT format via an auxiliary carry bit, formulated as a rank-2 MPO enforcing binary addition locally. This increases the bond dimension by only a factor of two (Riemann) or at most four with block-diagonal structure (trapezoidal weights), compared to roughly an elevenfold increase for the quantics Fourier transform. Trapezoidal endpoint corrections are incorporated as a three-term TT sum exploiting block-diagonal bond structure.

Gaussian benchmarks

Controlled scalar benchmarks using Gaussian propagators isolate the TCI error from quadrature error. The time difference parametrization reaches machine precision at the smallest bond dimension ($\chi = 16$); the time decomposition needs $\chi = 48$ and quantics time $\Sigma$0, while both quantics difference variants saturate at higher errors ($\Sigma$1 fused, $\Sigma$2 interleaved) near $\Sigma$3. Notably, the required bond dimension is independent of grid size $\Sigma$4 for all methods, indicating stable low-rank structure. Going from OCA to TOA requires approximately a fourfold increase in bond dimension across all methods — a systematic trend confirmed later in the self-consistent calculations.

Runtime analysis on a laptop CPU shows that direct evaluation scales as $\Sigma$5, while TCI-based evaluation is nearly linear and dominated by sampling cost. Rook pivoting provides large speedups over full pivoting without accuracy loss (e.g., a speedup of 326.9 for TCI time at $\Sigma$6). Already at $\Sigma$7, the logarithmic sampling cost of the QTCI difference variants outweighs their larger bond dimensions. The benchmarks also demonstrate that trapezoidal endpoint corrections are essential: plain Riemann summation produces leading short-time errors that can violate the Hermiticity of the lesser self-energy and decay only as $\Sigma$8.

Equilibrium and nonequilibrium DMFT validation

In equilibrium DMFT for the half-filled Hubbard model on the Bethe lattice ($\Sigma$9, metallic regime where higher orders matter), the RMSE between successive hybridization functions decreases exponentially with iteration count and then saturates at a plateau controlled by solver accuracy — analogous to stochastic solvers. Separate scans over bond dimension $n$0 and grid size $n$1 show that decomposition and quadrature errors are controlled independently. All three tested parametrizations recover the converged OCA solution within fewer than about twenty DMFT iterations. A stringent diagnostic — recovery of the fluctuation-dissipation relation, i.e., the logarithmic ratio of lesser and greater spectral components approaching $n$2 — is satisfied by converged results, with Riemann quadrature showing systematic frequency-dependent deviations while trapezoidal results converge cleanly. Practically, a serial OCA Bethe iteration with fused QTCI at $n$3 takes seconds to a minute on a laptop, versus several minutes for non-quantics TCI at $n$4.

For nonequilibrium steady states, the solver is applied to photodoped Mott insulators using the fixed-distribution protocol of Küenzel et al. (2024). Comparing NCA, OCA, and TOA at fixed excitation density $n$5, the Mott gap remains stable, OCA and TOA agree well in the undoped Hubbard-band regions, and TOA still improves the holon/doublon quasiparticle peaks over OCA — evidence for overall convergence of the diagrammatic series, consistent with prior frequency-domain QTCI studies. Residual discrepancies with inchworm QMC may stem from statistical noise and the effective cutoff time in the Monte Carlo reference data.

Extension to retarded interactions and EDMFT

The solver is extended to impurity models with retarded density-density interactions via a double expansion: strong coupling in the fermionic hybridization combined with a weak-coupling expansion in the bosonic propagator $n$6, following Golez, Eckstein, and Werner (2015). Diagram counts grow rapidly (e.g., 27,648 Keldysh-enumerated TOA diagrams without symmetry reduction, reduced to 3,456 by spin and endpoint symmetries), underscoring why such calculations were previously inaccessible.

Square-lattice EDMFT results at $n$7, $n$8 show rapid fermionic convergence but more pronounced bosonic sensitivity to expansion order. In the photodoped case, NCA substantially underestimates the photodoped carrier density relative to OCA/TOA (e.g., $n$9 vs. $2n-1$0/$2n-1$1 at $2n-1$2), while OCA and TOA are nearly identical — providing a controlled estimate of the truncation error. Increasing $2n-1$3 progressively fills the Mott gap through dynamical screening, a trend NCA overestimates in stability of the insulating phase. Photodoping induces low-energy charge spectral weight in the bosonic spectrum, which is more order-sensitive than the fermionic response. Since the local screened interaction feeds directly into the GW contribution in GW+EDMFT, this controlled treatment of the bosonic channel is a prerequisite for quantitative nonequilibrium GW+EDMFT simulations.

Limitations and open questions

Several caveats are stated explicitly. The quantics difference decompositions saturate at finite TCI errors ($2n-1$4–$2n-1$5), attributed to pivot-search difficulties in sparse tensors; whether this can be overcome by improved ergodicity handling remains open. Higher-order quadrature weights were implemented only for the difference-based schemes, and interleaved quantics layouts suffer from ergodicity problems and larger bond dimensions despite better asymptotic sampling complexity. The steady-state Dyson equation requires either sufficiently fast decay of $2n-1$6 within the simulation window or an artificial pseudo-particle bath broadening $2n-1$7; the quasiparticle peak depends noticeably on $2n-1$8, and annealing to zero broadening was deliberately left unexplored. In EDMFT, convergence is demonstrated only up to third order at specific parameters, and the double-expansion diagram proliferation suggests that fourth-order calculations would require further algorithmic advances. Finally, the optimal choice between non-quantics and quantics schemes depends on the required time cutoff $2n-1$9, and a fully QTT-contained DMFT cycle avoiding dense evaluation of $\sigma_j \in \{\pm 1\}$0 on large grids is proposed but not implemented.

Conclusion

This work establishes a controlled, systematically improvable framework for nonequilibrium strong-coupling impurity solvers by combining TCI with a direct real-time QTT convolution algorithm. The comparison of four parametrizations shows that time-difference representations dominate, with the quantics variant offering logarithmic scaling in grid size and the non-quantics variant lower bond dimensions and superior quadrature flexibility. The resulting solver makes third-order nonequilibrium calculations routine — OCA iterations take seconds on a laptop — and delivers the first controlled truncation-error estimates for impurity models with retarded interactions, removing a long-standing obstacle to quantitative nonequilibrium GW+EDMFT.