Efficient evaluation of the polarization function in dynamical mean-field theory

TL;DR

This paper presents a fast and accurate method for evaluating the polarization function in dynamical mean-field theory, reducing computational costs and improving the description of critical behavior in correlated electron systems.

Contribution

The authors introduce an exact decomposition of the DMFT polarization function, enabling efficient calculation of the lattice response with minimized cutoff errors and improved critical behavior analysis.

Findings

The new method reduces frequency summation complexity.

It accurately captures critical behavior near antiferromagnetic transitions.

Application to Hubbard models demonstrates improved results.

Abstract

The dynamical susceptibility of strongly correlated electronic systems can be calculated within the framework of the dynamical mean-field theory (DMFT). The required measurement of the four-point vertex of the auxiliary impurity model is however costly and restricted to a finite grid of Matsubara frequencies, leading to a cutoff error. It is shown that the propagation of this error to the lattice response function can be minimized by virtue of an exact decomposition of the DMFT polarization function into local and nonlocal parts. The former is measured directly by the impurity solver, while the latter is given in terms of a ladder equation for the Hedin vertex that features an unprecedentedly fast decay of frequency summations compared to previous calculation schemes, such as the one of the dual boson approach. At strong coupling the local approximation of the TRILEX approach is viable, but vertex corrections to the polarization should be dropped on equal footing to recover the correct prefactor of the effective exchange. In finite dimensions the DMFT susceptibility exhibits spurious mean-field criticality, therefore, a two-particle self-consistent and frequency-dependent correction term is introduced, similar to the Moriya-$\lambda$ correction of the dynamical vertex approximation. Applications to the two- and three-dimensional Hubbard models on the square and cubic lattices show that the expected critical behavior near an antiferromagnetic instability is recovered.

Equations

$$X_{\mathbf{q}}(\omega )=\chi (\omega )+{\tilde{X}}_{\mathbf{q}}(\omega ),$$

$$H=$$

$$G_k=[ı\nu -{\epsilon }_{\mathbf{k}}+\mu -\Sigma (\nu ){]}^{-1},$$

$$S_{\text{AIM}}=$$

$$\sum _{\mathbf{k}}{G}_k=g_{\nu },$$

$$g_{\nu {\nu }^{\prime }\omega }^{(4),\alpha }=$$

$$g_{\nu \omega }^{(3),\alpha }=$$

$$f_{\nu {\nu }^{\prime }\omega }^{\alpha }=$$

$${\bar{\lambda }}_{\nu \omega }^{\alpha }=$$

$${\lambda }_{\nu \omega }^{\alpha }=$$

$$f_{\nu {\nu }^{\prime }\omega }^{\alpha }=$$

$$X_q^{\alpha }=\frac{2{\Pi }_q^{\alpha }}{1-U^{\alpha }{\Pi }_q^{\alpha }},$$

$${\Pi }_q^{\alpha }=$$

$${\tilde{X}}_{\nu }^0(q)=\sum _{\mathbf{k}}{\tilde{G}}_{k+q}{\tilde{G}}_k.$$

$${\Lambda }_{\nu q}^{i,\alpha }=$$

$${\tilde{X}}_q^{\alpha }=2\sum _{\nu }{\Lambda }_{\nu q}^{\alpha }{\tilde{X}}_{\nu }^0(q){\bar{\lambda }}_{\nu \omega }^{\alpha },$$

$${\tilde{X}}_q/2(\text{or}{\tilde{\Pi }}_q)=$$

$$+$$

$$f_{\nu {\nu }^{\prime }\omega }^{\alpha }=$$

$$\underset{|\nu |\to \infty }{\lim }{f}_{\nu {\nu }^{\prime }\omega }^{\alpha }=$$

$$f_{\nu {\nu }^{\prime }\omega }^{i,\alpha }=$$

$${\Pi }_q^{(2),\alpha }=$$

$${\Pi }_q^{(1),\alpha }=$$

$$=$$

$$X^{\text{sp}}(\mathbf{q},\omega =0)\widetilde{=}-\frac{2}{2T-I_{\mathbf{q}}},$$

$$I_{\mathbf{q}}^{(2)}=$$

$$I_{\mathbf{q}}^{(1)}=$$

$$X_q^{\alpha }\to {\mathcal{X}}_q^{\alpha }=$$

$$\sum _{\mathbf{q}}{\mathcal{X}}_q^{\alpha }={\chi }_{\omega }^{\alpha },$$

$$\sum _q{\mathcal{X}}_q^{\text{ch}}=$$

$$\sum _q{\mathcal{X}}_q^{\text{sp}}=$$

$$c^{(i),\alpha }(\bar{n})=\sum _{n,n^{\prime }=-\bar{n}}^{\bar{n}-1}{g}^2({\nu }_n)f^{(i),\alpha }({\nu }_n,{\nu }_{n^{\prime }}^{\prime },\omega =0)g^2({\nu }_{n^{\prime }}^{\prime }),$$

$$G_{k{k}^{\prime }q}^{(4),\alpha }$$

Full text

I Introduction

The dynamical mean-field theory (DMFT) is a powerful non-perturbative approach to strong local correlations in the Hubbard model Georges et al. (1996). Although in widespread use, many aspects of the DMFT are still under investigation, which is fueled to large extent by persistent algorithmic advances in the solution of its auxiliary Anderson impurity model Gull et al. (2011). These improvements allow insights into the two-particle level of the DMFT approximation Rohringer et al. (2012), which is also the elemental precursor for its diagrammatic extensions Rohringer et al. (2018).

A basic application for DMFT at the two-particle level is the calculation of the dynamical susceptibility, which allows to study, for example, phase transitions Georges et al. (1996), the electron energy loss spectrum van Loon et al. (2014a), nuclear relaxation rate Boehnke and Lechermann (2012), and Goldstone excitations Geffroy et al. (2019). The DMFT susceptibility is furthermore an integral part of the ladder dynamical vertex approximation Toschi et al. (2007); Galler et al. (2017). Calculation of this correlation function however requires knowledge of the impurity vertex function, which is often evaluated by means of improved estimators for continuous-time quantum Monte-Carlo (CTQMC) solvers Hafermann et al. (2012); Gunacker et al. (2016). The further development of improved estimators is highly desirable, as they allow to efficiently calculate the DMFT susceptibility in multi-orbital settings, see for example Refs. Boehnke and Lechermann (2012); Boehnke et al. (2018); Geffroy et al. (2019). Recently, progress has been reported in the measurement of the vertex function within the exact diagonalization method Tanaka (2019).

The role of the improved estimators in CTQMC methods is to minimize the statistical noise of the Monte Carlo measurement, which for fixed run-time greatly increases with the number of dynamic degrees of freedom (Matsubara frequencies) of the measured quantity. A further numerical error is introduced because the measurement of the impurity vertex function is restricted to a finite grid of Matsubara frequencies. In order to obtain a gauge invariant lattice response function in DMFT it is necessary to account for an infinite number, that is, a ladder of vertex corrections Hafermann et al. (2014). For each vertex correction the value of the impurity vertex at all frequencies enters the calculation, and therefore due to the finite Matsubara grid a cutoff error arises. Consequently, the numerical error of the DMFT response function may not only be minimized by an improved Monte Carlo measurement but also by reduction of the cutoff error. A straightforward way to do this is to account for the asymptotics of the vertex function Kuneš (2011); Wentzell et al. (2016); Kaufmann et al. (2017); Tagliavini et al. (2018).

A further option for improvement, the subject of this work, is to use the numerically exact impurity solver to sum local diagrams exactly. For concreteness, within the dual boson approach and in a calculation scheme by Pruschke et al. the DMFT susceptibility $X$ is written as the sum of local and nonlocal parts Pruschke et al. (1996); Rubtsov et al. (2012); Hafermann et al. (2014),

[TABLE]

where $\mathbf{q}$ is the lattice momentum and $\omega$ the (bosonic) Matsubara frequency. The local part, the impurity susceptibility $\chi$, depends only on one frequency and is calculated directly by the impurity solver, which in effect sums all local two-particle diagrams that taken together yield $\chi$. Moreover, even though DMFT neglects nonlocal correlations, the lattice susceptibility takes local vertex corrections at different lattice sites into account, which give rise to the nonlocal term $\tilde{X}$ (see also Fig. 1 of Ref. van Loon et al. (2016)). The dual boson formula (1) is numerically efficient because the statistical and cutoff errors attached to the impurity vertex only affect the nonlocal term, not $\chi$. Even when the asymptotic behavior of the impurity vertex is neglected it allows the analytical continuation of the susceptibility to the real axis Hafermann et al. (2014); van Loon et al. (2014a). It is however desirable to preserve numerical resources, hence further improvements are welcome.

In this work it will be shown that the concept of breaking down the DMFT susceptibility into simpler diagrammatic pieces can be taken to a further level by separating exactly the diagrams from the vertex function that are irreducible with respect to the bare Hubbard interaction $U$. This leads to a decomposition of the polarization function $\Pi$, which is $U$-irreducible, into local and nonlocal parts Stepanov et al. (2016a), $\Pi_{\mathbf{q}}(\omega)=\pi(\omega)+\tilde{\Pi}_{\mathbf{q}}(\omega)$, analogous to equation (1). The nonlocal part $\tilde{\Pi}$ is obtained via an efficient ladder equation for the Hedin three-leg vertex Hedin (1965). The lattice polarization $\Pi_{\mathbf{q}}(\omega)$ in turn encapsulates all non-trivial information about the two-particle spectrum.

The Hedin vertex also plays a central role in the TRILEX approach Ayral and Parcollet (2016). In this method nonlocal vertex corrections to the Hedin vertex are neglected, and therefore the calculation of the four-point vertex function of the impurity model is not necessary. However, this approximation can be introduced in different ways, for example, within the dual boson formalism it accounts for more vertex corrections than in TRILEX Stepanov et al. (2016a). It is shown in this work that for large interaction these additional vertex corrections decide about the prefactor of the effective exchange coupling. Both TRILEX and dual boson account for a nonlocal self-energy, however, this work focuses on approximations to the polarization function.

Lastly, a further aspect is considered in the application of the efficient formula for the polarization: In finite dimensions the DMFT susceptibility violates the Pauli principle and suffers from a spurious mean-field instability in two dimensions. It has been shown previously that the Mermin-Wagner theorem is satisfied in the renormalized ladder dual fermion approach Otsuki et al. (2014) or after introduction of the Moriya-$\lambda$ correction to the DMFT susceptibility Katanin et al. (2009). In three dimensions both approaches renormalize the criticality of the underlying dynamical mean-field starting point Rohringer et al. (2011); Hirschmeier et al. (2015). Furthermore, in the ladder dynamical vertex approximation the satisfaction of local charge and spin sum rules by the Moriya-$\lambda$ correction is crucial to ensure the proper asymptotic behavior of the electronic self-energy Katanin et al. (2009); Rohringer and Toschi (2016). Similar to the Moriya-$\lambda$ and two-particle self-consistent approach (TPSC) Y.M. Vilk and A.-M.S. Tremblay (1997), in this work the mean-field artifacts of the DMFT susceptibility are removed by virtue of a frequency-dependent correction term that is fixed by a two-particle self-consistent constraint. It is shown that this approach satisfies the Mermin-Wagner theorem and predicts the same criticality of the half-filled three-dimensional Hubbard model as the Moriya-$\lambda$.

The paper is organized as follows: The Hubbard Hamiltonian, the DMFT approximation, and the Anderson impurity model are briefly recollected in Sec. II. The reducible and irreducible vertices of the impurity model are defined in Sec. III. The efficient formula for the DMFT polarization is presented in Sec. IV and compared to the dual boson formula. A two-particle self-consistent modification of the DMFT susceptibility and a TRILEX-like approximation are introduced in Sec. V and applied in Sec. VI. The conclusions follow in Sec. VII. A self-contained derivation of the ladder equation for the Hedin vertex is provided in the Appendices A-D.

II Hubbard Hamiltonian and DMFT approximation

The Hamiltonian of the paramagnetic two- or three-dimensional Hubbard model is given as,

[TABLE]

where ${t}_{ij}$ is the nearest neighbor hopping between lattice sites $i,j$, its absolute value ${t}=1$ is the unit of energy. $c,c^{\dagger}$ are the construction operators, $\sigma=\uparrow,\downarrow$ the spin index. $U$ is the Hubbard repulsion between the densities $n_{\sigma}=c^{\dagger}_{\sigma}c_{\sigma}$.

In the DMFT approximation the self-energy $\Sigma$ of Green’s function is local,

[TABLE]

where $k=(\mathbf{k},\nu)$ comprises lattice momentum and fermionic Matsubara frequency, $\varepsilon_{\mathbf{k}}$ is the dispersion, $\mu$ is the chemical potential. $\Sigma(\nu)$ is the self-energy of an auxiliary Anderson impurity model (AIM) that is solved numerically exactly. The action of the AIM reads,

[TABLE]

Here $\Delta_{\nu}$ denotes the hybridization function, $\omega$ is a bosonic Matsubara frequency, summations $\sum_{\nu},\sum_{\omega}$ imply multiplication with the temperature $T$. $c^{*},c$ are Grassmann numbers. In DMFT the hybridization function is fixed self-consistently according to the constraint,

[TABLE]

where $g$ is the numerically exact local Green’s function of the AIM. Note that summation over $\mathbf{k}$ implies division by the number of lattice sites $N$.

III Impurity vertices

The calculation of the dynamical susceptibility requires knowledge of higher correlation functions of the impurity. Directly measured by the solver are the susceptibility, $\chi^{\alpha}_{\omega}=-\langle{\rho^{\alpha}_{-\omega}\rho^{\alpha}_{\omega}}\rangle+\beta\langle n\rangle\langle n\rangle\delta_{\omega}\delta_{\alpha,\text{ch}}$, the four-point,

[TABLE]

and the three-point function,

[TABLE]

where $s^{\alpha}$ are the Pauli matrices ($\alpha=\text{ch},\text{sp}$), $\rho^{\text{ch}}=n_{\uparrow}+n_{\downarrow}$ and $\rho^{\text{sp}}=n_{\uparrow}-n_{\downarrow}$ are the charge and spin densities.

III.1 Reducible vertices

One defines the four- and three-point vertices $f$ and $\bar{\lambda}$,

[TABLE]

Although numerically unfavorable $\bar{\lambda}$ can in principle also be obtained by attaching legs to $f$ from the right and adding $1$, $\bar{\lambda}^{\alpha}_{\nu\omega}=1+\sum_{\nu^{\prime}}f^{\alpha}_{\nu\nu^{\prime}\omega}g_{\nu^{\prime}}g_{\nu^{\prime}+\omega}$, therefore $\bar{\lambda}$ is a right-sided three-leg vertex van Loon et al. (2018a), the left-sided one $\lambda$ is obtained by attaching the legs from the left or via the symmetry relation, $\bar{\lambda}^{\alpha}_{\nu\omega}=\lambda^{\alpha}_{\nu+\omega,-\omega}$.

III.2 $U^{\alpha}$-irreducible vertices

In order to make the later calculation of the DMFT lattice correlation functions efficient the impurity vertices are decomposed following Hertz and Edwards Hertz and Edwards (1973):

The diagram $a)$ in Fig. 1 shows that when the three-leg vertex $\lambda^{\alpha}$ is expanded diagrammatically one may encounter, in going from left to right, an insertion of the bare interaction $U^{\alpha}$, where $U^{\text{ch}}=+U$ or $U^{\text{sp}}=-U$. The Hubbard interaction is just a constant, and the incoming impurity Green’s function lines on the left of $U^{\alpha}$ can thus be contracted, the same is case for the out-going lines.

On the left of $U^{\alpha}$ there hence arises a contribution to the $U^{\alpha}$-irreducible polarization $\pi^{\alpha}$ of the impurity [related to the susceptibility via $\chi^{\alpha}_{\omega}=2\pi^{\alpha}_{\omega}/(1-U^{\alpha}\pi^{\alpha})$], whereas on the right of $U^{\alpha}$ begins once again an expansion of the three-leg vertex. As shown algebraically in Appendix A, one thus separates diagrams from $\lambda$ that are once or manifold $U^{\alpha}$-reducible,

[TABLE]

where $\lambda^{i}$ is the $U^{\alpha}$-irreducible three-leg vertex – the Hedin vertex – of the impurity.

Let us perform this procedure also for the four-point vertex $f$, as depicted in Fig. 1 b). $f$ obviously contains one part $f^{i}$ that is irreducible, whereas in the remaining terms one finds at least one insertion $U^{\alpha}$. At this point the incoming lines may be closed and a right-sided Hedin vertex $\bar{\lambda}^{i}$ arises on the left of $U^{\alpha}$. In fact, also on the right of $U^{\alpha}$ the lines may be closed, which means that a true four-point contribution does not arise in the $U^{\alpha}$-reducible diagrams. For this reason the whole of the reducible diagrams may be split into the three- and two-point objects $\bar{\lambda}^{i},\lambda^{i}$ and $\pi$, respectively,

[TABLE]

where $w^{\alpha}_{\omega}={U^{\alpha}}/({1-U^{\alpha}\pi^{\alpha}_{q}})$ is the screened interaction of the impurity [cf. Fig. 1 c)].

The equations (8) and (9) are valuable because they separate RPA-like diagrams from the vertices $\lambda$ and $f$, which are absorbed into the geometric series in Fig. 1 c), the screened interaction $w$ 111 The polarization $\pi$ can indeed be interpreted as the ‘self-energy’ of the screened interaction $w$, analogous to the Dyson equation $g=g^{0}/(1-g^{0}\Sigma)$, and the bare interaction $U$ assumes the role of the bare Green’s function $g^{0}$. On the other hand, $U$ also corresponds to the two-particle self-energy of the RPA approximation Mahan (2000); Vilk and Tremblay (1996), one may therefore refer to the diagrams in Fig. 1 c) as ‘RPA-like’. . Similar relations are also valid for the Hubbard model (2), see Ref. Held et al. (2011) and Appendix A. The characteristic triangle-wiggle-triangle diagram in Fig. 1 b) is typically large when the corresponding susceptibility $\chi^{\alpha}$ is large, since then $U^{\alpha}\pi^{\alpha}\approx 1$.

One should note that in the reducible contribution $\bar{\lambda}^{i}_{\nu\omega}w_{\omega}\lambda^{i}_{\nu^{\prime}\omega}$ in Eq. (9) the dependence on $\nu$ and $\nu^{\prime}$ is separated. Therefore, this term is comprised in the lowest order of a singular value decomposition of $f$ Otsuki et al. (2019).

IV Efficient formula

The goal is to calculate the dynamical susceptibility in the DMFT approximation, see also definition (53),

[TABLE]

where $q=(\mathbf{q},\omega)$ and $\Pi^{\alpha}_{q}$ is the lattice polarization. The form of equation (10) resembles the RPA susceptibility, however, the polarization $\Pi$ is similar to the Lindhard function only in the weak coupling limit, while for intermediate and large coupling $\Pi$ is strongly renormalized by the frequency dependence of the DMFT self-energy $\Sigma(\nu)$ and vertex corrections Hafermann et al. (2014). The latter can be taken into account in the following efficient way:

It is shown in Appendix C that in DMFT the polarization can be decomposed into local and nonlocal parts Stepanov et al. (2016a),

[TABLE]

The nonlocal corrections are denoted as $\tilde{\Pi}$, analogous to the dual boson formula for the susceptibility (1) and $\tilde{X}^{0}_{\nu}(q)$ is a nonlocal bubble,

[TABLE]

Here, $\tilde{G}_{k}=G_{k}-g_{\nu}$ is the nonlocal DMFT Green’s function, which decays with the frequency as ${1}/{\nu^{2}}$, and $\Lambda^{i,\alpha}$ is the left-sided lattice Hedin vertex. Equation (11) is depicted diagrammatically in Fig. 2 b).

We now come to the main result, which is a nonlocal ladder equation for the Hedin vertex in the DMFT approximation Katanin et al. (2009). It is shown in Appendix D that,

[TABLE]

which is depicted in Fig. 2 a). Note that $f^{i}$ is the $U^{\alpha}$-irreducible four-leg vertex of the impurity model, it is the only true four-point object needed in the calculation.

IV.1 Comparison to dual boson formula

It will now be shown that the formula (11) for the polarization is numerically more efficient than the dual boson formula (1). To this end, let us recall that in the latter case the nonlocal corrections are given as [see Appendix C, Refs. van Loon et al. (2014a); Hafermann et al. (2014), and Fig. 2 c)],

[TABLE]

similar to $\tilde{\Pi}$ in Eq. (11), except that the $U$-reducible three-leg vertices $\Lambda,\lambda$ are in place of the Hedin vertices $\Lambda^{i},\lambda^{i}$ (and the factor $2$). Furthermore, the vertex $\Lambda$ of the lattice is given by the same ladder equation (13) [see also Fig. 2 a)], albeit the label ’$i$’ needs to be omitted, and there is hence a complete formal analogy in the calculation of $\Pi$ and $X$.

Let us compare the first four-point vertex contribution to the nonlocal correction terms $\tilde{X}$ and $\tilde{\Pi}$, by expanding the ladder equations for the three-leg vertices $\Lambda$ and $\Lambda^{i}$, respectively, see also Eq. (13),

[TABLE]

where the flavor label $\alpha$ was omitted for readability.

Typically the calculation of the impurity three-leg vertices $\lambda^{(i)}$ is more efficient than that of the four-leg vertices $f^{(i)}$, in the latter case one likes to minimize the domain of measurement for $\nu,\nu^{\prime},\omega$. The question is therefore how the cutoff error in the four-point corrections that arise in the second line of Eq. (15) affects the calculation. It is useful to analyze the convergence of the term that is written out in the second line of Eq. (15), let us consider first the limit $|\nu|\rightarrow\infty$ while $\nu^{\prime}$ and $\omega$ are kept constant:

According to Eq. (8) the decay of the vertices $\lambda_{\nu\omega}$ and $\lambda^{i}_{\nu\omega}$ with the frequency $\nu$ is the same except for a prefactor $[1-U^{\alpha}\pi^{\alpha}_{\omega}]^{-1}$, therefore, the difference in the three-leg vertices does not lead to a different convergence of the $\nu$-summations in $\tilde{X}$ and $\tilde{\Pi}$. Also in both cases the nonlocal bubble $\tilde{X}^{0}_{\nu}(q)$ defined in Eq. (12) decays as $1/\nu^{4}$. However, the vertices $f$ and $f^{i}$ behave differently, which follows from an observation in Ref. Rohringer and Toschi (2016): In the limit $|\nu|\rightarrow\infty$ all diagrams contributing to $f$ that depend on $\nu$ have decayed, and hence asymptotically this vertex is given by the diagrams that do not depend on $\nu$ at all. According to the argument in the reference these diagrams are all $U$-reducible, one can write for fixed $\nu^{\prime}$,

[TABLE]

Factoring out $U^{\alpha}$ one identifies the reducible three-leg vertex $\lambda$ [see below Eq. (7)], therefore,

[TABLE]

In the last step Eq. (7) and $w^{\alpha}=U^{\alpha}/(1-U^{\alpha}\pi^{\alpha}_{\omega})$ were used. Let us now compare to the exact relation between the vertices $f$ and $f^{i}$ in Eq. (9). The asymptotic limit of $f$ in Eq. (17) is given exactly by the asymptotic limit of the $U$-reducible diagrams $\bar{\lambda}^{i}_{\nu\omega}\,w_{\omega}\,\lambda^{i}_{\nu^{\prime}\omega}$ (note that $\bar{\lambda}^{i}_{\nu\omega}\rightarrow 1$ for $|\nu|\rightarrow\infty$). This is not surprising in view of the observation of Ref. Rohringer and Toschi (2016) that only $U$-reducible diagrams can be independent of $\nu$. As a result, the irreducible vertex $f^{i}_{\nu\nu^{\prime}\omega}$ decays to zero for $|\nu|\rightarrow\infty$ and fixed $\nu^{\prime}$,

[TABLE]

For this reason the four-point corrections in Eq. (15) decay by at least one power of $\nu$ faster for $\tilde{\Pi}$ than for $\tilde{X}$, which is the central observation of this work.

A comprehensive discussion of the asymptotics of $f$ can be found in Ref. Wentzell et al. (2016), where it is also shown that in the double limit $|\nu|,|\nu^{\prime}|\rightarrow\infty$ one needs to consider separately the two cases $\nu-\nu^{\prime}=\text{const}$ and $\omega+\nu+\nu^{\prime}=\text{const}$, that is, the elements of $f$ near the main and secondary diagonal. However, as regards the scope of this work these cases can be ignored, because then the nonlocal bubbles in Eq. (15) decay as $1/\nu^{4}$ and $1/(\nu^{\prime})^{4}$, respectively, leading to a still faster decay than when only one frequency is large. In summary, in the dual boson formula each four-point correction comes with a factor $\tilde{X}^{0}_{\nu}(q)f_{\nu\nu^{\prime}\omega}$, which decays like the nonlocal bubble as $1/\nu^{4}$ due to the constant background of $f$, whereas in the efficient calculation scheme the corrections enter as $\tilde{X}^{0}_{\nu}(q)f^{i}_{\nu\nu^{\prime}\omega}$, which decays at least as $1/\nu^{5}$ by virtue of the combined decay of nonlocal bubble and vertex $f^{i}$.

V TRILEX-like approximation and two-particle self-consistency

This section considers an optimal truncation of the vertex corrections to the Hedin vertex and a two-particle self-consistent constraint on the DMFT susceptibility.

V.1 TRILEX-like approximation

Despite all optimizations it may be unfeasible to take four-point corrections to the Hedin vertex into account, for example, in multi-orbital settings, cf. Appendix E. In this case one may consider to neglect vertex corrections in Eq. (13), $\Lambda^{i}\approx\lambda^{i}$, which is the philosophy of the TRILEX approach. However, there are two ways to introduce this approximation: Firstly, in the efficient formula (11) the local approximation to the Hedin vertex leads to,

[TABLE]

which corresponds to replacing the full triangle in Fig. 2 b) with a bare triangle and one is left with two bare triangles. Secondly, a more direct way to apply the approximation is to insert it into the relation $\Pi^{\alpha}_{q}=\sum_{k}\Lambda^{i,\alpha}_{\nu q}G_{k}G_{k+q}$, which is equivalent to equation (11) when vertex corrections are kept (cf. Appendix C). Nevertheless, $\Lambda^{i}\approx\lambda^{i}$ leads to a different approximation,

[TABLE]

In the second line the nonlocal bubble (12) was introduced using the relation $\sum_{\mathbf{k}}G_{k}G_{k+q}=\tilde{X}^{0}_{\nu}(q)+g_{\nu}g_{\nu+\omega}$ and the exact impurity polarization was identified, $\pi^{\alpha}_{\omega}=\sum_{\nu}\lambda^{i,\alpha}_{\nu\omega}g_{\nu}g_{\nu+\omega}$. Equation (20) corresponds to the way the local approximation is introduced in the TRILEX approach Ayral and Parcollet (2016), it has only one bare triangle.

The obvious question is whether the first or the second option is a more viable way to truncate the vertex corrections. This question can be decided by considering the strong coupling limit, which shows that only $\Pi^{(2)}$ in equation (19) correctly describes the effective exchange: It is shown in Appendix F that for very large coupling $U\gg t,T$ the static DMFT spin susceptibility of the half-filled Hubbard model takes the form, see also Ref. Otsuki et al. (2019),

[TABLE]

where $T$ is the temperature and $I_{\mathbf{q}}$ is the effective exchange. The Appendix shows further that the approximations (19) and (20) yield different expressions for $I_{\mathbf{q}}$,

[TABLE]

respectively, where $t$ is the hopping, $\gamma_{\mathbf{q}}$ depends on the dispersion of the lattice, for the square lattice $\gamma_{\mathbf{q}}=\cos(q_{x})+\cos(q_{y})$. In this case $I_{\mathbf{q}}$ is a nearest neighbor interaction, it inherits this property from $\varepsilon_{\mathbf{k}}=-2t\gamma_{\mathbf{k}}$.

It is instructive to evaluate the impurity quantities that determine $I_{\mathbf{q}}$ in the atomic limit where the hybridization function $\Delta$ of DMFT vanishes. Fig. 3 shows that for $U\gg T$ one has $-I_{\mathbf{q}}^{(2)}/\gamma_{\mathbf{q}}\rightarrow 4t^{2}/U$, whereas $-I_{\mathbf{q}}^{(1)}/\gamma_{\mathbf{q}}\rightarrow t^{2}/U$. In combination with equation (21) this implies that only $I^{(2)}$ recovers the correct Néel temperature of the half-filled Hubbard model on the square lattice $[\mathbf{q}=(\pi,\pi)]$ in the Heisenberg mean-field limit, $T_{N}=\frac{4t^{2}}{U}$, whereas $I^{(1)}$ is off by a factor $4$.

Apparently, the vertex corrections at each lattice site need to be treated on an equal footing because the effective exchange is a coupling between equivalent nearest neighbors. Therefore, approximation (19) is used in the applications. The fact that it recovers the effective exchange implies that four-point vertex corrections to the efficient formula (11) can be neglected in the limit $U\gg t,T$, which is confirmed numerically further below.

V.2 Two-particle self-consistency

The DMFT susceptibility $X^{\text{sp}}$ in Eq. (10) may diverge in two dimensions, in violation of the Mermin-Wagner theorem, and it shows the mean-field critical behavior near an antiferromagnetic instability in three dimensions Rohringer et al. (2011); Hirschmeier et al. (2015). As discussed in the context of the two-particle self-consistent (TPSC) approach, these drawbacks are due to the violation of local sum rules Y.M. Vilk and A.-M.S. Tremblay (1997). In order to alleviate the mean-field artifacts a frequency-dependent correction is introduced,

[TABLE]

where $\Pi$ is the DMFT polarization (11). The correction term $\mathcal{U}^{\alpha}_{\omega}$ is fixed by the self-consistency condition,

[TABLE]

thereby $\mathcal{X}$ yields the same kinetic and potential energy as the impurity model of DMFT Krien et al. (2017). Furthermore, the local sum rules are satisfied,

[TABLE]

which are a manifestation of the Pauli principle ($n^{2}_{\sigma}=n_{\sigma}$, cf. Ref. Krien et al. (2017)). Note that $\langle n\rangle=2\sum_{\nu}g_{\nu}$ and $\langle n_{\uparrow}n_{\downarrow}\rangle$ are the density and double occupancy of the impurity model (4).

The boundedness of $\langle n_{\uparrow}n_{\downarrow}\rangle$ in Eq. (27) prevents the divergence of $\mathcal{X}^{\text{sp}}_{q}$ in two dimensions for $T>0$, because it would lead to the logarithmic divergence of the two-dimensional integral $\sum_{\mathbf{q}}$ on the left-hand-side 222 This does not directly imply satisfaction of the Mermin-Wagner theorem, because it has to be shown in practice that a solution $\mathcal{U}_{\omega}$ exists that satisfies Eq. (25).. For dimensions $d>2$ Eq. (27) allows magnetic instabilities for $T>0$, because then the integral $\sum_{\mathbf{q}}$ over the divergent integrand remains finite Y.M. Vilk and A.-M.S. Tremblay (1997). In the limit $d\rightarrow\infty$ the constraint (25) is satisfied by the DMFT susceptibility (10) and hence $\mathcal{U}$ is zero in this limit, as expected. Finally, $\mathcal{X}$ preserves the feature $\imath\omega\mathcal{X}^{\alpha}_{\mathbf{q}=\mathbf{0},\omega}=0$ that is satisfied by the conserving DMFT polarization $\Pi$ in the nominator of Eq. (24). The two-particle spectrum described by $\mathcal{X}$ is therefore ungapped, as required by the global conservation law Krien et al. (2017) 333 Despite the ungapped spectrum the Ward identity is nevertheless violated, because due to the correction $\mathcal{U}$ the static homogeneous limit of $\mathcal{X}$ is inconsistent with the one-particle level of the DMFT approximation Krien et al. (2018). .

For all these reasons the correction $\mathcal{U}_{\omega}$ in Eq. (24) and the constraint (25) appear as suitable in order to remove the mean-field artifacts from the DMFT susceptibility (10). Note that the self-consistency (25) does not lead to a feedback on the impurity model of DMFT, which would in general invalidate the conserving features of the polarization Krien et al. (2017). The correction $\mathcal{U}_{\omega}$ is similar to the constant Moriya-$\lambda$ correction Katanin et al. (2009), it can however not be interpreted straightforwardly as a renormalization of the correlation length, nor is it a retarded interaction. Instead, one may interpret $\mathcal{U}_{\omega}$ as an effective vertex correction to the susceptibility, which takes diagrams beyond DMFT into account that are needed to satisfy the constraint (25). This interpretation is consistent with the TPSC approach Y.M. Vilk and A.-M.S. Tremblay (1997), whose non-perturbative features follow due to effective vertex corrections to the RPA susceptibility. The latter renormalize the mean-field criticality of the RPA Daré et al. (1996). Due to the similarities equation (24) and the constraint (25) are referred to in this work as a two-particle self-consistent dynamical mean-field (TPSC-DMF) approach to the susceptibility.

VI Numerical results

The decomposition of the impurity vertex function in Sec. III, the efficient evaluation of the polarization in Sec. IV, and the TPSC-DMF approach in Sec. V are applied to the two- and three-dimensional Hubbard models (2) at half-filling. In the calculations firstly the DMFT cycle of Sec. II was completed, then the four- and three-point correlation functions (6) and (7) of the AIM were evaluated, where a CTQMC solver based on the ALPS libraries Bauer et al. (2011) with improved estimators Hafermann et al. (2012) was used. The polarization was then evaluated according to Sec. IV, then the TPSC-DMF susceptibility was obtained according to Sec. V. The implementation is based on the dual boson code by E.G.C.P van Loon and H. Hafermann van Loon et al. (2014b).

VI.1 Impurity vertex function

Figure 9 further below shows a phase diagram of the three-dimensional Hubbard model. In this subsection we focus on this model and the values $U/t=6$ and $U/t=14$ of the interaction, which correspond in DMFT to a bad metal and to an insulator with local moments Hirschmeier et al. (2015), respectively, the temperature is set to $T/t\approx 0.4$.

For the metallic regime ($U/t=6$) the left panels of Fig. 4 show the impurity spin vertex function $f^{\,\text{sp}}_{\nu\nu^{\prime}\omega}$ in the static limit $\omega_{0}=0$ and for $\omega_{3}=6\pi T$. In most directions $f^{\,\text{sp}}$ decays with increasing $\nu,\nu^{\prime}$ to a constant, however, it also shows two persistent structures with shapes $+$ and $\times$, see also Ref. Rohringer et al. (2012). For finite $\omega_{3}$ these patterns are shifted along the diagonal.

Important in this work is the exact decomposition $f=f^{i}+\bar{\lambda}^{i}w\lambda^{i}$ discussed in Sec. III. The part that is given by the impurity Hedin vertex $\lambda^{i}$ and by the screened interaction $w$ is shown for $\alpha=\text{sp}$ in the center panels of Fig. 4. This object merely shows a $+$ pattern, while the right panels show the $U^{\text{sp}}$-irreducible vertex $f^{i,\text{sp}}$, which features the $\times$ shape. This correspondence is also there in the charge channel and in different parameter regimes (not shown).

As proven in Sec. IV.1, the irreducible vertex $f^{i}$ does not have a constant background, and the one of the reducible vertex $f$ originates from the term $\bar{\lambda}^{i}w{\lambda}^{i}$. The magnitude of this term is determined by the quantity $[1-U^{\alpha}\pi^{\alpha}_{\omega}]^{-1}$, see equation (17), which can be very large near a quantum critical point $U^{\alpha}\pi^{\alpha}_{\omega=0}\approx 1$, where the impurity susceptibility $\chi^{\alpha}$ is large. The magnitude of $f$ compared to $f^{i}$ at large frequencies therefore depends on the physical regime. For a quantitative comparison Fig. 5 shows the ratio $f^{i}(\nu,\nu^{\prime},\omega)/f(\nu,\nu^{\prime},\omega)$ for fixed $\nu$ and $\omega$ along the $\nu^{\prime}$-direction. The left panels of Fig. 5 show the metallic regime $U/t=6$, where the charge and spin susceptibilities $\chi^{\text{ch}}$ and $\chi^{\text{sp}}$ are both of non-negligible magnitude, and hence $f$ is very large compared to $f^{i}$ at high frequencies. On the other hand, $\chi^{\text{ch}}$ is very small in the insulating regime $U/t=14$, and there is no big difference between $f^{\,\text{ch}}$ and $f^{i,\text{ch}}$, see in particular third panel on the right of Fig. 5. Instead, in this regime the static spin susceptibility $\chi^{\text{sp}}(\omega_{0})$ dominates, leading to the fast decay of $f^{i,\text{sp}}/f^{\,\text{sp}}$ visible in the top right panel.

VI.2 Convergence of frequency summations

Let us observe the faster convergence of Matsubara summations when the irreducible vertex $f^{i}$ is used instead of $f$. For this it is useful to consider the quantity,

[TABLE]

which determines the vertex corrections to the static impurity susceptibility $\chi^{\alpha}_{\omega=0}$ (polarization $\pi^{\alpha}_{\omega=0}$), for finite $\bar{n}$ subjected to a cutoff error. A meaningful measure for convergence is $\epsilon^{(i),\alpha}(\bar{n})=|1-c^{(i),\alpha}(\bar{n}-1)/c^{(i),\alpha}(\bar{n})|$.

Figure 6 shows the ratio $\epsilon^{i}(\bar{n})/\epsilon(\bar{n})$ as function of the cutoff $\bar{n}$ for the cases discussed in Sec. VI.1. Clearly, the summation over $f^{i}$ excels in all cases, having both the numerically smaller error $\epsilon^{i}(\bar{n})<\epsilon(\bar{n})$ and the better scaling with $\bar{n}$. Irregular behavior sets in for large $\bar{n}$ when the Monte Carlo noise exceeds the cutoff error. Surprisingly, the improvement is even sizable in the charge channel for $U/t=14$, where due to the tiny susceptibility $\chi^{\text{ch}}(\omega=0)$ the reducible vertex $f^{\text{ch}}$ has only a small constant background, a worst case scenario. Nevertheless, for $\bar{n}=16$ (i.e., a $32\times 32$ grid) the respective error $\epsilon^{i,\text{ch}}(\bar{n}=16)$ is ten times smaller than $\epsilon^{\text{ch}}(\bar{n}=16)$, see right panel of Fig. 6. In the physically more relevant spin channel this ratio is on the order of one hundred. One should note that in equation (15) summations converge even faster, thanks to the nonlocal bubble $\tilde{X}^{0}$. The improvement of $f^{i}$ over $f$ in the second line of equation (15) is comparable to or better than the example in Fig. 6 (not shown).

VI.3 Two-particle self-consistent susceptibility

The TPSC-DMF susceptibility $\mathcal{X}$ is calculated according to Sec. V. Firstly, it is verified for the Hubbard model on the square lattice that $\mathcal{X}^{\text{sp}}(\mathbf{Q},\omega_{0})$ obeys the exponential scaling with temperature required by the Mermin-Wagner theorem, where $\mathbf{Q}=(\pi,\pi)$. This is shown in the top panel of Fig. 7 for $U/t=8$, at low temperature this corresponds in the DMFT approximation to a strongly correlated Fermi liquid (when paramagnetism is enforced). With increasing $\beta=\frac{1}{T}$ the DMFT susceptibility $X^{\text{sp}}$ quickly diverges, whereas the effective vertex correction $\mathcal{U}$ prevents that the same happens to $\mathcal{X}^{\text{sp}}=[X^{\text{sp},-1}-{\mathcal{U}}/{2}]^{-1}$. For large $\beta$ the correlation length $\xi$ eventually exceeds any fixed system size. Finite-size effects are noticeable when $\xi$ is of order of the half linear system size, then the self-consistent calculation of $\mathcal{U}$ becomes inaccurate (arrows).

The bottom panel of Fig. 7 shows $\mathcal{X}^{\text{sp}}$ in the Brillouin zone for the largest lattice size $256\times 256$ and the lowest considered temperature $T=1/7$. The figure demonstrates simultaneously features of the Mermin-Wagner theorem and of the conservation law: On the one hand the static susceptibility $\mathcal{X}^{\text{sp}}(\mathbf{q}\approx\mathbf{Q},\omega_{0})$ shows the required Lorentzian (Ornstein-Zernike) form Rohringer and Toschi (2016), while on the other hand $\mathcal{X}^{\text{sp}}(\mathbf{q}\approx\mathbf{0},\omega_{1})\propto|\mathbf{q}|^{2}$, which is required by global spin conservation Hafermann et al. (2014).

The top panels of Fig. 8 show the effective vertex correction $\mathcal{U}^{\alpha}(\omega_{m})$ as function of the Matsubara index $m$ and as function of frequency $\omega_{m}$. It is $\mathcal{U}^{\text{sp}}(\omega_{0})>0$, which is required in order for $[-U+\mathcal{U}^{\text{sp}}(\omega_{0})]\Pi^{\text{sp}}(\mathbf{Q},\omega_{0})<1$, preventing the divergence of $\mathcal{X}^{\text{sp}}$ [cf. Eq. (24), note that $\mathcal{X},\Pi<0$]. The temperature dependence of the static spin component $\mathcal{U}^{\text{sp}}(\omega_{0})$ is drawn in the bottom panel of Fig. 8, it is consistent with a smooth crossover from a high temperature regime above the Néel temperature $T_{N}\approx 0.4$ of DMFT into a low temperature regime, which is located roughly below $T_{X}\sim 0.25$. Below this temperature the finite size effects documented in the top panel of Fig. 7 indicate a fast increase of the correlation length, consistent with a renormalized classical regime Y.M. Vilk and A.-M.S. Tremblay (1997). A change in the temperature dependence of $\mathcal{U}^{\text{sp}}(\omega_{0})$ here is plausible, because the momentum integration $\sum_{\mathbf{q}}\mathcal{X}(\mathbf{q},\omega_{0})$ that enters the TPSC-DMF self-consistency (25) is increasingly dominated by the Lorentzian centered at the M point, see bold red line in bottom panel of Fig. 7, whereas at high temperature also other parts of the Brillouin zone contribute. The magnitude of $T_{X}$ corresponds very well to TPSC results at smaller interaction Y.M. Vilk and A.-M.S. Tremblay (1997).

The corrections $\mathcal{U}^{\text{sp}}(\omega>0)$ to the dynamical susceptibility are not affected by $T_{N}$. Indeed, the dashed red line in the bottom panel of Fig. 7 exemplifies that the dynamical susceptibility $\mathcal{X}(\mathbf{q},\omega>0)$ remains flat even far below the Néel temperature of DMFT, which is therefore not a special point. Due to the different temperature dependence of its static and dynamic components $\mathcal{U}^{\text{sp}}(\omega)$ develops a kink and $\mathcal{U}^{\text{sp}}(\omega_{1})-\mathcal{U}^{\text{sp}}(\omega_{0})$ changes sign near $T_{X}$, see bottom panel of Fig. 8. In contrast, $\mathcal{U}^{\text{sp}}(\omega_{2})-\mathcal{U}^{\text{sp}}(\omega_{1})$ is largely independent of temperature over a wide range, although it does show a downturn at very low temperature. Weak temperature dependence of $\mathcal{U}^{\text{sp}}(\omega>0)$ was also observed in the three-dimensional case discussed in the following section.

Also the effective vertex correction $\mathcal{U}^{\text{ch}}(\omega)<0$ of the charge channel is drawn in Fig. 8. The bottom panel shows that its static component is significant only in a region around the Néel temperature of DMFT. Interestingly, it seems therefore that static charge correlators of DMFT, such as the compressibility, remain asymptotically unrenormalized at low temperature. On the other hand, the top left panel of Fig. 8 shows that the dynamic part $\mathcal{U}^{\text{ch}}(\omega>0)$ is mostly on the order of half the Hubbard interaction $U$, indeed a very large correction.

As function of $\omega$ both $\mathcal{U}^{\text{ch}}(\omega)$ and $\mathcal{U}^{\text{sp}}(\omega)$ approach a constant, reminiscent of the Moriya-$\lambda$ correction and of the self-consistent dual boson approach Stepanov et al. (2016b). The sign of these corrections is consistently the opposite of $U^{\text{ch}}=+U$ and $U^{\text{sp}}=-U$, respectively, which may be interpreted as a screening. Due to the frequency dependence of $\mathcal{U}(\omega)$ the criticality of static quantities does not affect dynamic ones. This is different from TPSC and Moriya-$\lambda$, where the same self-consistent correction enters the susceptibility at all frequencies equally.

VI.4 Criticality in three dimensions

A further benchmark for the TPSC-DMF susceptibility is to consider criticality when a spontaneous phase transition is indeed allowed, as is the case in the half-filled three-dimensional Hubbard model. Figure 9 shows the Néel temperature predicted by the ladder dual fermion approach (LDFA) and by the Moriya-$\lambda$-corrected DMFT susceptibility Rohringer et al. (2018). The figure also shows the phase boundary predicted by the TPSC-DMF susceptibility, where $\mathcal{X}^{\text{sp},-1}(\mathbf{Q},\omega_{0})$, $\mathbf{Q}=(\pi,\pi,\pi)$ was fitted with the function $a(T-T_{c})^{-\gamma}$ in order to obtain the critical temperature $T_{c}$ and the critical exponent $\gamma$. The fit interval needs to be bounded from above by the high-$T$ mean-field regime and from below by finite size effects. The upper bound was determined as in Ref. Rohringer et al. (2011), the lower bound is the temperature where the correlation length $\xi$ exceeds $1/6$ of the linear system size of the $16\times 16\times 16$ lattice, as in Ref. Hirschmeier et al. (2015). The boundary obtained by fitting $a,T_{c}$, and $\gamma$ is in excellent agreement with the Moriya-$\lambda$ correction.

The maximum of $T_{c}$ at $U/t=10$ marks the crossover from the bad metal to the insulating regime Hirschmeier et al. (2015). It was found that already at this point the three-dimensional Hubbard model exhibits the Heisenberg universality class Rohringer et al. (2011), where $\gamma\approx 1.4$. Consistent with this the fit of $\mathcal{X}$ for $U/t\geq 10$ yields an exponent of roughly $1.35$, which compares to the mean-field exponent $1$ of DMFT. In this regime $T_{c}$ was also estimated with $\gamma$ assumed to be known from the Heisenberg model, see blue circles in Fig. 9, which leads to an even better agreement with the Moriya-$\lambda$ correction, it therefore seems that the TPSC-DMF approach predicts the same critical behavior 444 The similar results are a consequence of similar self-consistency conditions. The Moriya-$\lambda$ is fixed by the local sum rules (26) and (27), whose left-hand-sides are in general dominated by the static term $\omega=0$ near a phase transition, in this case $\mathcal{U}^{\alpha}(\omega=0)\approx\lambda^{\alpha}_{\text{Moriya}}$. .

Lastly, $T_{c}$ was also determined when vertex corrections to the Hedin vertex are neglected, $\Lambda^{i}\approx\lambda^{i}$. This approximation is applied as in equation (19), $\Pi\approx\Pi^{(2)}$, for the reasons explained in Sec. V.1. Note that once again the constraint $\mathcal{X}_{\text{loc}}=\chi$ is satisfied by self-consistent adjustment of $\mathcal{U}$ in Eq. (24). In fact, also this approximation clearly deviates from the mean-field criticality near the transition and for $U/t\geq 10$ is well-described by the Heisenberg critical exponent, as shown in the right panel of Fig. 9. Without vertex corrections $T_{c}$ lies reasonably close to the result with the vertex corrections (left panel, yellow and red lines) but the deviation depends on the physical regime. For large coupling the vertex corrections have negligible influence on $T_{c}$, which confirms the analytical result of Sec. V.1, but they play an important role in the region where DMFT predicts a bad metal.

VII Conclusions

An efficient method to evaluate the DMFT susceptibility was presented by making use of the Hedin three-leg vertex. Vertex corrections to the latter arise in the form of a four-point vertex $f^{i}$ of the Anderson impurity model that is irreducible with respect to the bare interaction $\pm U$. This vertex has no constant background, in contrast to the full impurity vertex $f$. Furthermore, the ladder equation for the Hedin vertex is formulated in terms of nonlocal Green’s functions, as in the dual fermion approach Rubtsov et al. (2008). The combination of the fast decay of the nonlocal Green’s functions with the decay of the irreducible vertex $f^{i}$ leads to a faster convergence of frequency summations compared to the dual fermion and dual boson approaches Rubtsov et al. (2012). As a result, the measurement of the four-point vertex can be restricted to a smaller frequency window. The efficient calculation scheme can be generalized to multi-orbital Hubbard models and symmetry-broken phases (see Appendix E), furthermore, it may be possible to incorporate it into the dual fermion and dual boson formalisms Rubtsov et al. (2008, 2012); Stepanov et al. (2016a).

The efficient calculation scheme implicitly takes vertex asymptotics into account, which were discussed, for example, in Refs. Kuneš (2011); Wentzell et al. (2016); Kaufmann et al. (2017); Tagliavini et al. (2018). In the implementation it is nevertheless not necessary to consider the large frequency limits explicitly, because the contributions to the DMFT susceptibility that originate from the constant background of the reducible vertex $f$ are handled in an exact way. The main difference to the previously presented approaches to reduce the cutoff error by taking vertex asymptotics into account is that a diagrammatic decomposition of $f$ is employed that is exact for all frequencies, leading to a particularly simple calculation scheme. The cutoff error may be reduced further by taking the asymptotic behavior of the irreducible vertex $f^{i}$ into account.

The mean-field instability of the DMFT susceptibility was removed by introduction of a frequency-dependent correction $\mathcal{U}(\omega)$ that is fixed by adjusting the local susceptibility to the impurity, $\mathcal{X}_{\text{loc}}(\omega)=\chi(\omega)$. This approach ensures an ungapped two-particle spectrum and the expected critical behavior in two dimensions in agreement with the Mermin-Wagner theorem, reminiscent of the two-particle self-consistent (TPSC) approach that is based on the Hartree/RPA approximation Y.M. Vilk and A.-M.S. Tremblay (1997). Indeed, the temperature dependence of $\mathcal{U}(\omega=0)$ shows a crossover to a renormalized classical regime, a hallmark effect of the TPSC approach Y.M. Vilk and A.-M.S. Tremblay (1997). In the half-filled three-dimensional Hubbard model the criticality of the approach is consistent with the similar Moriya-$\lambda$ correction used in the dynamical vertex approximation Katanin et al. (2009), which leads to a renormalized correlation length.

The interpretation of $\mathcal{U}(\omega)$ is however different as a somewhat intransparent vertex correction beyond DMFT, it is therefore necessary to consider the domain of validity of the approach: To do this for the weak coupling limit, one may recall that the TPSC approach requires that the Hartree approximation provides a reasonable description of the Fermi surface nesting Y.M. Vilk and A.-M.S. Tremblay (1997). However, in the half-filled two-dimensional Hubbard model on the square lattice a pseudogap opens at low temperature due to antiferromagnetic fluctuations Vilk and Tremblay (1996); Schäfer et al. (2015); van Loon et al. (2018b). In this case neither the Hartree approximation nor DMFT provide a good starting point, because they predict a homogeneous Fermi surface with strong nesting. On the other hand, even when the feedback of the pseudogap on the two-particle spectrum is taken into account it leads to similar results as the Moriya-$\lambda$-corrected DMFT susceptibility Tanaka (2019); Rohringer et al. (2018). In the large coupling limit the self-consistency $\mathcal{X}_{\text{loc}}=\chi$ imposes the unscreened local moment of a Mott insulator by construction, although in reality it may be screened due to short-ranged correlations. Two-particle self-consistency can therefore impose a bias towards the physics of the impurity model, furthermore, when it makes a feedback on the impurity model it can violate conservation laws Krien et al. (2017), which was therefore avoided. In the future it may be investigated whether the $\mathcal{U}(\omega)$ correction yields a similar feedback on the single-particle spectrum as the Moriya-$\lambda$ correction Rohringer and Toschi (2016) and whether it can be generalized to the multi-orbital case Galler et al. (2017). A further perspective is to consider the effect of the frequency dependence of $\mathcal{U}(\omega)$ on the two-particle spectrum.

Finally, it was shown that for large coupling vertex corrections to the Hedin vertex play a minor role for the Néel temperature of the half-filled three-dimensional Hubbard model. This strengthens the case for a local approximation to the Hedin vertex in this regime, as in the TRILEX approach Ayral and Parcollet (2016). However, it was found that at the level of DMFT the polarization diagram of TRILEX underestimates the prefactor of the effective exchange with energy scale $t^{2}/U$. The correct prefactor is obtained when the local approximation to the Hedin vertex is applied to the efficient formula for the polarization, which corresponds to the dual boson approach Stepanov et al. (2016a). This formula treats vertex corrections at each lattice site on an equal footing.

During the completion of this work a manuscript was preprinted Otsuki et al. (2019) that derives a strong coupling form of the DMFT spin susceptibility with an effective exchange cutoff. Here this quantity was expressed in terms of local Hedin vertex and polarization of the impurity model. The latter remain finite at zero temperature, the effective exchange is therefore well-defined in this limit. The calculation of the spin susceptibility in the Mott phase at zero temperature is an unsolved problem Georges et al. (1996); Guerci et al. (2018); Krien et al. (2018).

Acknowledgements.

I thank the anonymous referees for constructive comments that improved this work. F.K. thanks E.G.C.P. van Loon and A. Valli for their reading of the manuscript, and A.I. Lichtenstein, E.G.C.P. van Loon, and H. Hafermann for long-time support, and M. Capone, A. Toschi, K. Held, E.A. Stepanov, J. Otsuki, A. Katanin, and L. Fanfarillo for fruitful discussions. F.K. acknowledges support by MIUR through the PRIN 2015 program (Prot. 2015C5SEJJ001) and the SISSA/CNR project ”Superconductivity, Ferroelectricity and Magnetism in bad metals” (Prot. 232/2015).

Appendix A $U^{\alpha}$-irreducible vertices

It is shown how diagrams that are reducible with respect to the bare interaction $U^{\alpha}$ can be separated from the three-leg vertex $\Lambda$ and from the vertex function $F$, following an approach of Hertz and Edwards Hertz and Edwards (1973). The relations in this section of the appendix are formally exact for the paramagnetic Hubbard model (2), for the Anderson impurity model (4) capital letters may be replaced by small letters ($\Lambda\rightarrow\lambda,F\rightarrow f$, and so on) and four-momenta are replaced by frequencies [$k=(\mathbf{k},\nu)\rightarrow\nu,q=(\mathbf{q},\omega)\rightarrow\omega$]. Generalizations to more general lattice and impurity models are briefly discussed in Appendix E.

A.1 Correlation functions

The four-point function is defined as,

[TABLE]

where definitions are as in the main text. It is convenient to define the generalized susceptibility,

[TABLE]

the latter can be represented in terms of a ladder equation $\hat{X}=\hat{X}^{0}+\hat{X}^{0}\hat{\Gamma}\hat{X}$, where $\Gamma$ is the two-particle self-energy and all quantities denote matrices in the labels $k,k^{\prime}$ and $X^{0}_{kk^{\prime}}=N\beta G_{k}G_{k+q}\delta_{kk^{\prime}}$ is the bubble. Matrix multiplication implies a factor $(N\beta)^{-1}$, the labels $q,\alpha$ are suppressed.

A.2 $U^{\alpha}$-irreducible generalized susceptibility

The goal is to separate the diagrams from $\hat{X}$ that are reducible with respect to $U^{\text{ch}}=+U$ and $U^{\text{sp}}=-U$, respectively. To this end, one defines $\hat{\Gamma}^{i}=\hat{\Gamma}-\hat{\Gamma}^{0}$, where ${\Gamma}^{0,\alpha}_{kk^{\prime}}=U^{\alpha}$ is the bare two-particle self-energy. The ladder equation for $\hat{X}$ can therefore be written as,

[TABLE]

which implies super-matrix inversion with respect to $k,k^{\prime}$. Let us now define the $\Gamma^{0}$-irreducible generalized susceptibility $\hat{\Pi}$,

[TABLE]

There are no diagrams in $\Pi$ that can be separated into two parts by removing a single vertex $\Gamma^{0}$ [in the sense of Fig. 10]. Subtracting Eq. (30) from (29) eliminates $\Gamma^{i}$ and $X^{0}$,

[TABLE]

In explicit notation this relation simplifies (the label $\alpha$ remains dropped),

[TABLE]

where $\Gamma^{0}=\pm U$, summations imply $(N\beta)^{-1}$.

A.3 Three-leg vertices and polarization

$X$ and $\Pi$ will now be related to the left- and right-sided three-leg vertices $\Lambda^{(i)}$ and $\bar{\Lambda}^{(i)}$, using the definitions,

[TABLE]

where in the second line the $\Gamma^{0}$-irreducible (Hedin) three-leg vertex $\Lambda^{i}$ was introduced and $X^{0}_{kq}=G_{k+q}G_{k}$ is the bubble. The reducible and irreducible three-leg vertices are related via Eq. (32), which is seen by summation over $k^{\prime}$,

[TABLE]

Finally, dividing by $X^{0}_{kq}$ and identifying the susceptibility, $X_{q}=2\sum_{kk^{\prime}}X_{kk^{\prime}q}$, one arrives at the simple relation,

[TABLE]

where the label $\alpha$ was reintroduced. In the second line the polarization was defined,

[TABLE]

By summing Eq. (35) over $k$ one sees that,

[TABLE]

Note that in contrast to the susceptibility $X_{q}$ a factor $2$ does not occur [see above Eq. (37)]. Similar to Eq. (38) one derives in an analogous way the relation for the left-sided three-leg vertex,

[TABLE]

A.4 Four-leg vertices and screened interaction

Next, also the vertex function $F$ will be expressed in terms of a $\Gamma^{0}$-irreducible counterpart $F^{i}$. To do this, the following relation between the generalized susceptibility $X$ and $F$ will be used,

[TABLE]

Inserting these relations into Eq. (32), and using once again Eqs. (33) and (34) leads to,

[TABLE]

Finally, dividing by $X^{0}_{kq}X^{0}_{k^{\prime}q}$ and using Eq. (41), the reducible vertex $F$ can be expressed in terms of the irreducible vertices $F^{i}$ and $\Lambda^{i}$,

[TABLE]

where the label $\alpha$ was reintroduced and the screened interaction $W$ is defined as,

[TABLE]

For the impurity model one makes in Eqs. (41), (44), and (45) the replacements $F\rightarrow f$, $\Lambda\rightarrow\lambda$, $W\rightarrow w$, and $\Pi\rightarrow\pi$, leading to Eqs. (8) and (9) in the main text.

Appendix B Ladder equation for the reducible three-leg vertex

Ladder equations for the reducible and irreducible three-leg vertices $\Lambda$ and $\Lambda^{i}$ are derived in the DMFT approximation, where the two-particle self-energy is approximated with the one of the impurity model (4), $\Gamma^{\alpha}_{kk^{\prime}q}=\gamma^{\alpha}_{\nu\nu^{\prime}\omega}$ Georges et al. (1996); Krien et al. (2017). In this case the Bethe-Salpeter equation for the lattice vertex function $F$ reads,

[TABLE]

where it was used that for a local two-particle self-energy $\Gamma$ the vertex function does not depend on the momenta $\mathbf{k},\mathbf{k}^{\prime}$. $X^{0}_{\nu}(q)=\sum_{\mathbf{k}}G_{k}G_{k+q}$ denotes the bubble of DMFT Green’s functions (3).

By $\nu,\nu^{\prime}$-matrix inversion one obtains from Eq. (46) in a short notation, $\hat{\gamma}^{\alpha,-1}_{\omega}=\hat{F}^{\alpha,-1}_{q}+\hat{X}^{0}(q)$, where $X^{0}_{\nu\nu^{\prime}}(q)=\beta X^{0}_{\nu}(q)\delta_{\nu\nu^{\prime}}$. Similarly, there exists an impurity Bethe-Salpeter equation, $\hat{\gamma}^{\alpha,-1}_{\omega}=\hat{f}^{\alpha,-1}_{\omega}+\hat{\chi}^{0}(\omega)$, where $f$ denotes the impurity vertex function and $\chi^{0}_{\nu\nu^{\prime}}(\omega)=\beta g_{\nu}g_{\nu+\omega}\delta_{\nu\nu^{\prime}}$. Thereby, $\gamma$ is eliminated in favor of $f$, leading to the exact reformulation of Eq. (46),

[TABLE]

where $\tilde{X}^{0}_{\nu}(q)=\sum_{\mathbf{k}}(G_{k}G_{k+q}-g_{\nu}g_{\nu+\omega})$ is the nonlocal bubble, see also Ref. Hafermann et al. (2014).

In order to arrive at an analogous ladder equation for the three-leg vertex $\Lambda$, Eq. (47) is multiplied by $G_{k^{\prime}}G_{k^{\prime}+q}$, summed over $k^{\prime}$, and $1$ is added on both sides,

[TABLE]

On the left-hand-side (LHS) arises the right-sided three-leg vertex, $\bar{\Lambda}_{\nu q}=1+\sum_{k^{\prime}}F_{\nu\nu^{\prime}q}G_{k^{\prime}}G_{k^{\prime}+q}$, on the right-hand-side (RHS) $\sum_{\mathbf{k}}G_{k}G_{k+q}=\tilde{X}^{0}_{\nu}(q)+g_{\nu}g_{\nu+\omega}$ is inserted,

[TABLE]

On the RHS one identifies the right-sided impurity three-leg vertex $\bar{\lambda}_{\nu\omega}=1+\sum_{\nu^{\prime}}f_{\nu\nu^{\prime}\omega}g_{\nu^{\prime}}g_{\nu^{\prime}+\omega}$, and $\sum_{\nu^{\prime\prime}}f^{\alpha}_{\nu\nu^{\prime\prime}\omega}\tilde{X}^{0}_{\nu^{\prime\prime}}(q)$ is factored out,

[TABLE]

The term in brackets is again $\bar{\Lambda}$, leading to the ladder equation for the right-sided three-leg vertex,

[TABLE]

The analogous ladder equation for the left-sided three-leg vertex $\Lambda$ follows from the symmetry of the impurity vertex, $f_{\nu\nu^{\prime}\omega}=f_{\nu^{\prime}+\omega,\nu+\omega,-\omega}$,

[TABLE]

Appendix C Efficient formulae for susceptibility and polarization

Efficient formulae for the susceptibility and polarization are derived. The susceptibility may be calculated from the reducible three-leg vertex $\Lambda$ as,

[TABLE]

In the DMFT approximation $\Lambda$ does not depend on $\mathbf{k}$, hence, $X^{\alpha}_{q}=2\sum_{\nu}\Lambda^{\alpha}_{\nu q}X^{0}_{\nu}(q)$, where $X^{0}_{\nu}(q)=\sum_{\mathbf{k}}G_{k}G_{k+q}$. This relation will be rewritten as the sum of impurity susceptibility $\chi$ and nonlocal corrections $\tilde{X}$.

To do this, the bubble $X^{0}$ is expressed in terms of the nonlocal bubble $\tilde{X}^{0}$ and the impurity bubble $g_{\nu}g_{\nu+\omega}$, $X^{0}_{\nu}(q)=\tilde{X}^{0}_{\nu}(q)+g_{\nu}g_{\nu+\omega}$, furthermore, Eq. (52) is substituted for the three-leg vertex $\Lambda$,

[TABLE]

Four terms arise, the impurity susceptibility can be identified, $\chi^{\alpha}_{\omega}=2\sum_{\nu}\lambda^{\alpha}_{\nu\omega}g_{\nu}g_{\nu+\omega}$. Furthermore,

[TABLE]

where the right-sided impurity three-leg vertex was identified, $\bar{\lambda}^{\alpha}_{\nu^{\prime}\omega}=1+\sum_{\nu}f^{\alpha}_{\nu^{\prime}\nu\omega}g_{\nu}g_{\nu+\omega}$ [its trivial part $1$ is canceled by the second term on the RHS of Eq. (55)]. Using these relations in Eq. (54) leads to,

[TABLE]

Using the ladder equation (52) for $\Lambda$ it is seen that the second line cancels the third, hence,

[TABLE]

which is the dual boson formula (1) Rubtsov et al. (2012); Hafermann et al. (2014).

A similar relation will be derived for the polarization $\Pi$. To do this, let us invoke the local analogue of Eq. (38),

[TABLE]

where $\bar{\lambda},\bar{\lambda}^{i}$, and $\pi$ are the three-leg vertices and the polarization of the impurity. The latter is related to $\chi$ analogous to Eq. (39),

[TABLE]

Using Eqs. (58), (59) for the impurity quantities, and Eqs. (41), (39) for the lattice quantities in Eq. (57) leads to,

[TABLE]

Multiplication by $1-U^{\alpha}\Pi^{\alpha}_{q}$ and ${1-U^{\alpha}\pi^{\alpha}_{\omega}}$ leads to the desired relation (11) for the polarization,

[TABLE]

Again, compared to Eq. (57) a factor $2$ does not occur.

Appendix D Ladder equation for the Hedin vertex

Equation (52) is now reformulated for the $U^{\alpha}$-irreducible three-leg vertex $\Lambda^{i}$. To do this, Eq. (41) and its local analogue ${\lambda}^{\alpha}_{\nu\omega}={{\lambda}^{i,\alpha}_{\nu\omega}}/({1-U^{\alpha}\pi^{\alpha}_{\omega}})$ are inserted into the ladder equation (52) for $\Lambda$,

[TABLE]

both sides were multiplied by a factor $1-U^{\alpha}\Pi^{\alpha}_{q}$. On the RHS appears the reducible impurity vertex function $f$, which will be eliminated in favor of its irreducible counterpart $f^{i}$ using the local analogue of Eq. (44),

[TABLE]

where $w$ is the screened interaction of the impurity,

[TABLE]

Inserting Eq. (63) into Eq. (62) leads to,

[TABLE]

Using Eqs. (64) and (45) the fraction on the RHS can be expressed as $\frac{w^{\alpha}_{\omega}}{W^{\alpha}_{q}}$. Furthermore, Eq. (61) can be used to identify in the second line, $\sum_{\nu^{\prime}}\Lambda^{i,\alpha}_{\nu^{\prime}q}\tilde{X}^{0}_{\nu^{\prime}}(q)\bar{\lambda}^{i,\alpha}_{\nu^{\prime}\omega}=\Pi^{\alpha}_{q}-\pi^{\alpha}_{\omega}$. Eq. (65) thus becomes,

[TABLE]

Using the relation (64) between $w$ and $\pi$, and the relation (45) between $W$ and $\Pi$ leads to the desired ladder equation (13) for the Hedin vertex,

[TABLE]

Appendix E General bare interaction

In the Hedin formalism the bosons arise because Green’s function lines are contracted at a bare interaction vertex that does not depend on fermionic momentum-energies $k=(\mathbf{k},\nu)$, see Sec. III.2. This requirement allows for much more general interaction Hamiltonians than considered here.

In particular, the Appendices A-D (i.e., the efficient calculation of the DMFT polarization) can be generalized to multi-orbital systems and/or symmetry-broken phases. In these cases a matrix-valued bare interaction of the form $U^{ab}$ enters the Bethe-Salpeter equation, where $a=(m_{1}m_{2}\sigma_{1}\sigma_{2})$ is a superindex of two orbital and two spin indices 555 For the single-band Hubbard model $U^{\uparrow\uparrow\downarrow\downarrow}=U^{\downarrow\downarrow\uparrow\uparrow}=-U^{\downarrow\uparrow\uparrow\downarrow}=-U^{\uparrow\downarrow\downarrow\uparrow}=U$, the other elements are zero. , see also Ref. Galler et al. (2017). As in Appendix A.2 one removes the bare part from the two-particle self-energy, $\hat{\Gamma}^{i,ab}=\hat{\Gamma}^{ab}-\hat{\Gamma}^{0,ab}$, where $\Gamma^{0,ab}_{kk^{\prime}}=U^{ab}$. One then derives the crucial equation (32), which becomes a matrix relation with respect to the superindices, it serves as the vantage point for the remaining calculations.

On an equal footing it seems possible to introduce a TPSC-DMF prescription (25), $\sum_{\mathbf{q}}\mathcal{X}^{ab}_{q}=\chi^{ab}_{\omega}$, which is fixed by an effective vertex correction $\mathcal{U}^{ab}(\omega)$ 666 It is unclear whether a generalization of TPSC-DMF to symmetry-broken phases inherits thermodynamic consistency at second order critical points from the DMFT approximation Krien (2018), which may be clarified in future work. .

Finally, it is possible to generalize Appendix A to a nonlocal and/or retarded interaction. However, only the RPA-like vertex $U(\mathbf{q},\omega)$ can be separated from the Bethe-Salpeter equation, not the Fock exchange $U(\mathbf{k}^{\prime}-\mathbf{k},\nu^{\prime}-\nu)$, since it depends on the fermionic variables. Appendices B-D rely on the DMFT approximation where interaction of lattice and impurity need to be equivalent.

Appendix F Strong coupling limit

This appendix considers phase transitions of the half-filled Hubbard model in the strong coupling limit $U\gg T,t$. Static impurity quantities carry a label ‘[math]’, e.g., $\pi(\omega=0)=\pi_{0}$, furthermore $q_{0}=(\mathbf{q},\omega=0)$.

F.1 DMFT

Near an instability of the static DMFT spin susceptibility $X^{\text{sp}}(q_{0})=2/[\Pi^{\text{sp},-1}(q_{0})-U^{\text{sp}}]$ one has for the polarization $\Pi^{\text{sp}}(q_{0})=\pi^{\text{sp}}_{0}+\tilde{\Pi}^{\text{sp}}(q_{0})\approx\frac{1}{U^{\text{sp}}}$. On the other hand, for strong coupling and at half-filling DMFT predicts a Mott insulator with $\chi^{\text{sp}}_{0}\propto-\beta$. Using $\pi^{\text{sp}}_{0}=\frac{1}{2}\chi^{\text{sp}}_{0}/[1+U^{\text{sp}}\frac{1}{2}\chi^{\text{sp}}_{0}]$ and $U\beta\gg 1$ it follows that $\pi^{\text{sp}}_{0}\approx\frac{1}{U^{\text{sp}}}$. Hence, $|\pi^{\text{sp}}_{0}|\gg|\tilde{\Pi}^{\text{sp}}(q_{0})|$ and one can expand,

[TABLE]

In the second line $\chi^{\text{sp}}_{0}=2/(\pi^{\text{sp},-1}_{0}-U^{\text{sp}})$ was used. Defining the effective exchange as $I_{\mathbf{q}}=-\pi^{\text{sp},-2}_{0}\tilde{\Pi}^{\text{sp}}(q_{0})$ one arrives at the strong coupling form of the DMFT spin susceptibility Otsuki et al. (2019). For very large interaction the local moment is fully developed and $\chi^{\text{sp},-1}_{0}\approx-T$, leading to equation (21) in the main text.

F.2 TRILEX-like approximation

Let us consider the approximations (19) and (20) for the polarization $\Pi$. Both expressions contain the nonlocal bubble $\tilde{X}^{0}_{\nu}(q)=\sum_{\mathbf{k}}G_{k}G_{k+q}-g_{\nu}g_{\nu+\omega}$, which can be simplified in the strong coupling limit using similar steps as in Ref. Otsuki et al. (2019). For small hybridization $\Delta\approx 0$ one can expand Green’s function $G_{k}\approx g_{\nu}+g_{\nu}\varepsilon_{\mathbf{k}}g_{\nu}$, hence,

[TABLE]

where it was used that $\sum_{\mathbf{k}}\varepsilon_{\mathbf{k}}=\sum_{\mathbf{k}}\varepsilon_{\mathbf{k}+\mathbf{q}}=0$. For the dispersion $\varepsilon_{\mathbf{k}}=-2t\gamma_{\mathbf{k}}$ of the $d$-dimensional hypercubic lattice with $\gamma_{\mathbf{k}}=\sum_{i=1}^{d}\cos(k_{i})$ one has $\sum_{\mathbf{k}}\varepsilon_{\mathbf{k}}\varepsilon_{\mathbf{k}+\mathbf{q}}=2t^{2}\gamma_{\mathbf{q}}$. Using this and Eq. (70) yields for the nonlocal part of (19) and (20), respectively,

[TABLE]

Inserting into Eq. (69) leads to the expressions $I^{(2)}$ and $I^{(1)}$ for the effective exchange in equations (22) and (23).