SciPost Submission Page
Automated evaluation of imaginary time strong coupling diagrams by sum-of-exponentials hybridization fitting
by Zhen Huang, Denis Golež, Hugo U. R. Strand, Jason Kaye
This is not the latest submitted version.
This Submission thread is now published as
Submission summary
| Authors (as registered SciPost users): | Zhen Huang · Jason Kaye |
| Submission information | |
|---|---|
| Preprint Link: | https://arxiv.org/abs/2503.19727v3 (pdf) |
| Date submitted: | Aug. 18, 2025, 7:14 p.m. |
| Submitted by: | Zhen Huang |
| Submitted to: | SciPost Physics |
| Ontological classification | |
|---|---|
| Academic field: | Physics |
| Specialties: |
|
| Approach: | Computational |
Abstract
We present an efficient separation of variables algorithm for the evaluation of imaginary time Feynman diagrams appearing in the bold pseudo-particle strong coupling expansion of the Anderson impurity model. The algorithm uses a fitting method based on AAA rational approximation and numerical optimization to obtain a sum-of-exponentials expansion of the hybridization function, which is then used to decompose the diagrams. A diagrammatic formulation of the algorithm leads to an automated procedure for diagrams of arbitrary order and topology. We also present methods of stabilizing the self-consistent solution of the pseudo-particle Dyson equation. The result is a low-cost and high-order accurate impurity solver for quantum embedding methods using general multi-orbital hybridization functions at low temperatures, appropriate for low-to-intermediate expansion orders. In addition to other benchmark examples, we use our solver to perform a dynamical mean-field theory study of a minimal model of the strongly correlated compound Ca$_2$RuO$_4$, describing the anti-ferromagnetic transition and the in- and out-of-plane anisotropy induced by spin-orbit coupling.
Author indications on fulfilling journal expectations
- Provide a novel and synergetic link between different research areas.
- Open a new pathway in an existing or a new research direction, with clear potential for multi-pronged follow-up work
- Detail a groundbreaking theoretical/experimental/computational discovery
- Present a breakthrough on a previously-identified and long-standing research stumbling block
Author comments upon resubmission
Dear Editor,
Please find enclosed our resubmission to SciPost Physics. We have carefully addressed all concerns raised by the reviewers and made the corresponding modification to the manuscript.
We sincerely thank you for your time and consideration.
Best regards,
Zhen Huang, Denis Golež, Hugo U. R. Strand, Jason Kaye
Please find enclosed our resubmission to SciPost Physics. We have carefully addressed all concerns raised by the reviewers and made the corresponding modification to the manuscript.
We sincerely thank you for your time and consideration.
Best regards,
Zhen Huang, Denis Golež, Hugo U. R. Strand, Jason Kaye
List of changes
- Page 1 [modified]: We refer to higher-order expansions of this type as bold pseudo-particle strong coupling, or hybridization, expansions.
- Page 2 [added]: Another proposed approach to mitigating the sign problem, for equilibrium real time calculations, uses the bold strong coupling expansion in conjunction with diagrammatic Monte Carlo \cite{haule23}.
- Page 2 [modified]: Several works have used explicit expressions of the propagator along with symbolic computation, as in algorithmic Matsubara integration \cite{taheridehkordi19, elazab22, burke25, assi24}, or analytical evaluation of imaginary time integrals combined with diagrammatic Monte Carlo \cite{vucicevic20,PhysRevResearch.3.023082,kovacevic25,tupitsyn21}.
- Page 3 [modified]: We define the kernel \begin{equation} \label{eq:kdef} K(\tau,\omega) = -\frac{e^{-\omega \tau}}{1 + e^{-\beta \omega}}\end{equation} for $\tau \in (0,\beta)$, and by the anti-periodicity property $K(\tau,\omega) = -K(\beta + \tau, \omega) = -K(-\tau, -\omega)$ for $\tau \in (-\beta, 0)$. We will represent the hybridization function as \begin{equation} \label{eq:deltasoe} \Delta_{\nu \lambda} (\tau) \approx \sum_{l=1}^p \Delta_{\nu \lambda l} K(\tau, \omega_l). \end{equation}
- Page 3 [modified]: Sec.~\ref{sec:hybfit} presents an algorithm to obtain a compact SOE approximation \eqref{eq:deltasoe}, with illustrative examples for several hybridization functions. Sec.~\ref{sec:algorithm} describes the diagram evaluation algorithm in detail. Finally, Sec.~\ref{sec:results} presents several benchmark examples, including timing results.
- Page 4 [modified]: \begin{equation} \label{eq:barycentric} r^{(k)}(z) = \frac{n(z)}{d(z)}=\sum_{j=1}^k \frac{w_j f_j}{z-z_j} /\sum_{j=1}\DIFdelbegin \DIFdel{^m }\DIFdelend \DIFaddbegin \DIFadd{^k }\DIFaddend \frac{w_j}{z-z_j}.\end{equation}
- Page 6 added: sum of $\delta$-functions (left), semicircle (middle), and sum of Gaussians (right)
- Page 8 [modified]: multiplying by $-1$, and swapping the $\omega_l \leq 0$ and $\omega_l > 0$ summation indices to avoid overflow.
- Page 12[added]: The crossing of the blue and orange lines in Fig.~\ref{fig:dimer_time}(b) is a consequence of the different expansion order convergence rates of the high temperature calculation (at $\beta t = 2$) and the lower temperature calculations (e.g., $\beta t = 16$). We attribute the difference to the presence of thermal excitations across the model's spectral gap (of size $\sim t$) when $\beta = 2t$. At the lower temperatures, these excitations are exponentially suppressed.
- Page 15 [added]: These phases are particularly well-suited to our method, as the gapped spectrum of the hybridization function leads to exponential convergence with diagrammatic order, which is even enhanced at lower temperatures. In contrast, the correlated metallic phase, characterized by a continuous hybridization spectrum, poses greater challenges, requiring significantly higher expansion orders to achieve convergence.
- Page 15 [added]: This approach would enable a more thorough exploration of how the expansion truncation error behaves in systems with many orbitals, e.g. full d or f orbital shells. Extension of our method to real time diagrammatics is also possible \cite{paprotzki25,kim2025},
- The lines of figure 1-6 have been adjusted to be thinner.
Current status:
Has been resubmitted
Reports on this Submission
Report #3 by Nayuta Takemori (Referee 1) on 2025-9-3 (Invited Report)
- Cite as: Nayuta Takemori, Report on arXiv:2503.19727v3, delivered 2025-09-03, doi: 10.21468/SciPost.Report.11856
Report
I thank the authors for their detailed replies in the revised manuscript. Most of my previous concerns have been addressed satisfactorily. I only have one remaining point (Major point 3) :
Regarding point 3:
In my previous comment, I referred to “low-frequency errors” and “self-energy”, which may have caused a misunderstanding of the authors. What I intended to point out is based on the behaviour visible in the τ-space Green’s functions shown in Fig. 5 and Fig. 7. Since these functions are obtained from Matsubara frequency data via the sum-of-exponentials representation, my concern is that the deviations near the endpoints (τ \sim 0 and τ \sim β), may originate from inaccuracies in the low-frequency Matsubara components of the hybridization fitting especially for the magnetic ordered state. I understand that there should be no significant problem when the system is gapped and exponential convergence is observed, but I appreciate if the authors can comment on the situation in more general cases (e.g., metallic case)?
Regarding point 3:
In my previous comment, I referred to “low-frequency errors” and “self-energy”, which may have caused a misunderstanding of the authors. What I intended to point out is based on the behaviour visible in the τ-space Green’s functions shown in Fig. 5 and Fig. 7. Since these functions are obtained from Matsubara frequency data via the sum-of-exponentials representation, my concern is that the deviations near the endpoints (τ \sim 0 and τ \sim β), may originate from inaccuracies in the low-frequency Matsubara components of the hybridization fitting especially for the magnetic ordered state. I understand that there should be no significant problem when the system is gapped and exponential convergence is observed, but I appreciate if the authors can comment on the situation in more general cases (e.g., metallic case)?
Recommendation
Ask for minor revision
Report
I am fully satisfied with the text. However I still think that is would be much better to publish the paper along with the code. This will not cause a big delay if the code is realised this fall.
Recommendation
Publish (meets expectations and criteria for this Journal)
Category:
remark
answer to question
Thanks for the suggestion!
We have made the code public, see github.
The code has a detailed documentation, which includes installation guide, user-friendly tutorials and reference manual.
Report #1 by Jan von Delft (Referee 3) on 2025-8-20 (Invited Report)
Report
I am satisfied with the new version of this paper and recommend publication at this time, without waiting for the announced code release this fall.
Recommendation
Publish (easily meets expectations and criteria for this Journal; among top 50%)
Author: Zhen Huang on 2025-09-19 [id 5836]
(in reply to Report 3 by Nayuta Takemori on 2025-09-03)We thank the referee for the clarification, but still have some confusion about the updated remarks.
Fig. 7 shows physical observables (magnetizations and fillings) at different temperatures, rather than the full τ-space Green's functions. While the data is computed from Green's function values at the $\tau$ endpoints, this figure corresponds to a numerical experiment for a system beyond previous approaches, and there is no exact reference for comparison. Therefore, it does not provide a direct indication of error.
In Fig. 5, the errors near the endpoints are not particularly large. In particular, for βt = 2, 16, and 1024, the maximum errors occur in the interior of [0, β]. According to Fourier theory, inaccuracies in the low-frequency Matsubara components would produce uniform errors for all $\tau$ in $[0, \beta]$, which is not observed in our data. Consequently, the deviations at the $\tau$ endpoints are not due to low-frequency fitting errors. Rather, the errors for a calculation of a given diagram order stem from truncation of higher-order diagrams, and as we have demonstrated for both gapped and metallic cases, these errors can be systematically reduced by including higher-order contributions.