# Dual Bethe-Salpeter equation for the multi-orbital lattice susceptibility
within dynamical mean-field theory

Erik G. C. P. van Loon  NanoLund and Division of Mathematical Physics,
Department of Physics, Lund University, Lund, Sweden  Hugo U. R. Strand
School of Science and Technology, Örebro University, SE-701 82 Örebro, Sweden
Institute for Molecules and Materials, Radboud University, 6525 AJ Nijmegen,
the Netherlands

###### Abstract

Dynamical mean-field theory describes the impact of strong local correlation
effects in many-electron systems. While the single-particle spectral function
is directly obtained within the formalism, two-particle susceptibilities can
also be obtained by solving the Bethe-Salpeter equation. The solution requires
handling infinite matrices in Matsubara frequency space. This is commonly
treated using a finite frequency cut-off, resulting in slow linear
convergence. We show that decomposing the two-particle response in local and
non-local contributions enables a reformulation of the Bethe-Salpeter equation
inspired by the dual boson formalism. The re-formulation has a drastically
improved cubic convergence with respect to the frequency cut-off, facilitating
the calculation of susceptibilities in multi-orbital systems considerably. The
dual Bethe-Salpeter equation uses the fully reducible vertex which is free
from vertex divergences. We benchmark the approach on several systems
including the spin susceptibility of strontium ruthenate Sr2RuO4, a strongly
correlated Hund’s metal with three active orbitals. We propose the dual Bethe-
Salpeter equation as a new standard for calculating two-particle response
within dynamical mean-field theory.

##  I Introduction

Electronic correlations are present in many interesting quantum materials,
including Mott insulators [1], ruthenates [2], iron pnictides [3] and cuprates
[4]. In these materials, the interaction between electronic quasiparticles is
strong, leading to the breakdown of the independent electron picture. The
theoretical description of this quantum many-body problem is hard and
generally requires approximations. The Dynamical mean-field theory (DMFT) [5,
6] approximation is one prominent example, which quantitatively explains many
experimentally observed phenomena [2], although there are also situations
where it is insufficient [7, 8].

The properties of a system of correlated electrons can be discussed in terms
of expectation values of the form $\langle
c^{\dagger}_{\alpha}c^{\phantom{\dagger}}_{\beta}\rangle$, $\langle
c^{\dagger}_{\alpha}c^{\phantom{\dagger}}_{\beta}c^{\dagger}_{\gamma}c^{\phantom{\dagger}}_{\delta}\rangle$
and so on, where $c_{\alpha}$, $c_{\alpha}^{\dagger}$ stand for the
annihilation and creation operator of an electron and $\alpha$ is a combined
time $\tau$, space $\mathbf{R}$, and spin-orbital $a$ label
$\alpha\equiv(a,\mathbf{R},\tau)$. Here $\langle
c^{\dagger}_{\alpha}c^{\phantom{\dagger}}_{\beta}\rangle$ is called a single-
particle quantity, $\langle
c^{\dagger}_{\alpha}c^{\phantom{\dagger}}_{\beta}c^{\dagger}_{\gamma}c^{\phantom{\dagger}}_{\delta}\rangle$
a two-particle quantity, and so on. With the increase in the number of
operators and labels in the correlation function, these many-particle
quantities become increasingly complicated. In practice, expectation values
with more than two pairs of operators are often (but not always [9, 10]) out
of reach, but both the single-particle and two-particle quantities are
experimentally accessible. A simplification of the structure is possible for
translational invariant systems in equilibrium where only relative time and
position matters. Whence, the number of space-time (or momentum-frequency)
labels in the correlation function can be reduced by one. Another name for the
two-particle quantitiy is the generalized susceptibilitiy $\chi$, see Fig. 1,
since these correlation functions describe the linear response of the
electronic system to an external field, according to the fluctuation-
dissipation theorem.

Figure 1: Generalized susceptibility $\chi$ (top), schematic scattering
experiment (bottom).

Experimentally, the susceptibilities are accessible via techniques such as
inelastic neutron scattering (INS) [11], electron energy loss spectroscopy
(EELS) [12] or resonant inelastic x-ray scattering (RIXS) [13]. Roughly
speaking, these experiments consist of shooting a probing particle with a
known energy $E_{\text{in}}$ and momentum $\mathbf{k}_{\text{in}}$ at the
correlated material of interest, and measuring the distribution $S$ over
energy $E_{\text{out}}$ and momentum $\mathbf{k}_{\text{out}}$ of the outgoing
probing particle. The energy and momentum difference between the ingoing and
outgoing particle is transferred to the sample, i.e. the energy
$\omega=E_{\text{in}}-E_{\text{out}}$ and momentum
$\mathbf{q}=\mathbf{k}_{\text{in}}-\mathbf{k}_{\text{out}}$ are absorbed by
the correlated material, and large magnitudes of $S$ indicate that the
material supports (collective) excitations at $(\omega,\mathbf{q})$. Thus, the
scattering amplitude $S$ measured by these spectroscopic techniques is
proportional to the susceptibility $\chi(\omega,\mathbf{q})$ of the material,
up to further scattering matrix element and interaction effects, see Fig. 1.

Thus, the experimental spectrum $S$ only depends on the transferred energy and
momentum $(\omega,\mathbf{q})$, while the generalized electronic
susceptibility $\chi$ is a four-index object, which depends on three
frequencies and momenta due to time-space translation symmetry. To get to the
experimentally relevant susceptibility, two of the frequencies and momenta of
the generalized susceptibility need to be traced out. To understand the
meaning of this trace, consider the momentum labels
$\mathbf{q},\mathbf{k},\mathbf{k}^{\prime}$ of the generalized susceptibility.
These labels correspond to an electronic transition
$\mathbf{k}\rightarrow\mathbf{k}+\mathbf{q}$ and another electronic transition
$\mathbf{k}^{\prime}+\mathbf{q}\rightarrow\mathbf{k}^{\prime}$, see Fig. 1. In
the theoretical description, it is possible to consider/calculate the
likelihood of such a transition as a function of both $\mathbf{k}$ and
$\mathbf{k}^{\prime}$, as well as $\mathbf{q}$, but in the experiment, only
the momentum transfer $\mathbf{q}$ is observable and there is no way to
determine the momenta of the electronic states involved in the transition, and
the trace takes the form

$S(E_{\text{out}},\mathbf{k}_{\text{out}})\propto\chi(\omega,\mathbf{q})=\sum_{\mathbf{k}\mathbf{k}^{\prime}}\sum_{\nu\nu^{\prime}}\chi(\omega,\mathbf{q};\nu,\mathbf{k};\nu^{\prime},\mathbf{k}^{\prime})\,.$

Thus, there is a tremendous loss of information when going from the full
generalized susceptibility to a particular traced-out version that is actually
experimentally observable.

In this work, we will show that this loss of information can actually be used
to our advantage when we want to calculate the susceptibility in dynamical
mean-field theory [6]. The usual formulation of the Bethe-Salpeter equation in
Dynamical Mean-Field Theory proceeds by first calculating the generalized
susceptibility of the material, based on the generalized susceptibility of the
auxiliary impurity model, and only taking the trace over momenta and fermionic
frequencies at the end. The disadvantage of this is that this formulation of
the Bethe-Salpeter equation converges only _linearly_ with respect to the
number of internal frequencies used in the calculation, especially for local
scattering processes. Some techniques have been developed to correct for the
known high-frequency asymptotics [14, 15, 16, 17, 18], but converging with
respect to the frequency grid remains the main challenge for reaching
experimentally relevant and physically interesting temperatures.

An alternative formulation, based on the dual boson approach [19], takes
advantage of the distinction between local and non-local fluctuations in
Dynamical Mean-Field Theory. By using the local single-frequency and two-
frequency susceptibilities of the auxiliary impurity model, it is possible to
efficiently take into account all local scattering processes, regardless of
their interal frequency. The resulting dual Bethe-Salpeter equation (DBSE)
then only has to account for non-local scattering processes [20, 21, 22, 23],
which converge _cubically_ with respect to the number of fermionic frequencies
used in the calculation. This is a substantial improvement that enables
calculations at lower temperature or in more complex systems, since the
computational effort scales with the number of orbitals and the momentum
discretization, as well as the number of fermionic frequencies.

Previous applications [24, 25] of this dual boson based approach were
restricted to single-orbital models. Here, we present the generalization to
multi-orbital systems. We provide a publicly available open-source solver
implementation of the dual Bethe-Salpeter equation in the two particle
response function toolbox (TPRF), based on the toolbox for interacting quantum
systems (TRIQS [26]).

##  II Theory

We will consider lattice models in thermodynamic equilibrium at finite
temperature, using the imaginary time/Matsubara formalism. We start with a
brief overview of the susceptibility and how it is usually calculated in DMFT,
including our notation for describing it. Some familiarity with these concepts
is assumed, we refer the interested reader to the literature [6, 27, 28, 7,
29] for more details.

###  II.1 Notation

The single-particle Green’s function is defined as
$G_{a\bar{b}}(\tau,\mathbf{R})=\left\langle
c^{\phantom{\dagger}}_{\tau,a,\mathbf{R}}c^{\dagger}_{0,b,\mathbf{0}}\right\rangle$,
where $\tau\in[0,\beta)$ is the imaginary time, $a$ and $b$ are combined spin-
orbital labels and $\mathbf{R}$ is a lattice vector. We use $\bar{\cdot}$ to
denote indices that appear on creation operators, which becomes very relevant
for many-particle Green’s functions. The Fourier transform of
$G_{a\bar{b}}(\tau,\mathbf{R})$ is $G_{a\bar{b}}(\nu,\mathbf{k})$, where $\nu$
is a fermionic Matsubara frequency and $\mathbf{k}$ is the momentum in the
first Brillouin Zone.

Similarly, the generalized susceptibility
$\chi_{\bar{a}b\bar{c}d}(\omega,\mathbf{q})$ is the Fourier transform of

$\displaystyle\chi_{\bar{a}b\bar{c}d}(\tau,\mathbf{R})=$
$\displaystyle\left\langle
c^{\dagger}_{\tau,\mathbf{R},a}c^{\phantom{\dagger}}_{\tau,\mathbf{R},b}c^{\dagger}_{0,0,c}c^{\phantom{\dagger}}_{0,0,d}\right\rangle$
$\displaystyle-\left\langle
c^{\dagger}_{\tau,\mathbf{R},a}c^{\phantom{\dagger}}_{\tau,\mathbf{R},b}\right\rangle\left\langle
c^{\dagger}_{0,0,c}c^{\phantom{\dagger}}_{0,0,d}\right\rangle$
$\displaystyle=$ $\displaystyle\left\langle
c^{\dagger}_{\tau,\mathbf{R},a}c^{\phantom{\dagger}}_{\tau,\mathbf{R},b}c^{\dagger}_{0,0,c}c^{\phantom{\dagger}}_{0,0,d}\right\rangle$
$\displaystyle-
G_{b\bar{a}}(\beta,\mathbf{0})G_{d\bar{c}}(\beta,\mathbf{0})\,.$ (1)

Here, $\omega$ is a bosonic Matsubara frequency. For clarity, we use
$\mathbf{q}$ for the bosonic momentum and $\mathbf{k}$ for fermionic momenta,
both take values in the first Brillouin Zone.

Note that the order of barred operators differs between the single-particle
Green’s function and the susceptibility. The notation of the single-particle
Green’s function is conventional. For the susceptibility there are many
notations found in the literature, we follow the notation of TPRF [30], see
also Appendix A.

The spin-orbital labels can be contracted to go from the generalized
susceptibility to a particular physical susceptibility,
$\chi^{AB}(\omega,\mathbf{q})=\sum_{abcd}A_{ab}B_{cd}\chi_{\bar{a}b\bar{c}d}(\omega,\mathbf{q})$,
where $A$ and $B$ are matrices in spin-orbital space. As an example, the
$S^{z}$ component of the spin susceptibility is obtained as
$\chi^{S^{z}S^{z}}(\omega,\mathbf{q})=\sum_{abcd}\hat{\sigma}^{z}_{ab}\hat{\sigma}^{z}_{cd}\chi_{\bar{a}b\bar{c}d}(\omega,\mathbf{q})$,
where $\hat{\sigma}^{z}$ is the appropriate Pauli matrix.

The static susceptibility $\chi^{AB}(\omega=0,\mathbf{q}=0)$ can be
interpreted in linear response theory. Define the operators
$\hat{A}=\sum_{ab}A_{ab}c^{\dagger}_{a}c_{b}$ and
$\hat{B}=\sum_{cd}B_{cd}c^{\dagger}_{c}c_{d}$, which we assume to be
Hermitian. Then, the application of a spatially homogeneous field
$\lambda\hat{B}$, in the sense that the Hamiltonian changes as
$\hat{H}\rightarrow\hat{H}-\lambda\hat{B}$, leads to a change in
$\langle\hat{A}\rangle$, where the zero field response is
$d\langle\hat{A}\rangle/d\lambda\big{|}_{\lambda=0}=\chi^{AB}(\omega=0,\mathbf{q}=0)$.
This definition can be generalized to spatially inhomogeneous fields (finite
$\mathbf{q}$) by considering enlarged unit cells [31], while finite $\omega$
corresponds to time dependent applied fields.

###  II.2 Bethe-Salpeter equation in Dynamical Mean-Field Theory

In dynamical mean-field theory, the susceptibility can be calculated using the
Bethe-Salpeter equation [6]. Here, we introduce the equations without orbital,
frequency and momentum labels to emphasize the structure of the equations,
explicit formulas including labels are given later. To distinguish the
physical susceptibility $\chi$ of the system from the auxiliary impuritiy
quantities in DMFT, it is commonly called the _lattice_ suceptibility. In
general, the Bethe-Salpeter equation for $\chi$ takes the form

$\chi=\chi^{0}+\chi^{0}\Gamma\chi,$ (2)

where the bubble $\chi^{0}=G\ast G$ denotes the independent propagation of a
particle-hole pair according to the Green’s function $G$ given by the DMFT
approximation. $\Gamma$ is the irreducible vertex, which is assumed to be
local in the DMFT approximation. Both $G$ and $\Gamma$ can be extracted from
the solution of DMFT’s auxiliary impurity problem. The Bethe-Salpeter Equation
(2) is a self-consistent equation for $\chi$, which can be inverted to give
$\chi$ as a function $\chi^{0}$ and $\Gamma$.

In DMFT, it is also possible to discuss the susceptibility of the auxiliary
impurity problem, $X$, which has it’s own Bethe-Salpeter equation

$X=X^{0}+X^{0}\Gamma X.$ (3)

Here, the impurity bubble $X^{0}=g\ast g$ is defined in analogy with the
lattice bubble, and $g$ is the Green’s function of the impurity model.
Importantly, the same vertex $\Gamma$ is used for both BSEs in the DMFT
approximation.

###  II.3 Locality

Locality is a central concept in DMFT. The DMFT self-consistency condition
guarantees that $g$ is the local part of $G$, which is often expressed as
$g=G(\mathbf{R}=0)=\frac{1}{N_{k}}\sum_{\mathbf{k}}G(\mathbf{k})$. This
implies that $X^{0}$ is the local part of $\chi^{0}$ as well, since the bubble
is calculated as a product in real space. However, $X$ is not the local part
of $\chi$ [22], since the Bethe-Salpeter equation (2) allows for processes
where the electron-hole pair moves to a different site and back. Physically,
$X$ corresponds to the linear response of the impurity model while keeping the
dynamical mean-field fixed, whereas $\chi$ corresponds to linear response
where the mean-field is self-consistently adjusted.

Although $X$ is not identical to (the local part of) $\chi$, it might still be
a much better starting point for the calculation of $\chi$ than $\chi^{0}$,
since $X$ already contains local vertex corrections in the form of Eq. (3). We
should stress that $X$ can be extracted directly from the impurity model,
without needing to use the identity Eq. (3). Thus, the central idea of the
dual Bethe-Salpeter equation [19, 20] is to express $\chi$ as $X$ and some
corrections. Finally, we note that this idea has also been used to decompose
the polarization [32].

###  II.4 Dual Bethe-Salpeter equation: idea

As mentioned above, it makes sense to pull apart the lattice bubble $\chi^{0}$
into two parts, the impurity bubble $X^{0}$ and a remainder
$\tilde{\chi}^{0}$,

$\chi^{0}=X^{0}+\tilde{\chi}^{0}.$ (4)

We are looking for an analoguous division of $\chi$,

$\chi=X+\overline{\chi}$ (5)

in terms of the impurity susceptibility $X$ and a correction term
$\overline{\chi}$ which should be determined. Insertion of the Bethe-Salpeter
equations [Eqs. (2) and (3)] gives

$\displaystyle\overline{\chi}$
$\displaystyle=\chi-X=\frac{\chi^{0}}{1-\chi^{0}\Gamma}-\frac{X^{0}}{1-X^{0}\Gamma}$
(6)
$\displaystyle=\frac{X^{0}+\tilde{\chi}^{0}}{1-(X^{0}+\tilde{\chi}^{0})\Gamma}-\frac{X^{0}}{1-X^{0}\Gamma}$
(7)
$\displaystyle=\frac{1}{1-X^{0}\Gamma}\frac{\tilde{\chi}^{0}}{1-\tilde{\chi}^{0}\frac{\Gamma}{1-\Gamma
X^{0}}}\frac{1}{1-X^{0}\Gamma}\,.$ (8)

The result can be further simplified by writing it in terms of the impurity
reducible vertex $F$, defined as

$X=X^{0}+X^{0}FX^{0}\,,\quad F=\Gamma+\Gamma X^{0}F$ (9)

which gives the relations

$\frac{\Gamma}{1-\Gamma X^{0}}=F\,,\quad\frac{1}{1-X^{0}\Gamma}=1+X^{0}F\equiv
L.$ (10)

Here, we have introduced $L$ to compactify the notation further.

Hence, the DMFT+DBSE approach gives the lattice susceptibility $\chi$ in terms
of the impurity susceptibility $X$ and the dual correction term

$\chi=X+L\frac{\tilde{\chi}^{0}}{1-\tilde{\chi}^{0}F}L\equiv
X+L\tilde{\chi}L\,.$ (11)

Note that this formulation of the lattice susceptibility only uses the
reducible vertex $F$ and is therefore free from the issues related to
divergences that occur in $\Gamma$ [33, 34, 35, 36].

Both formulations of the Bethe-Salpeter equation, Eqs. (2) and (11), are
matrix equations, and the matrix multiplication involves a formally infinite
sum over fermionic frequencies. In practice, this sum has to be cut off at
some $N_{\nu}$. In Sec. III, we will show that the Eq. (11) converges more
rapidly with respect to $N_{\nu}$, making it superior for practical
implementations. Before studying the convergence, however, we will first
introduce how the calculation is done in practice, with all relevant orbital
and frequency labels and prefactors.

###  II.5 Impurity correlators

The impurity susceptibility $X_{\bar{a}b\bar{c}d}(\omega)$ is the Fourier
transform to Matsubara frequency of Eq. (1) for the impurity model, i.e.,
without a spatial label $\mathbf{R}$.

Similarly, for the fermion-boson vertex $L$, we start with the impurity
correlation function $g^{(3)}_{\bar{a}b\bar{c}d}(\omega,\nu)=\left\langle
c^{\dagger}_{\nu,a}c^{\phantom{\dagger}}_{\nu+\omega,b}(c^{\dagger}_{c}c^{\phantom{\dagger}}_{d})_{\omega}\right\rangle$,
where $\nu$ is a fermionic and $\omega$ a bosonic Matsubara frequency. $L$ is
the amputated connected component of $g^{(3)}$, i.e.,

$-\frac{1}{T}g_{e\bar{a}}(\nu)g_{b\bar{f}}(\nu+\omega)L_{\bar{e}f\bar{c}d}(\omega,\nu)\\\
=g^{(3)}_{\bar{a}b\bar{c}d}(\omega,\nu)-\frac{1}{T}g_{b\bar{a}}(\nu)g_{d\bar{c}}(\tau=0^{+})\delta_{\omega,0}.$
(12)

Here, $g$ is the single-particle Green’s function of the impurity model. To
solve Eq. (12) for $L$, it is convenient to express it in terms of the
impurity bubble,

$X^{0}_{\bar{a}b\bar{c}d}(\omega,\nu,\nu^{\prime})=-\frac{1}{T}\delta_{\nu\nu^{\prime}}g_{d\bar{a}}(\nu)g_{b\bar{c}}(\nu+\omega),$
(13)

so that Eq. (12) becomes

$X^{0}_{\bar{a}b\bar{f}e}(\omega,\nu,\nu)L_{\bar{e}f\bar{c}d}(\omega,\nu)\\\
=g^{(3)}_{\bar{a}b\bar{c}d}(\omega,\nu)-\frac{1}{T}g_{b\bar{a}}(\nu)g_{d\bar{c}}(\tau=0^{+})\delta_{\omega,0}.$
(14)

This is a tensorial equation in orbital space while being diagonal in both the
fermionic and the bosonic Matsubara frequency. $L$ is obtained by
multiplication with the inverse of $X^{0}$, TPRF’s implementation of these
operations is used and this is the reason for the factor $-1/T$ in our
definition of $L$.

The fermion-boson vertex has some useful properties, which have been discussed
in detail for single-orbital DMFT [37]. First, the impurity susceptibility $X$
can be obtained from $g^{(3)}$ by performing a trace over $\nu$. Secondly, $L$
describes the linear response of the self-energy to an external field,

$-\frac{1}{T}\sum_{cd}B_{\bar{d}c}L_{\bar{a}b\bar{c}d}(\omega=0,\nu)=B_{\bar{a}b}+\frac{\partial\Sigma_{b\bar{a}}(\nu)}{\partial\lambda}.$
(15)

Thirdly, for a non-interacting impurity model,
$L_{\bar{a}b\bar{c}d}(\omega,\nu)=-T\delta_{\bar{a}d}\delta_{\bar{c}b}$, is
constant, since
$-Tg^{(3)}_{\bar{a}b\bar{c}d}(\omega,\nu)=g_{d\bar{a}}(\nu)g_{b\bar{c}}(\nu+\omega)-g_{b\bar{a}}(\nu)g_{d\bar{c}}(\tau=0^{+})\delta_{\omega,0}$.
This expression for $L$ is consistent with the linear response of Eq. (15). It
can also be shown more generally that $L$ asymptotically goes to a constant in
the limit of large $\nu$, keeping $\omega$ fixed. Finally, there is the
symmetry relation

$\displaystyle
L_{\bar{a}b\bar{c}d}(\omega,\nu-\omega/2)=L_{\bar{b}a\bar{d}c}(-\omega,\nu+\omega/2).$
(16)

This formula is obtained by reversing the arrow of all fermionic lines (time-
reversal), and re-arranging the diagram to bring it back into the customary
form.

From our definition in Eq. (12), $L$ has units of energy since $T$ has units
of energy, $g(\nu)$ has units of inverse energy and $g(\tau)$ is
dimensionless. For the linear response, it is convenient to define the
operators $\hat{A}$ and $\hat{B}$ in a dimensionless way (this includes the
important examples of $\hat{N}$ and $\hat{S}^{i}$), so $\lambda$ has units of
energy and the right-hand side of Eq. (15) is dimensionless.

Another ingredient for the dual Bethe-Salpeter equation is the reducible
impurity vertex, $F_{\bar{a}b\bar{c}d}$, which is the connected, amputated
two-particle correlation function of the impurity model. This object is used
in many DMFT-based theories [7], and a variety of notations is used in the
literature. We use the same convention as TPRF for the particle-hole channel,
i.e.,

$g^{(4)}_{\bar{a}b\bar{c}d}(\omega,\nu,\nu^{\prime})=\\\
g_{e\bar{a}}(\nu)g_{b\bar{f}}(\nu\\!+\\!\omega)F_{\bar{e}f\bar{g}h}(\omega,\nu,\nu^{\prime})g_{g\bar{c}}(\nu^{\prime}\\!+\\!\omega)g_{d\bar{h}}(\nu^{\prime})\\\
+\beta\delta_{0,\omega}g_{b\bar{a}}(\nu)g_{d\bar{c}}(\nu^{\prime})-\beta\delta_{\nu,\nu^{\prime}}g_{d\bar{a}}(\nu)g_{b\bar{c}}(\nu+\omega)\,.$
(17)

###  II.6 Dual Bethe-Salpeter equation in multi-orbital systems

The lattice Green’s function $G$ in DMFT is defined by

$\delta_{a\bar{c}}=\left[i\nu\delta_{a\bar{b}}-t_{a\bar{b}}(\mathbf{k})-\Sigma_{a\bar{b}}(\nu)\right]G_{b\bar{c}}(\nu,\mathbf{k}),$
(18)

where $\Sigma$ is the self-energy of the auxiliary impurity model. We define
the non-local single-particle Green’s function $\tilde{G}$ as

$\tilde{G}_{a\bar{b}}(\nu,\mathbf{R})=\begin{cases}G_{a\bar{b}}(\nu,\mathbf{R})&\text{
if }\mathbf{R}\neq 0\\\ 0&\text{ if }\mathbf{R}=0\end{cases},$ (19)

where $G_{a\bar{b}}(\nu,\mathbf{R})$ is the Fourier transform of
$G_{a\bar{b}}(\nu,\mathbf{k})$ to real space. Corresponding to the non-local
Green’s function $\tilde{G}$, there is a “bubble” $\tilde{\chi}^{0}$,

$\tilde{\chi}^{0}_{\bar{a}b\bar{c}d}(\omega,\nu,\nu^{\prime},\mathbf{R})=-\frac{1}{T}\delta_{\nu\nu^{\prime}}\tilde{G}_{d\bar{a}}(\nu,\mathbf{R})\tilde{G}_{b\bar{c}}(\nu+\omega,-\mathbf{R}).$
(20)

This object is Fourier transformed back to momentum space, giving
$\tilde{\chi}^{0}_{\bar{a}b\bar{c}d}(\omega,\nu,\nu^{\prime},\mathbf{q})$.

From $\tilde{\chi}^{0}$ and $F$, the dual or non-local Bethe-Salpeter equation

$\displaystyle\tilde{\chi}_{\bar{a}b\bar{c}d}(\omega,\nu,\nu^{\prime},\mathbf{q})=\tilde{\chi}^{0}_{\bar{a}b\bar{c}d}(\omega,\nu,\nu^{\prime},\mathbf{q})+$
(21)
$\displaystyle\,\,\,\,\tilde{\chi}^{0}_{\bar{a}b\bar{f}e}(\omega,\nu,\nu_{1},\mathbf{q})F_{\bar{e}f\bar{g}h}(\omega,\nu_{1},\nu_{2})\tilde{\chi}_{\bar{h}g\bar{c}d}(\omega,\nu_{2},\nu^{\prime},\mathbf{q}),$

can now be solved for $\tilde{\chi}$ by inversion.

Finally, to get the physical susceptibility, we have to combine the impurity
susceptibility with the non-local correction, adding fermion-boson vertices at
the end,

$\displaystyle\chi_{\bar{a}b\bar{c}d}(\omega,\mathbf{q})=X_{\bar{a}b\bar{c}d}(\omega)$
$\displaystyle\,\,\,\,+L_{\bar{f}e\bar{a}b}(-\omega,\nu+\omega)\tilde{\chi}_{\bar{e}f\bar{g}h}(\omega,\nu,\nu^{\prime},\mathbf{q})L_{\bar{h}g\bar{c}d}(\omega,\nu^{\prime}+\omega)$
(22)

Using the symmetry of $L$, this can be rewritten as

$\displaystyle\chi_{\bar{a}b\bar{c}d}(\omega,\mathbf{q})=X_{\bar{a}b\bar{c}d}(\omega)$
$\displaystyle\,\,\,\,+L_{\bar{e}f\bar{b}a}(\omega,\nu)\tilde{\chi}_{\bar{e}f\bar{g}h}(\omega,\nu,\nu^{\prime},\mathbf{q})L_{\bar{h}f\bar{c}d}(\omega,\nu^{\prime}+\omega)$
(23)

Equation (23) has the advantage that it only uses a single bosonic frequency
$\omega$, and also performs better in terms of frequency windows when $\omega$
is large, since $\nu+\omega$ is replaced by $\nu$. On the other hand, Equation
(22) is obviously invariant under an orbital basis transformation, while this
is not obvious to see in Equation (22).

The computationally non-trivial step in this procedure is the inversion of the
dual Bethe-Salpeter equation in Eq. (21). It corresponds to a geometric series
with kernel $\tilde{\chi}^{0}F$, which is divergent if one of the eigenvalues
of the kernel is equal to $+1$, signalling a phase transition. Thus, analysis
of this kernel can give valuable insight into the leading electronic
fluctuations [36].

$L_{\bar{a}b\bar{c}d}(\omega,\nu)=$$b$$\nu+\omega$$>$$a$$\nu$$<$$>$$<$$d$$c$$\omega$

$>$$<$$\tilde{\chi}$$f$$g$$\nu\phantom{{}^{\prime}}+\omega$$\nu^{\prime}+\omega$$e$$h$$\nu$$\nu^{\prime}$$>$$<$$>$$<$$a$$b$$d$$c$$\omega,\mathbf{q}$$\omega,\mathbf{q}$

Figure 2: Diagrammatic representation of elements of the dual Bethe-Salpeter
equation.

###  II.7 Comparison of dual and usual Bethe-Salpeter equation

The usual Bethe-Salpeter equation in DMFT is based on the Green’s function $G$
instead of $\tilde{G}$, and the vertex $\Gamma$ instead of $F$. The
susceptibility is then obtained directly as the output of the Bethe-Salpeter
equation, i.e., without needing Eq. (22). Diagrammatically, it corresponds to
the repeated scattering of electron-hole pairs propagating with the DMFT
Green’s function and interacting with the irreducible vertex $\Gamma$.

Comparing the two formulations of the Bethe-Salpeter equation, the dual boson
based approach can be interpreted as a resummation of the original Bethe-
Salpeter equation where terms in the geometric series are grouped together if
subsequent vertices $\Gamma$ lie at the same lattice site. This occurs in the
following way: If all vertices $\Gamma$ involve only a single lattice site,
the diagram is a part of $X$. Otherwise, all vertices before the first change
of site are grouped into the $L$ on the left and all vertices after the last
change of site are grouped into the $L$ on the right. In between, all vertices
$\Gamma$ that occur between two changes of site are grouped into a single $F$.
To avoid double-counting of diagrams, the dual bubble $\tilde{\chi}^{0}$ is
used to enforce a change of sites between every vertex.

This construction provides for a one-to-one correspondence of the original
Bethe-Salpeter series and the new, dual Bethe-Salpeter with all elements
expressed in terms of $\Gamma$ and bubbles.

###  II.8 Implementation details

After converging an ordinary DMFT self-consistency loop using TRIQS and the
CTHYB solver [38], we use w2dynamics [39] to measure the two-particle
correlation functions with one, two and three frequencies. Due to different
conventions the output of the w2dynamics solver is then transformed to the
TPRF two-particle convention, see Appendix A for details. After post-
processing the two-particle correlation functions to get $\chi^{\text{imp}}$,
$L$ and $F$, the dual Bethe-Salpeter equation is evaluated and
$\chi(\omega,\mathbf{q})$ is extracted as in Sec. II.6. This evaluation is the
main contribution of the present manuscript.

The DMFT+DBSE calculation is compute bound by the measurements of the impurity
correlation functions (using QMC). The dual Bethe-Salpeter equation solver is
memory bound, since we use fast Fourier transforms to evaluate the lattice
bubble suceptibility $\chi^{0}$ as direct products in real space and imaginary
time ($\mathbf{R}$, $\tau$). However, the DBSE equation is diagonal in
$\omega$ and $\mathbf{q}$ and can be evaluated one $(\omega,\mathbf{q})$ at a
time. All computationally non-trivial parts of the DBSE calculation are
performed using the two-particle response function toolbox (TPRF) with
routines implemented in C++, using hybrid OpenMP+MPI parallelization.

##  III Convergence with respect to the number of fermionic frequencies

An important benefit of the dual version of the susceptibility calculation is
the improved convergence with respect to the number of fermionic frequencies
$N_{\nu}$. Whereas the traditional Bethe-Salpeter equation scales as
$N_{\nu}^{-1}$, the new implementation generically scales as $N_{\nu}^{-3}$,
as shown below.

###  III.1 The non-interacting dimer

Consider the Hubbard dimer in the limit $U=0$, i.e., two atoms with a hopping
amplitude $t$. In that case, the Hamiltonian has two eigenvalues
$\epsilon(k)=\pm t$ and all calculations can be done straightforwardly. Define
$\operatorname{sgn}(k)=\operatorname{sgn}(\epsilon_{k})$, so
$\epsilon_{k}=t\operatorname{sgn}(k)$.

$\displaystyle G(i\nu,k)$
$\displaystyle=\frac{1}{i\nu-t\operatorname{sgn}(k)},$ (24)
$\displaystyle\frac{1}{N_{k}}\sum_{k}G(i\nu,k)$
$\displaystyle=\frac{1}{2}\left(\frac{1}{i\nu-t}+\frac{1}{i\nu+t}\right),$
(25) $\displaystyle=-\frac{i\nu}{\nu^{2}+t^{2}}$ (26)
$\displaystyle\tilde{G}(i\nu,k)$
$\displaystyle=\operatorname{sgn}(k)\frac{t}{\nu^{2}+t^{2}}.$ (27)

The non-local bubble susceptibility at $q=0$ and $i\omega=0$ is

$\displaystyle\tilde{\chi}^{0}(\omega=0,\nu,\nu^{\prime})$
$\displaystyle=-\beta\delta_{\nu\nu^{\prime}}\frac{1}{N_{k}}\sum_{k}\tilde{G}(\nu,k)\tilde{G}(\nu,k)$
(28)
$\displaystyle=-\beta\delta_{\nu\nu^{\prime}}\frac{t^{2}}{(\nu^{2}+t^{2})^{2}}$
(29)
$\displaystyle=-\beta\delta_{\nu\nu^{\prime}}\frac{t^{2}}{\nu^{4}}\frac{1}{(1+t^{2}/\nu^{2})^{2}}.$
(30)

From the final expression, it is clear that this frequency-diagonal matrix
element decays as $(i\nu)^{-4}$. Given a fermionic frequency box of size
$N_{\nu}$, the estimated total contribution from frequencies outside of the
box scales as $\int_{N_{\nu}}^{\infty}dx/x^{4}\propto N_{\nu}^{-3}$. Thus, we
find cubic convergence.

For the dimer, it is easy to verify that the other value of $q$ has the same
magnitude and only differs in the sign. Thus, it also has the same scaling. A
finite bosonic frequency $\omega$ does not change the asymptotic scaling
(since it is irrelevant when $\nu\gg\omega$), although it does change when the
asymptotic scaling sets in.

The expectation is that this scaling also holds up in interacting systems. In
that case, one has to attach $L$ to both ends of the diagram, but $L$ is
asymptotically constant in $\nu$, and one needs to insert vertex corrections
to $\tilde{\chi}^{0}$ via a geometric equation, but these are not expected to
change the scaling of $\tilde{\chi}^{0}$.

Note that the first term in the dual formula, the impurity susceptibility, is
measured exactly in the impurity solver without introducing a fermionic
frequency box.

For completeness, the conventional evaluation of the susceptibility in the
dimer, at $\omega=0$ and $q=0$, involves

$\displaystyle\chi^{0}(\nu,\nu^{\prime})$
$\displaystyle=-\beta\delta_{\nu\nu^{\prime}}\frac{1}{N_{k}}\sum_{k}G(\nu,k)G(\nu,k)$
(31)
$\displaystyle=-\beta\delta_{\nu\nu^{\prime}}\frac{1}{2}\left[\frac{1}{(i\nu-t)^{2}}+\frac{1}{(i\nu+t)^{2}}\right]$
(32)
$\displaystyle=-\beta\delta_{\nu\nu^{\prime}}\frac{-\nu^{2}+t^{2}}{(\nu^{2}+t^{2})^{2}}.$
(33)

This matrix element decays as $(i\nu)^{-2}$, so the total cutoff error scales
as $N_{\nu}^{-1}$, i.e., linearly.

The difference between $\chi^{0}$ and $\tilde{\chi}^{0}$ is given by the
impurity susceptibility, which satisfies

$\displaystyle X$
$\displaystyle=-\beta\sum_{\nu}\left(\frac{1}{N_{k}}\sum_{k}G(i\nu,k)\right)^{2}$
(34) $\displaystyle=-\beta\sum_{\nu}\frac{-\nu^{2}}{(\nu^{2}+t^{2})^{2}}.$
(35)

However, this fermionic frequency sum is not performed explicitly in the fast
evaluation of the susceptibility. Instead, $X$ is measured directly in the
Monte Carlo.

###  III.2 The non-interacting Hubbard model

Figure 3: Magnetic susceptibility at $\mathbf{q}=0$, $\omega=0$ in a two-
orbital square lattice tight-binding model ($t=1$, $U=0$, $\mu\approx-1.84$)
as a function of the number of frequencies $N_{\nu}$. The dual Bethe-Salpeter
equation (DBSE) converges as $1/N_{\nu}^{3}$, whereas the usual Bethe-Salpeter
equation (BSE) converges only as $1/N_{\nu}$, as is visible in the log-log
plot on the right.

The scaling in the non-interacting Hubbard model is the same as in the
corresponding non-interacting dimer, with the difference that the momentum
sums cannot be done analytically. The essential point is that
$G(i\nu,\mathbf{k})\sim 1/i\nu+t_{k}/(i\nu)^{2}+\ldots$, so the leading
contribution is local ($\mathbf{k}$ independent) and gets removed from
$\tilde{G}$, which decays as $1/(i\nu)^{2}$ for large $\nu$. Since the bubble
contains two Green’s functions, $\tilde{\chi}$ decays as $\nu^{-4}$ compared
to the $\nu^{-2}$ scaling of $\chi$. Finally, the integral over all remaining
frequencies leads to $N_{\nu}^{-3}$ versus $N_{\nu}^{-1}$ scaling. A numerical
illustration of the scaling is given in Figure 3.

As a reference, Figure 3 uses the exact $\chi(\mathbf{q}=0,\omega=0)=dm/dh$
obtained using finite differences from the magnetization $m(h)$, given by the
Fermi-Dirac distribution for non-interacting fermions. To achieve consistency,
the same momentum space discretization should be used in the BSE and the
linear response. Note that the non-interacting model with orbital-diagonal
hopping is actually just a sum over independent single-orbital models, so the
orbital physics is trivial, but it is a useful test of the general
implementation.

###  III.3 Atomic limit

In the atomic limit, $t_{\mathbf{k}}=0$, we have $\chi(\omega,q)=X(\omega)$
and the new formulation gives the exact result regardless of the number of
fermionic frequencies. This is consistent with the observation that the
leading term in the non-local bubble is asymptotically of the form
$t_{\mathbf{k}}/(i\nu)^{2}$, which vanishes in the atomic limit. On the other
hand, the conventional formulation of the Bethe-Salpeter equation does depend
on the finite box size in the usual linear way, since the local Bethe-Salpeter
equation of the Hubbard atom is non-trivial [40].

##  IV Benchmarks

To benchmark the dual BSE formulation, we apply it to two strongly correlated
models, namely the single-band Hubbard model on the square lattice and the
three-band effective low-energy model for strontium ruthenate (Sr2RuO4).

###  IV.1 The single-band Hubbard model

Figure 4: Static magnetic susceptibility $\chi(\mathbf{0},0)$ ($\mathbf{q}=0$,
$\omega=0$) in the square lattice Hubbard model as a function of the number of
frequencies $N_{\nu}$. The dual BSE result (DBSE) (blue) is compared with the
linear response result from self-consistent DMFT in an applied field
extrapolated to zero field (DMFT+Field) (orange and inset).

The single band Hubbard model on the square lattice is a canonical system for
the study of strong correlations. Here we benchmark the static spin
susceptibility, computed using dual BSE, against the extrapolated zero field
response from self consistent DMFT calculations in an applied magnetic field.
The calculations are performed in the strongly correlated regime with nearest
neighbour hopping $t=1$ and local Hubbard interaction $U=10$ in the high
temperature paramagnetic state at inverse temperature $\beta=1$.

Figure 4 shows the convergence of the DBSE static susceptibility (blue
markers) using fermionic frequency windows in the range $4\leq N_{\nu}\leq
30$, which gives four digits of accuracy and a quantitative agreement with the
applied field result (orange circle). In fact, the discrepancy between the
extrapolated DBSE result (blue circle) and the applied field (orange circle)
can be attributed to MC erros and the extrapolation to zero field in the
linear response calculation (as shown in the inset).

The expected cubic scaling with the size of the frequency window $N_{\nu}$ is
clearly visible in the DBSE result, which allows us to reach a higher accuracy
than the conventional BSE with linear convergence.

###  IV.2 Three-band Hubbard model for strontium ruthenate

Figure 5: Upper panel: Static spin susceptibility
$\chi_{S_{z}S_{z}}(\mathbf{q})$ of strontium ruthenate Sr2RuO4 at $T=464$ K
along the pseudo-tetragonal high-symmetry path ($\Gamma$-$X$-$M$-$\Gamma$).
Results for $N_{\nu}=4$, 8, and 16 fermionic frequencies are shown for the
DMFT+BSE method (dashed lines) and our dual BSE method (DMFT+DBSE) (solid
lines). Lower panel: Convergence at the $X$ point as a function of $1/N_{\nu}$
for $4\leq N_{\nu}\leq 30$. Extrapolatins to the infinite frequency limit
$N_{\nu}\rightarrow\infty$ are also shown (lines), using the $N_{\nu}^{-3}$
scaling (solid line) for DMFT+DBSE results (circles). Due to the small
$N_{\nu}$ values the DMFT+BSE results (squares) do not scale as $N_{\nu}^{-1}$
(dashed line). However, a $N_{\nu}^{-2}$ corrected extrapolation (dotted line)
agrees quantitatively with DMFT+DBSE.

While computing the lattice susceptibility within DMFT for the single band
Hubbard model to higher accuracy is a feat, the main motivation for
implementing the DBSE is to push the boundaries of applicability of the method
to material realistic models in order to connect to experimental
spectroscopies. This generally requires DMFT treatment of multi-orbital
models. For example, the correlated electrons in transition metal compounds
can often be described by effective low-energy Hubbard models with two to five
orbitals per unit cell.

One such strongly correlated material is strontium ruthenate (Sr2RuO4) which
has a substantial momentum-dependence in its spin susceptibility due to the
presence of incommensurate spin fluctuations, as seen in INS experiments [41,
42] and reproduced by DMFT+BSE calculations [43, 44, 31]. The combination of
orbital and momentum structure is responsible for the observed phenomena, with
both strong local fluctuations and sharp features in momentum space. This
makes Sr2RuO4 suitable as a benchmark for our dual Bethe-Salpeter equation
implementation. Here, we focus on the frequency convergence of the
implementation, the physics of the spin susceptibility in this material is
discussed at length in Ref. 31.

Following Ref. 31, we construct a low-energy model of Sr2RuO4 from the three
bands crossing the Fermi level with Ru-4d t2g symmetry. The corresponding
Wannier model was fit using Wannier90 [45] on a $10^{3}$ grid to the density
functional theory band structure calculated with GPAW [46, 47], using the PBE
exchange correlation functional [48], with a $20^{3}$ k-point grid, $600$ eV
plane wave cutoff, and $25$ meV Methfessel-Paxton smearing [49], using the
experimental crystal structure of Sr2RuO4 at $100\,$K [50]. The interaction in
the low energy three band t2g model was modelled using the rotationally
invariant Kanamori interaction, with Hubbard $U=2.3$ eV and Hund’s $J=0.4$ eV
and the self-consistent dynamical mean field calculations were performed at
$T=464$ K using TRIQS [26] and TRIQS/cthyb [38] on a $16^{3}$ k-point grid.
The two-particle quantities were measured using W2Dynamics [51] worm sampling
[52] using $4\cdot 10^{9}$ samples and a total compute budget of 110’000 core
hours on the EuroHPC Joint Undertaking system LUMI. The static spin-
suseptibility $\chi_{S_{z}S_{z}}(\mathbf{q},\omega=0)$ was finally computed on
the level of DMFT+BSE and DMFT+DBSE using the Two-Particle Response Function
toolbox TPRF [30] using a $48^{3}$ k-point grid.

Figure 5 shows the resulting susceptibility along a high-symmetry path through
the Brillouin Zone. Results are shown for both the traditional implementation
of the Bethe-Salpeter equation (BSE) and our dual implementation (DBSE). In
both cases, relatively small frequency boxes are used ($N_{\nu}=4,8,16$). For
the conventional BSE, the slow $N_{\nu}^{-1}$ convergence is observed, and
extrapolation is needed to obtain qualitatively correct results. On the other
hand, the DBSE yields qualitative results already for the smallest box size of
$N_{\nu}=4$, which is a truly remarkable improvement. Enlarging the frequency
box, there is a small enhancement of the momentum-dependence, but no overall
qualitative change. This is in stark contrast with the conventional BSE, where
enlarging the fermionic frequency box increases the susceptibility across the
entire Brillouin Zone. The origin of this difference is clear: local spin
fluctuations are the major contribution to $\chi(\mathbf{q})$ and the dual
Bethe-Salpeter equation includes all of these in one go via the impurity
susceptibility, independent of the size of the fermionic frequency box.

For the computational efficiency of the approach, the cost of obtaining the
two-particle correlation functions of the impurity model is also important.
For this particular case, the worm sampling [52] used here (as implemented in
w2dynamics) is orders of magnitude more efficient than the partition function
sampling used in Ref. [31]. Together, the improved convergence of the dual
Bethe-Salpeter equation and the efficiency of the worm sampling make the
calculation of DMFT lattice susceptibilities in multi-orbital models a routine
task instead of a gargantuan undertaking.

##  V Conclusion

We have presented an efficient implementation of the DMFT susceptibility for
general multi-orbital systems based on the dual boson formalism. This
algorithm reduces the scaling with the inverse number of frequencies from
linear to cubic, leading to a dramatic speed-up of the calculation. A
comparison with linear response shows that convergence with respect to the
frequency box is much easier to achieve, leading to several digits of accuracy
in the resulting static susceptibility in our example, even before
extrapolation with respect to the size of the frequency box. Note that the
improved frequency scaling we obtain is due to the dual Bethe-Salpeter
reformulation in combination with a simple frequency cut-off. The improved
scaling originates from applying the frequency cut-off on a different set of
two-particle objects, with faster decaying high frequency asumptotics, than
the terms in the standard DMFT Bethe-Salpeter equation.

Regarding the general treatment of the Bethe-Salpeter equation, going beyond
the drastic hard frequency cut-off in Matsubara space, there are several
interesting approaches that incorporate the asymptotic frequency dependence in
different ways [14, 15, 16, 17, 18] or that use a compact representation of
the generic two frequency behaviour [53, 54]. Generally, the dual Bethe-
Salpeter equation would benefit from these refined approaches, further
improving the convergence.

The implementation of the dual Bethe-Salpeter equation is based on the TRIQS
toolbox (in particular: TPRF) and can in principle be used together with any
impurity solver that provides the required two-particle correlation functions.
An interface for the w2dynamics solver is provided with the implementation. We
have found the worm sampling of two-particle correlation functions in
w2dynamics to be very efficient for the susceptibilities studied here.
Compared to the usual Bethe-Salpeter equation, two additional impurity
correlation functions with zero and one fermionic frequencies, respectively.
These are generally similar in computational cost to the vertex with two
fermionic frequencies that is needed in both versions of the Bethe-Salpeter
equation. However, Monte Carlo sampling issues have been observed for
situations where the observables in the susceptibility do not commute with the
Hamiltonian [55].

This efficient implementation of the Bethe-Salpeter equation can be useful for
the study of magnetic [56, 31] and other susceptibilities [57] in systems with
multiple correlated orbitals in the unit cell. In particular, the better
scaling with respect to the number of fermionic Matsubara frequencies makes it
possible to reach lower temperatures.

###### Acknowledgements.

EvL acknowledges support from Gyllenstiernska Krapperupstiftelsen, Crafoord
Foundation and from the Swedish Research Council (Vetenskapsrådet, VR) under
grant 2022-03090. EvL also acknowledges support by eSSENCE, a strategic
research area for e-Science, grant number eSSENCE@LU 9:1. HURS acknowledges
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (Grant agreement No.
854843-FASTCORR). The computations were enabled by resources provided by the
National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the
Swedish National Infrastructure for Computing (SNIC) through the projects LU
2023/2-37, LU 2022/2-32, SNIC 2022/23-304, SNIC 2022/21-15, SNIC 2022/13-9,
SNIC 2022/6-113, and SNIC 2022/1-18, at Lunarc, PDC, NSC and CSC partially
funded by the Swedish Research Council through grant agreements no. 2022-06725
and no. 2018-05973.

##  Appendix A Conventions for two-particle quantities

As discussed in the main text, the two-particle Green’s functions of the
impurity problem are a necessary ingredient for the dual Bethe-Salpeter
equation. They are obtained from the impurity solver, in this case TRIQS/cthyb
or w2dynamics. There are small differences in the definition of these objects
between the two codes, which are essential for a correct implementation. This
appendix provides their relation.

###  A.1 Two-particle Green’s function in the particle-hole channel

The four-point response function in TRIQS/cthyb for the particle-hole channel
(PH) is

$g^{(4)ph}_{abcd}(\omega,\nu,\nu^{\prime})=\frac{1}{\beta}\iiiint_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}d\tau_{4}\\\
\exp\left[i\omega(\tau_{2}-\tau_{3})+i\nu(\tau_{2}-\tau_{1})+i\nu^{\prime}(\tau_{4}-\tau_{3})\right]\\\
\langle\mathcal{T}c^{\dagger}_{a}(\tau_{1})c_{b}(\tau_{2})c^{\dagger}_{c}(\tau_{3})c_{d}(\tau_{4})\rangle$
(36)

In the w2dynamics documentation, the four-point response function for the
particle-hole channel (PH) is defined as

$\tilde{g}^{(4)ph}_{abcd}(\omega,\nu,\nu^{\prime})=\iiiint_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}d\tau_{4}\\\
\exp\left[i\omega(\tau_{1}-\tau_{4})+i\nu(\tau_{1}-\tau_{2})+i\nu^{\prime}(\tau_{3}-\tau_{4})\right]\\\
\langle\mathcal{T}c_{a}(\tau_{1})c^{\dagger}_{b}(\tau_{2})c_{c}(\tau_{3})c^{\dagger}_{d}(\tau_{4})\rangle\\\
=\iiiint_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}d\tau_{4}\\\
\exp\left[i\omega(\tau_{2}-\tau_{3})+i\nu(\tau_{2}-\tau_{1})+i\nu^{\prime}(\tau_{4}-\tau_{3})\right]\\\
\langle\mathcal{T}c^{\dagger}_{b}(\tau_{1})c_{a}(\tau_{2})c^{\dagger}_{d}(\tau_{3})c_{c}(\tau_{4})\rangle=\beta
g^{(4)ph}_{badc}(\omega,\nu,\nu^{\prime})$ (37)

Note that the bosonic frequency convention in W2Dynamics is controlled by
setting the parameter `WormPHConvention=0` to get the
$i\omega(\tau_{1}-\tau_{4})$ factor (setting `WormPHConvention=1` gives
$i\omega(\tau_{2}-\tau_{3})$ and is not directly compatible with
TRIQS/cthyb+TPRF).

In numerical comparisons, we have found that the actual output data of both
solver is related as

$g^{(4)ph}_{abcd}(\omega,\nu,\nu^{\prime})=\beta\tilde{g}^{(4)ph}_{badc}(\omega,\nu,\nu^{\prime}),$
(38)

which differs from Eq. (37) by a factor $\beta^{2}$. This is the relation used
in our implementation. It should be noted that factors of $\beta$ are not
always written explicitly in parts of the literature, which can make it harder
to translate equations from one notation to another.

###  A.2 Two-particle Green’s function in the particle-hole channel with
$\tau_{3}=\tau_{4}$

The two-particle Green’s function with $\tau_{3}=\tau_{4}$ following the
TRIQS/cthyb PH format from Eq. (36) is

$g^{(3)ph}_{abcd}(\omega,\nu)=\iiint_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}\\\
\times\exp[i\omega(\tau_{2}-\tau_{3})+i\nu(\tau_{2}-\tau_{1})]\\\
\times\langle\mathcal{T}c^{\dagger}_{a}(\tau_{1})c_{b}(\tau_{2})c^{\dagger}_{c}(\tau_{3})c_{d}(\tau_{3})\rangle.$
(39)

The sampled quantity as defined from W2Dynamics is

$\tilde{g}^{(3)ph}_{abcd}(\omega,\nu)=\iiint_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}\\\
\times\exp[i\omega(\tau_{2}-\tau_{3})+i\nu(\tau_{1}-\tau_{2})]\\\
\times\langle\mathcal{T}c_{a}(\tau_{1})c^{\dagger}_{b}(\tau_{2})c_{c}(\tau_{3})c^{\dagger}_{d}(\tau_{3})\rangle\\\
=\iiint_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}\exp[i\omega(\tau_{1}-\tau_{3})+i\nu(\tau_{2}-\tau_{1})]\\\
\times\langle\mathcal{T}c^{\dagger}_{b}(\tau_{1})c_{a}(\tau_{2})c^{\dagger}_{d}(\tau_{3})c_{c}(\tau_{3})\rangle\\\
=\iiint_{0}^{\beta}d\tau_{1}d\tau_{2}d\tau_{3}\exp[i\omega(\tau_{2}-\tau_{3})+i(\nu-\omega)(\tau_{2}-\tau_{1})]\\\
\times\langle\mathcal{T}c^{\dagger}_{b}(\tau_{1})c_{a}(\tau_{2})c^{\dagger}_{d}(\tau_{3})c_{c}(\tau_{3})\rangle=g^{(3)ph}_{badc}(\omega,\nu-\omega).$
(40)

Note that the bosonic frequency convention produces a shift in the fermionic
frequency. The implementation is not sensitive to the `WormPHConvention`
parameter. Hence both orbital label pair-permutations and a bosonic frequency
shift are required to transform from W2Dynamics to TRIQS format

$g^{(3)ph}_{abcd}(\omega,\nu)=\tilde{g}^{(3)ph}_{badc}(\omega,\nu+\omega).$
(41)

The fermionic frequency shift, which is somewhat inconvenient, can be avoided
by using time-reversal symmetry, as in Eq. (16),

$\displaystyle
g^{(3)ph}_{badc}(-\omega,\nu)=g^{(3)ph}_{abcd}(\omega,\nu-\omega)=\tilde{g}^{(3)ph}_{badc}(\omega,\nu).$
(42)

## References

  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).
  * Tamai  _et al._ [2019] A. Tamai, M. Zingl, E. Rozbicki, E. Cappelli, S. Riccò, A. de la Torre, S. McKeown Walker, F. Y. Bruno, P. D. C. King, W. Meevasana, M. Shi, M. Radović, N. C. Plumb, A. S. Gibbs, A. P. Mackenzie, C. Berthod, H. U. R. Strand, M. Kim, A. Georges, and F. Baumberger, High-resolution photoemission on Sr2RuO4 reveals correlation-enhanced effective spin-orbit coupling and dominantly local self-energies, Phys. Rev. X 9, 021048 (2019).
  * Coldea and Watson [2018] A. I. Coldea and M. D. Watson, The key ingredients of the electronic structure of FeSe, Annual Review of Condensed Matter Physics 9, 125 (2018), https://doi.org/10.1146/annurev-conmatphys-033117-054137 .
  * Phillips  _et al._ [2022] P. W. Phillips, N. E. Hussey, and P. Abbamonte, Stranger than metals, Science 377, eabh4273 (2022).
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, Correlated lattice fermions in $d=\infty$ dimensions, Phys. Rev. Lett. 62, 324 (1989).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Rohringer  _et al._ [2018] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held, Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory, Rev. Mod. Phys. 90, 025003 (2018).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard model: A computational perspective, Annual Review of Condensed Matter Physics 13, 275 (2022), https://doi.org/10.1146/annurev-conmatphys-090921-033948 .
  * Ribic  _et al._ [2017] T. Ribic, P. Gunacker, S. Iskakov, M. Wallerberger, G. Rohringer, A. N. Rubtsov, E. Gull, and K. Held, Role of three-particle vertex within dual fermion calculations, Phys. Rev. B 96, 235127 (2017).
  * Kappl  _et al._ [2022] P. Kappl, F. Krien, C. Watzenböck, and K. Held, Non-linear responses and three-particle correlators in correlated electron systems exemplified by the Anderson impurity model (2022).
  * Enderle [2014] M. Enderle, Neutrons and magnetism, Collection SFN 13 (2014).
  * Platzman and Wolff [1973] P. M. Platzman and P. A. Wolff,  _Waves and interactions in solid state plasmas_ , Vol. 13 (Academic Press New York, 1973).
  * Ament  _et al._ [2011] L. J. P. Ament, M. van Veenendaal, T. P. Devereaux, J. P. Hill, and J. van den Brink, Resonant inelastic x-ray scattering studies of elementary excitations, Rev. Mod. Phys. 83, 705 (2011).
  * Kuneš [2011] J. Kuneš, Efficient treatment of two-particle vertices in dynamical mean-field theory, Phys. Rev. B 83, 085102 (2011).
  * Boehnke  _et al._ [2011] L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann, and O. Parcollet, Orthogonal polynomial representation of imaginary-time Green’s functions, Phys. Rev. B 84, 075145 (2011).
  * Kaufmann  _et al._ [2017] J. Kaufmann, P. Gunacker, and K. Held, Continuous-time quantum Monte Carlo calculation of multiorbital vertex asymptotics, Phys. Rev. B 96, 035114 (2017).
  * Tagliavini  _et al._ [2018] A. Tagliavini, S. Hummel, N. Wentzell, S. Andergassen, A. Toschi, and G. Rohringer, Efficient Bethe-Salpeter equation treatment in dynamical mean-field theory, Phys. Rev. B 97, 235140 (2018).
  * Wentzell  _et al._ [2020] N. Wentzell, G. Li, A. Tagliavini, C. Taranto, G. Rohringer, K. Held, A. Toschi, and S. Andergassen, High-frequency asymptotics of the vertex function: Diagrammatic parametrization and algorithmic implementation, Phys. Rev. B 102, 085106 (2020).
  * Rubtsov  _et al._ [2012] A. Rubtsov, M. Katsnelson, and A. Lichtenstein, Dual boson approach to collective excitations in correlated fermionic systems, Annals of Physics 327, 1320 (2012).
  * Hafermann  _et al._ [2014] H. Hafermann, E. G. C. P. van Loon, M. I. Katsnelson, A. I. Lichtenstein, and O. Parcollet, Collective charge excitations of strongly correlated electrons, vertex corrections, and gauge invariance, Phys. Rev. B 90, 235105 (2014).
  * van Loon  _et al._ [2014a] E. G. C. P. van Loon, A. I. Lichtenstein, M. I. Katsnelson, O. Parcollet, and H. Hafermann, Beyond extended dynamical mean-field theory: Dual boson approach to the two-dimensional extended Hubbard model, Phys. Rev. B 90, 235135 (2014a).
  * van Loon  _et al._ [2015] E. G. C. P. van Loon, H. Hafermann, A. I. Lichtenstein, and M. I. Katsnelson, Thermodynamic consistency of the charge response in dynamical mean-field based approaches, Phys. Rev. B 92, 085106 (2015).
  * van Loon  _et al._ [2016] E. G. C. P. van Loon, F. Krien, H. Hafermann, E. A. Stepanov, A. I. Lichtenstein, and M. I. Katsnelson, Double occupancy in dynamical mean-field theory and the dual boson approach, Phys. Rev. B 93, 155162 (2016).
  * van Loon  _et al._ [2014b] E. G. C. P. van Loon, H. Hafermann, A. I. Lichtenstein, A. N. Rubtsov, and M. I. Katsnelson, Plasmons in strongly correlated systems: Spectral weight transfer and renormalized dispersion, Phys. Rev. Lett. 113, 246407 (2014b).
  * van Loon  _et al._ [2018a] E. G. C. P. van Loon, M. Rösner, G. Schönhoff, M. I. Katsnelson, and T. O. Wehling, Competing Coulomb and electron–phonon interactions in NbS2, npj Quantum Materials 3, 10.1038/s41535-018-0105-4 (2018a).
  * Parcollet  _et al._ [2015] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio, and P. Seth, Triqs: A toolbox for research on interacting quantum systems, Computer Physics Communications 196, 398 (2015).
  * Pavarini  _et al._ [2014] E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein,  _DMFT at 25: Infinite dimensions: Lecture notes of the autumn school on correlated electrons 2014_ , Vol. 4 (Forschungszentrum Jülich, 2014).
  * Boehnke [2015] L. V. Boehnke,  _Susceptibilities in materials with multiple strongly correlated orbitals_ , Ph.D. thesis, Staats-und Universitätsbibliothek Hamburg Carl von Ossietzky (2015).
  * Krien [2018] F. Krien,  _Conserving dynamical mean-field approaches to strongly correlated systems_ , Ph.D. thesis, Staats-und Universitätsbibliothek Hamburg Carl von Ossietzky (2018).
  * Strand [2019] H. U. R. Strand, Two-particle response function tool-box (TPRF) for TRIQS, github.com/TRIQS/tprf (2019).
  * Strand  _et al._ [2019] H. U. R. Strand, M. Zingl, N. Wentzell, O. Parcollet, and A. Georges, Magnetic response of ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$: Quasi-local spin fluctuations due to Hund’s coupling, Phys. Rev. B 100, 125120 (2019).
  * Krien [2019] F. Krien, Efficient evaluation of the polarization function in dynamical mean-field theory, Phys. Rev. B 99, 235106 (2019).
  * Schäfer  _et al._ [2013] T. Schäfer, G. Rohringer, O. Gunnarsson, S. Ciuchi, G. Sangiovanni, and A. Toschi, Divergent precursors of the Mott-Hubbard transition at the two-particle level, Phys. Rev. Lett. 110, 246405 (2013).
  * Chalupa  _et al._ [2018] P. Chalupa, P. Gunacker, T. Schäfer, K. Held, and A. Toschi, Divergences of the irreducible vertex functions in correlated metallic systems: Insights from the Anderson impurity model, Phys. Rev. B 97, 245136 (2018).
  * Reitner  _et al._ [2020] M. Reitner, P. Chalupa, L. Del Re, D. Springer, S. Ciuchi, G. Sangiovanni, and A. Toschi, Attractive effect of a strong electronic repulsion: The physics of vertex divergences, Phys. Rev. Lett. 125, 196403 (2020).
  * van Loon  _et al._ [2020] E. G. C. P. van Loon, F. Krien, and A. A. Katanin, Bethe-Salpeter equation at the critical end point of the Mott transition, Phys. Rev. Lett. 125, 136402 (2020).
  * van Loon  _et al._ [2018b] E. G. C. P. van Loon, F. Krien, H. Hafermann, A. I. Lichtenstein, and M. I. Katsnelson, Fermion-boson vertex within dynamical mean-field theory, Phys. Rev. B 98, 205148 (2018b).
  * Seth  _et al._ [2016] P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, TRIQS/CTHYB: A continuous-time quantum Monte Carlo hybridisation expansion solver for quantum impurity problems, Computer Physics Communications 200, 274 (2016).
  * Wallerberger  _et al._ [2019a] M. Wallerberger, A. Hausoel, P. Gunacker, A. Kowalski, N. Parragh, F. Goth, K. Held, and G. Sangiovanni, w2dynamics: Local one- and two-particle quantities from dynamical mean field theory, Computer Physics Communications 235, 388 (2019a).
  * Thunström  _et al._ [2018] P. Thunström, O. Gunnarsson, S. Ciuchi, and G. Rohringer, Analytical investigation of singularities in two-particle irreducible vertex functions of the Hubbard atom, Phys. Rev. B 98, 235107 (2018).
  * Sidis  _et al._ [1999] Y. Sidis, M. Braden, P. Bourges, B. Hennion, S. NishiZaki, Y. Maeno, and Y. Mori, Evidence for incommensurate spin fluctuations in ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$, Phys. Rev. Lett. 83, 3320 (1999).
  * Steffens  _et al._ [2019] P. Steffens, Y. Sidis, J. Kulda, Z. Q. Mao, Y. Maeno, I. I. Mazin, and M. Braden, Spin fluctuations in ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$ from polarized neutron scattering: Implications for superconductivity, Phys. Rev. Lett. 122, 047004 (2019).
  * Boehnke  _et al._ [2018] L. Boehnke, P. Werner, and F. Lechermann, Multi-orbital nature of the spin fluctuations in Sr2RuO4, EPL 122, 57001 (2018).
  * Acharya  _et al._ [2019] S. Acharya, D. Pashov, C. Weber, H. Park, L. Sponza, and M. V. Schilfgaarde, Evening out the spin and charge parity to increase $T_{c}$ in Sr2RuO4, Communications Physics 2, 163 (2019).
  * Mostofi  _et al._ [2008] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, wannier90: A tool for obtaining maximally-localised Wannier functions, Comput. Phys. Commun. 178, 685 (2008).
  * Mortensen  _et al._ [2005] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Real-space grid implementation of the projector augmented wave method, Phys. Rev. B 71, 035109 (2005).
  * Enkovaara  _et al._ [2010] J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljungberg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen, T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G. A. Tritsaris, M. Vanin, M. Walter, B. Hammer, H. Häkkinen, G. K. H. Madsen, R. M. Nieminen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz, K. S. Thygesen, and K. W. Jacobsen, Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method, J. Phys.: Condens. Matter 22, 253202 (2010).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Methfessel and Paxton [1989] M. Methfessel and A. T. Paxton, High-precision sampling for Brillouin-zone integration in metals, Phys. Rev. B 40, 3616 (1989).
  * Vogt and Buttrey [1995] T. Vogt and D. J. Buttrey, Low-temperature structural behavior of ${\mathrm{Sr}}_{2}$${\mathrm{RuO}}_{4}$, Phys. Rev. B 52, R9843 (1995).
  * Wallerberger  _et al._ [2019b] M. Wallerberger, A. Hausoel, P. Gunacker, A. Kowalski, N. Parragh, F. Goth, K. Held, and G. Sangiovanni, w2dynamics: Local one- and two-particle quantities from dynamical mean field theory, Comput. Phys. Commun. 235, 388 (2019b).
  * Gunacker  _et al._ [2015] P. Gunacker, M. Wallerberger, E. Gull, A. Hausoel, G. Sangiovanni, and K. Held, Continuous-time quantum Monte Carlo using worm sampling, Phys. Rev. B 92, 155102 (2015).
  * Shinaoka  _et al._ [2018] H. Shinaoka, J. Otsuki, K. Haule, M. Wallerberger, E. Gull, K. Yoshimi, and M. Ohzeki, Overcomplete compact representation of two-particle green’s functions, Phys. Rev. B 97, 205111 (2018).
  * Wallerberger  _et al._ [2021] M. Wallerberger, H. Shinaoka, and A. Kauch, Solving the bethe-salpeter equation with exponential convergence, Phys. Rev. Res. 3, 033168 (2021).
  * van Loon and Strand [2023] E. G. C. P. van Loon and H. U. R. Strand, Larmor precession in strongly correlated itinerant electron systems (2023), arXiv:2303.03468 [cond-mat.str-el] .
  * Hausoel  _et al._ [2017] A. Hausoel, M. Karolak, E. Şaş$\iota$oğlu, A. Lichtenstein, K. Held, A. Katanin, A. Toschi, and G. Sangiovanni, Local magnetic moments in iron and nickel at ambient and earth’s core conditions, Nature communications 8, 16062 (2017).
  * Geffroy  _et al._ [2019] D. Geffroy, J. Kaufmann, A. Hariki, P. Gunacker, A. Hausoel, and J. Kuneš, Collective modes in excitonic magnets: Dynamical mean-field study, Phys. Rev. Lett. 122, 127601 (2019).