Strongly correlated electrons: Analytic mean-field theories with two-particle self-consistency

TL;DR

This paper develops an analytic mean-field framework incorporating two-particle self-consistency via parquet equations, enabling accurate modeling of strongly correlated electrons and avoiding unphysical artifacts.

Contribution

It introduces a general scheme for self-consistent approximations using parquet equations, ensuring quantum criticality is preserved in strong coupling regimes.

Findings

Applied to the single-impurity Anderson model demonstrating Kondo behavior.

Ensured compliance with Ward identity and Schwinger-Dyson equation.

Provided a static approximation capturing essential physics of strong correlations.

Abstract

A two-particle self-consistency is rarely part of mean-field theories. It is, however, essential for avoiding spurious critical transitions and unphysical behavior. We present a general scheme for constructing analytically controllable approximations with self-consistent equations for the two-particle vertices based on the parquet equations. We explain in detail how to reduce the full set of parquet equations not to miss quantum criticality in strong coupling. We further introduce a decoupling of convolutions of the dynamical variables in the Bethe-Salpeter equations to make them analytically solvable. We connect the self-energy with the two-particle vertices to satisfy the Ward identity and the Schwinger-Dyson equation. We discuss the role of the one-particle self-consistency in making the approximations reliable in the whole spectrum of the input parameters. We demonstrate the general construction on the simplest static approximation that we apply to the Kondo behavior of the single-impurity Anderson model.

Tables

Table 1. Table 1: Notation for 2P vertices used here and in our preceding papers compared to that used in recent papers on parquet equations as in review  Rohringer et al. ( 2018 ) , Ref.  Bickers ( 1991 ) , and one of the first papers on parquet equations in condensed mater  Roulet et al. ( 1969 ) . We skipped the spin indices but kept the channel index α = e ​ h , e ​ e , e ​ h ¯ 𝛼 𝑒 ℎ 𝑒 𝑒 ¯ 𝑒 ℎ \alpha=eh,ee,\overline{eh} .

Source	Completely	Channel	Channel	Full
	irreducible	irreducible	reducible	vertex
This paper	$I$	${\Lambda }^{\alpha }$	$K^{\alpha }$	$\Gamma$
Ref. Rohringer et al. (2018)	$\Lambda$	${\Gamma }^{\alpha }$	${\Phi }^{\alpha }$	$F$
Ref. Bickers (1991)	${\Lambda }^{irr}$	${\Gamma }_{\alpha }$	–	$\Gamma$
Ref. Roulet et al. (1969)	$R$	$I^{\alpha }$	${\gamma }^{\alpha }$	$\Gamma$

Equations

$$\begin{gathered} {\Gamma }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)={\Lambda }_{\sigma \bar{\sigma }}^{eh}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)-\frac{1}{N}\sum _{{\mathbf{k}}^{\prime \prime }}\frac{1}{\beta }\sum _{{\omega }_l}{\Lambda }_{\sigma \bar{\sigma }}^{eh}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime \prime },i{\omega }_l;\mathbf{q},i{\nu }_m) \\ \times G_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l)G_{\bar{\sigma }}({\mathbf{k}}^{\prime \prime }+\mathbf{q},i{\omega }_{m+l}){\Gamma }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m), \end{gathered}$$

$$\begin{gathered} {\Gamma }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)={\Lambda }_{\sigma \bar{\sigma }}^{ee}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)-\frac{1}{N}\sum _{{\mathbf{k}}^{\prime \prime }}\frac{1}{\beta }\sum _{{\omega }_l}{\Lambda }_{\sigma \bar{\sigma }}^{ee}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime \prime },i{\omega }_l;\mathbf{q}+{\mathbf{k}}^{\prime }-{\mathbf{k}}^{\prime \prime },i{\nu }_{m+n^{\prime }-l}) \\ \times G_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l)G_{\bar{\sigma }}(\mathbf{q}+\mathbf{k}+{\mathbf{k}}^{\prime }-{\mathbf{k}}^{\prime \prime },i{\omega }_{m+n+n^{\prime }-l}){\Gamma }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q}+\mathbf{k}-{\mathbf{k}}^{\prime \prime },i{\nu }_{m+n-l}). \end{gathered}$$

$$\begin{gathered} {\Gamma }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)={\Lambda }_{\sigma \bar{\sigma }}^{eh}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)+{\Lambda }_{\sigma \bar{\sigma }}^{ee}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m) \\ -I_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m). \end{gathered}$$

$$\begin{gathered} {\Lambda }_{\sigma \bar{\sigma }}^{eh}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)=U-\frac{1}{N}\sum _{\mathbf{Q}}\frac{1}{\beta }\sum _{{\nu }_l}{\Lambda }_{\sigma \bar{\sigma }}^{ee}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime }+\mathbf{q}+\mathbf{Q},i{\omega }_{n^{\prime }+m+l};-\mathbf{Q}-i{\nu }_l) \\ \times G_{\sigma }({\mathbf{k}}^{\prime }+\mathbf{q}+\mathbf{Q},i{\omega }_{n^{\prime }+m+l})G_{\bar{\sigma }}(\mathbf{k}-\mathbf{Q},i{\omega }_{n-l}){\Gamma }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime }+\mathbf{q}+\mathbf{Q},i{\omega }_{n^{\prime }+m+l},{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{k}-{\mathbf{k}}^{\prime }-\mathbf{Q},i{\nu }_{n-n^{\prime }-l}), \end{gathered}$$

$$\begin{gathered} {\Lambda }_{\sigma \bar{\sigma }}^{ee}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime \prime },i{\omega }_l;{\mathbf{q}}^{\prime }-\mathbf{k}-{\mathbf{k}}^{\prime \prime },i{\nu }_{m^{\prime }-n-l})G_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l)G_{\bar{\sigma }}({\mathbf{q}}^{\prime }-{\mathbf{k}}^{\prime \prime },i{\omega }_{m^{\prime }-l}){\Gamma }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};{\mathbf{q}}^{\prime }-{\mathbf{k}}^{\prime }-{\mathbf{k}}^{\prime \prime },i{\nu }_{m^{\prime }-n^{\prime }-l}) \\ \to K_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime \prime },i{\omega }_l;{\mathbf{q}}^{\prime }-\mathbf{k}-{\mathbf{k}}^{\prime \prime },i{\nu }_{m^{\prime }-n-l})G_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l)G_{\bar{\sigma }}({\mathbf{q}}^{\prime }-{\mathbf{k}}^{\prime \prime },i{\omega }_{m^{\prime }-l}){\Lambda }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};{\mathbf{q}}^{\prime }-{\mathbf{k}}^{\prime }-{\mathbf{k}}^{\prime \prime },i{\nu }_{m^{\prime }-n^{\prime }-l}), \end{gathered}$$

$$\begin{gathered} {\Lambda }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n;{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }})=U-\frac{1}{N}\sum _{\mathbf{Q}}\frac{1}{\beta }\sum _{{\nu }_l}{K}_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime }+\mathbf{Q},i{\omega }_{n^{\prime }}+i{\nu }_l;-\mathbf{Q},-i{\nu }_l) \\ \times G_{\sigma }({\mathbf{k}}^{\prime }+\mathbf{Q},i{\omega }_{n^{\prime }}+i{\nu }_l)G_{\bar{\sigma }}(\mathbf{k}-\mathbf{Q},i{\omega }_n-i{\nu }_l){\Lambda }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime }+\mathbf{Q},i{\omega }_{n^{\prime }}+i{\nu }_l,\mathbf{k}-\mathbf{Q},i{\omega }_n-i{\nu }_l). \end{gathered}$$

$$\begin{gathered} K_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)=-\frac{1}{N}\sum _{{\mathbf{k}}^{\prime \prime }}\frac{1}{\beta }\sum _{{\omega }_l}{\Lambda }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n;\mathbf{q}+{\mathbf{k}}^{\prime \prime },i{\omega }_{m+l})G_{\bar{\sigma }}({\mathbf{k}}^{\prime \prime }+\mathbf{q},i{\omega }_{m+l}) \\ \times G_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l)\left[K_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)+{\Lambda }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l;\mathbf{q}+{\mathbf{k}}^{\prime },i{\omega }_{m+n^{\prime }})\right] \end{gathered}$$

$$\begin{gathered} {\Lambda }_{\uparrow \downarrow }(\mathbf{k},i{\omega }_n;{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }}) \\ =\frac{U}{1+N^{-1}{\sum }_{\mathbf{Q}}{\beta }^{-1}{\sum }_{{\nu }_m}{K}_{\uparrow \downarrow }({\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }},\mathbf{k},i{\omega }_n;-\mathbf{Q},-i{\nu }_m)G_{\uparrow }(\mathbf{k}+\mathbf{Q},i{\omega }_{n+m})G_{\downarrow }({\mathbf{k}}^{\prime }-\mathbf{Q},i{\omega }_{n^{\prime }-m})}. \end{gathered}$$

$$\begin{gathered} \frac{1}{N}\sum _{{\mathbf{k}}^{\prime \prime }}\frac{1}{\beta }\sum _{{\omega }_l}\left[\beta N{\delta }_{\mathbf{k},{\mathbf{k}}^{\prime \prime }}{\delta }_{n,l}+{\Lambda }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime \prime }+\mathbf{q},i{\nu }_m+i{\omega }_l)G_{\bar{\sigma }}(\mathbf{q}+{\mathbf{k}}^{\prime \prime },i{\nu }_m+i{\omega }_l)G_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l)\right] \\ \times K_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)=-\frac{1}{N}\sum _{{\mathbf{k}}^{\prime \prime }}\frac{1}{\beta }\sum _{{\omega }_l}{\Lambda }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n;{\mathbf{k}}^{\prime \prime },i{\omega }_l)G_{\bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l)G_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l){\Lambda }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }}). \end{gathered}$$

$$\begin{gathered} {\Gamma }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m)={\Lambda }_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime }+\mathbf{q},i{\omega }_{n^{\prime }+m}) \\ +K_{\sigma \bar{\sigma }}(\mathbf{k},i{\omega }_n,{\mathbf{k}}^{\prime },i{\omega }_{n^{\prime }};\mathbf{q},i{\nu }_m). \end{gathered}$$

$$\begin{gathered} \Delta \Sigma (\mathbf{k},i{\omega }_n)=-\frac{1}{N}\sum _{{\mathbf{k}}^{\prime \prime }}\frac{1}{\beta }\sum _{{\omega }_l} \\ \Lambda (\mathbf{k},i{\omega }_n;{\mathbf{k}}^{\prime \prime },i{\omega }_l)\Delta G({\mathbf{k}}^{\prime \prime },i{\omega }_l), \end{gathered}$$

$$\begin{gathered} {\Sigma }_{\sigma }(\mathbf{k},i{\omega }_n)=\frac{U}{2}\left(n-\sigma m\right) \\ -\frac{U}{N^2}\sum _{{\mathbf{k}}^{\prime \prime },\mathbf{q}}\frac{1}{{\beta }^2}\sum _{{\omega }_l,{\nu }_m}{G}_{\sigma }({\mathbf{k}}^{\prime \prime },i{\omega }_l)G_{\bar{\sigma }}({\mathbf{k}}^{\prime \prime }+\mathbf{q},i{\omega }_l+i{\nu }_m) \\ \times {\Gamma }_{\sigma \bar{\sigma }}({\mathbf{k}}^{\prime \prime },i{\omega }_l,\mathbf{k},i{\omega }_n;\mathbf{q},i{\nu }_m)G_{\bar{\sigma }}(\mathbf{k}+\mathbf{q},i{\omega }_n+i{\nu }_m), \end{gathered}$$

$$n$$

$$m$$

$$\Sigma (\mathbf{k},{\omega }_{+})=\frac{1}{2}\left[{\Sigma }_{\uparrow }(\mathbf{k},i{\omega }_n)+{\Sigma }_{\downarrow }(\mathbf{k},i{\omega }_n)\right].$$

$$\begin{gathered} G_{\sigma }(\mathbf{k},\omega ) \\ =\frac{1}{\omega +\mu -ϵ(\mathbf{k})+\sigma \left[h-\Delta \Sigma (\mathbf{k},\omega )\right]-\Sigma (\mathbf{k},\omega )}, \end{gathered}$$

$${\Lambda }_{\sigma \bar{\sigma }}(i{\omega }_n,i{\omega }_l)=\frac{U}{1+{\mathcal{S}}_{\sigma \bar{\sigma }}(i{\omega }_l,i{\omega }_n)},$$

$$\begin{gathered} {\mathcal{S}}_{\sigma \bar{\sigma }}(i{\omega }_l,i{\omega }_n)=\frac{1}{\beta }\sum _{{\nu }_m}{K}_{\sigma \bar{\sigma }}(i{\omega }_l,i{\omega }_n;-i{\nu }_m) \\ \times G_{\sigma }(i{\omega }_{n+m})G_{\bar{\sigma }}(i{\omega }_{l-m}). \end{gathered}$$

$$\begin{gathered} \frac{1}{\beta }\sum _{{\omega }_l}[\beta {\delta }_{n,l}+{\Lambda }_{\sigma \bar{\sigma }}(i{\omega }_n,i{\nu }_m+i{\omega }_l)G_{\bar{\sigma }}(i{\nu }_m+i{\omega }_l) \\ {G}_{\sigma }(i{\omega }_l)]K_{\sigma \bar{\sigma }}(i{\omega }_l,i{\omega }_{n^{\prime }};i{\nu }_m)=-\frac{1}{\beta }\sum _{{\omega }_l}{\Lambda }_{\sigma \bar{\sigma }}(i{\omega }_n,i{\omega }_l) \\ \times G_{\bar{\sigma }}(i{\omega }_l)G_{\sigma }(i{\omega }_l){\Lambda }_{\sigma \bar{\sigma }}(i{\omega }_l,,i{\omega }_{n^{\prime }}). \end{gathered}$$

$${⟨X(i{\omega }_l)⟩}_l=\frac{1}{\beta }\sum _{{\omega }_l>0}X(i{\omega }_l).$$

$$\begin{gathered} \sum _{\alpha =\pm 1}{⟨X(i{\omega }_n,i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l})Y(i{\omega }_{\alpha l},i{\omega }_{n^{\prime }})⟩}_l \\ \to \sum _{\alpha =\pm 1}\{X(i{\omega }_n,0_{\alpha }){⟨G(i{\omega }_{\alpha l})G(i{\omega }_{\alpha l-m})Y(i{\omega }_{\alpha l-m},i{\omega }_{n^{\prime }})⟩}_l \\ +{⟨X(i{\omega }_n,i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l}⟩}_l)Y(0_{\alpha },i{\omega }_{n^{\prime }}) \\ -X(i{\omega }_n,i{0}_{\alpha }){⟨G(i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l})⟩}_{l}Y(0_{\alpha },i{\omega }_{n^{\prime }})\}, \end{gathered}$$

$$\begin{gathered} {⟨X(i{\omega }_n,i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l})Y(i{\omega }_{\alpha l},i{\omega }_{n^{\prime }})⟩}_l \\ \to X(i{\omega }_n,i{0}_{\alpha }){⟨G(i{\omega }_{\alpha l+m})G(i{\omega }_{\alpha l})⟩}_{l}Y(0_{\alpha },i{\omega }_{n^{\prime }}). \end{gathered}$$

$${\Lambda }_{\uparrow \downarrow }\equiv {\Lambda }_{\uparrow \downarrow }(0_{+},0_{-})=\frac{U}{1+{\mathcal{K}}_{\uparrow \downarrow }{X}_{\uparrow \downarrow }}$$

$$\begin{gathered} X_{\uparrow \downarrow }={\int }_{-\infty }^{\infty }\frac{dx}{\pi }\{\frac{\Re \left[G_{\uparrow }(x_{+})G_{\downarrow }(-x_{+})\right]}{\sinh (\beta x)} \\ \times \Im \left[\frac{1}{D_{\uparrow \downarrow }(-x_{+})}\right]-f(x)\Im \left[\frac{G_{\uparrow }(x_{+})G_{\downarrow }(-x_{+})}{D_{\uparrow \downarrow }(-x_{+})}\right]\}. \end{gathered}$$

$$\begin{gathered} D_{\uparrow \downarrow }({\Omega }_{+})=1+{\Lambda }_{\uparrow \downarrow }[{⟨G_{\downarrow }(x+{\Omega }_{+})\Im G_{\uparrow }(x_{+})⟩}_x \\ +{⟨G_{\uparrow }(x-{\Omega }_{+})\Im G_{\downarrow }(x_{+})⟩}_x]. \end{gathered}$$

$$\begin{gathered} {⟨G_s(x+\omega )G_{s^{\prime }}(x+{\omega }^{\prime })⟩}_x=-{\int }_{-\infty }^{\infty }\frac{dx}{\pi }f(x) \\ \times G_s(x+\omega )G_{s^{\prime }}(x+{\omega }^{\prime }). \end{gathered}$$

$${\mathcal{K}}_{\uparrow \downarrow }=-{\Lambda }_{\uparrow \downarrow }^2{⟨\Im \left[G_{\downarrow }(x_{+})G_{\uparrow }(x_{+})\right]⟩}_x$$

$$a=D_{\uparrow \downarrow }(0)=1+{\Lambda }_{\uparrow \downarrow }{g}_{\uparrow \downarrow }(+).$$

$$1>1-\frac{U-{\Lambda }_{\uparrow \downarrow }}{X_{\uparrow \downarrow }{\Lambda }_{\uparrow \downarrow }^2}=D_{\uparrow \downarrow }(0)>0.$$

$$\Delta \Sigma =\frac{\Lambda }{2}\left[{⟨\Im G_{\downarrow }(x_{+})⟩}_x-{⟨\Im G_{\uparrow }(x_{+})⟩}_x\right].$$

$$\begin{gathered} {\Sigma }_{\uparrow }({\omega }_{+})=\frac{U}{2}\left(n-\sigma m\right) \\ +U{\int }_{-\infty }^{\infty }\frac{dx}{\pi }\{f(x+\omega )\frac{{\varphi }_{\uparrow \downarrow }(x_{-})}{D_{\uparrow \downarrow }(x_{-})}\Im G_{\downarrow }(x_{+}+\omega ) \\ -b(x)G_{\downarrow }(x_{+}+\omega ){\Im }_x\left[\frac{{\varphi }_{\uparrow \downarrow }(x_{+})}{D_{\uparrow \downarrow }(x_{+})}\right]\}, \end{gathered}$$

Full text

I Introduction

Many-body systems are those in which particle interactions cannot be neglected. In particular, when they are strong, they cause critical fluctuations and may lead to qualitative changes and phase transitions in extended systems. The critical behavior must be treated self-consistently. Even simple models of correlated elementary objects are, however, unsolvable. We hence must resort to approximations, apart from a few limiting cases where exact partial solutions exist. The approximations are either numerical or semi-analytic. The numerical approach tends to offer unbiased approximations with all degrees of freedom left in play, while the analytic approximations are based on a reduction of complexity of interaction effects. The numerical calculations offer good quantitative predictions that can set trends in the dependence of the solution on the model input parameters. The latter schemes aspire to reproduce qualitative features of the exact solution. They, unlike the numerical methods, can address and control singularities directly.

The first successful effort to tame the non-analytic behavior at the critical point of the continuous phase transitions in classical statistical models is the Landau mean-field theory in which the free energy is expanded in the small order parameter around the critical point Landau (1937). This simplest local and static self-consistent approximation inspired further attempts to improve upon it by including dynamical corrections Brout (1959, 1960); Horwitz and Callen (1961). But these attempts failed to match consistently ordered and disordered phases at the critical point Englert (1963). These inconsistencies were removed later with a quite different improvement of the mean field by using scaling arguments and the renormalization group Wilson (1971a, b). A mean-field approximation remains nevertheless the starting point also for the renormalization-group construction in analytic treatments.

The concept of a local, comprehensive mean-field theory was revitalized by realizing that it can be obtained as an exact solution in the limit of high spatial dimensions for classical spin models Brout (1960); Fisher and Gaunt (1964); Thompson (1974); Sherrington and Kirkpatrick (1975), as well as for quantum itinerant models Metzner and Vollhardt (1989); Janiš (1989); Brandt and Mielsch (1989); Janiš (1991). In particular, the limit to high spatial dimensions in models of correlated and disordered electrons initiated a boom in the applications of the dynamical mean-field theory (DMFT) Georges et al. (1996); Kotliar et al. (2006). The major asset of the DMFT is its unbiased inclusion of quantum fluctuations missed in classical static, weak-coupling theories. It hence offers a reliable way to investigate the strong-coupling limit of correlated electron systems at low temperatures. Since the single-impurity Anderson model (SIAM) is contained as the first, non-self-consistent iteration to the DMFT, the advanced methods of SIAM have been used to derive impurity solvers for the the dynamical mean-field approximations.

The standard mean-field theories, including DMFT, introduce renormalizations only for one-particle quantities represented by order parameters or the self-energy. They contain no renormalizations of vertex functions and hence, there is no direct control of the singularities in the Bethe-Salpeter equations for the two-particle response functions. The attempts to include two-particle and vertex renormalizations in the perturbation expansion are presently mostly made via non-local corrections to the DMFT Rohringer et al. (2018) or loop corrections with the functional renormalization group approach Kopietz et al. (2010). Recently, a two-particle self-consistency resulting from the parquet approach was used to construct an analytic mean-field theory for strong coupling Janiš and Augustinský (2007, 2008); Janiš et al. (2017a, b). The solution was shown to reproduce qualitatively correctly the Kondo effect in the strong-coupling limit of SIAM at zero temperature. It can be viewed upon as a consistent generalization of the Hartree approximation to strong coupling which is free of the spurious transition to an ordered phase of the weak-coupling solution. An effective interaction, as the only two-particle object determined self-consistently, is a static approximation to the irreducible vertex in the electron-hole scattering channel which resembles the GW construction Hedin (1965); Aryasetiawan and Gunnarsson (1997).

Although popularity of two-particle approaches has increased in recent years Janiš (1998, 1999); Toschi et al. (2007); Rubtsov et al. (2008); Rohringer et al. (2011); Rubtsov et al. (2012); Rohringer et al. (2013); Valli et al. (2015); Hirschmeier et al. (2015); Ayral and Parcollet (2016); Kugler and Delft (2018); Del Re et al. (2019), complexity of two-particle vertices demands the application of heavy numerics to reach quantitative results unless further approximations are used Yang et al. (2009); Tam et al. (2013); Li et al. (2016). The numerical solutions do not allow for the identification of the relevant degrees of freedom and for the control of the non-analytic behavior in the critical regions of vertex functions offered by mean-field theories. One has to introduce specific simplifications making the use of two-particle functions in the construction of mean-field approximations effective.

A few attempts were also made to simplify approximations for two-particle functions and to derive a set of static effective interactions determined self-consistently Vilk and Tremblay (1997) and Kusunose (2010). The former approach determines effective interactions in the charge and magnetic channels from local sum rules, while the latter determines the irreducible vertices in all two-particle channels from the crossing symmetry. Neither of these approaches was able to determine the Kondo scale in SIAM analytically.

The approximate scheme developed in Refs. Janiš and Augustinský (2007, 2008); Janiš et al. (2017a, b) used three leading principles of simplifying the equations for the two-particle vertices to end up with an analytic, mean-field-like theory of quantum criticality. 1) The full scheme of the parquet equations was replaced by a reduced set, 2) the dependence of the irreducible vertex on its dynamical variables was suppressed, and 3) the standard construction of approximate theories of Baym and Kadanoff Baym and Kadanoff (1961); Baym (1962) had to be left to comply with the Ward identity and to keep the approximation conserving whereby two self-energies were introduced. These steps were performed pragmatically, goal-directed without fully clarifying their general meaning for the construction of analytic approximations with the self-consistent determination of two-particle vertex functions.

The aim of this paper is to present a systematic derivation of mean-field theories with a two-particle self-consistency based on the parquet equations with the necessary simplifications leading to semi-analytic conserving approximations. The resulting approximations are free of unphysical behavior and spurious phase transitions and are applicable in strong coupling, both in disordered as well as in ordered phases of models of correlated electrons. In particular, we explain why a reduction of the full set of parquet equations with the bare interaction as the fully two-particle irreducible vertex is needed to reach the quantum critical behavior in the Bethe-Salpeter equations. We further explain how to separate the relevant from irrelevant dynamical fluctuations near the singularities in the two-particle vertex and how to decouple convolutions of fermionic Matsubara frequencies in the Bethe-Salpeter equations to make them analytically solvable. But most importantly, we give the proper meaning to two self-energies understood now as parts of a single self-energy with even and odd symmetry with respect to the symmetry-breaking field controlling the critical behavior. The former self-energy is determined from the dynamical Schwinger-Dyson equation. The latter is coupled with the electron-hole irreducible vertex via the Ward identity linearized in the symmetry-breaking field. It plays the role of the order parameter. The full self-energy is then compatible with both, the Ward identity and the Schwinger-Dyson equation. We thus set a framework for systematic improvements of the impurity solver from Refs. Janiš and Augustinský (2007, 2008); Janiš et al. (2017a, b) to dynamical mean-field-like approximations with a two-particle self-consistency offering a reliable description of thermodynamic and spectral properties in the strong-coupling limit in different settings of impurity and bulk models of correlated electrons.

II Two-particle self-consistency: Reduced parquet equations

A controlled and reliable way to suppress spurious transitions of the weak-coupling mean-field approximation is to introduce a two-particle self-consistency where two-particle vertices are determined self-consistently from non-linear equations. Only integrable singularities survive there as real phase transitions. One possibility to reach a two-particle self-consistency is to replace the bare interaction in response functions by effective ones determined self-consistently Vilk and Tremblay (1997). A more systematic way to do this is to use the parquet equations introduced to condensed matter by De Dominicis and Martin Dominicis and Martin (1964a, b). The full set of parquet equations was used to solve the soft X-ray problem Roulet et al. (1969), to understand the local-moment formation Weiner (1970, 1971) and was also applied to SIAM Bickers and White (1991); Bickers and Scalapino (1992). Unfortunately, the parquet equations simultaneously self-consistent at the one and two-particle level have not brought much progress beyond the fluctuation exchange, where only one-particle self-consistency is kept Bickers and Scalapino (1989). Although the spurious phase transition of the Hartree theory was suppressed with the one-particle self-consistency no Kondo limit was reproduced with it Hamann (1969); Chen and Bickers (1992). There is a general problem with the full set of parquet equations with the bare interaction as the completely irreducible vertex. Quantum criticality and the Kondo effect in SIAM are completely missed Janiš (2006). One must either replace the bare interaction by a more complex vertex or reduce appropriately the parquet equations to be able to reach the Kondo limit. Our aim is to find an analytically tractable mean-field theory for strongly correlated electron systems for which purpose we reduce the parquet equations and simplify their complexity.

II.1 Reduction scheme of the parquet equations

II.1.1 Two-channel parquet equations

The full parquet scheme contains three scattering channels represented by three different Bethe-Salpeter equations for the full two-particle vertex. They contain sums of multiple scatterings of singlet electron-hole ($eh$) pairs, electron-electron pairs ($ee$), and triplet electron-hole pairs ($\overline{eh}$) Bickers (1991). It is, however, sufficient to use only the singlet electron-electron and electron-hole multiple scatterings, since the third, $\overline{eh}$ channel, shares the critical fluctuations with either of the two fundamental channels Janiš (1998). We use the Hubbard model and SIAM to study the strong-coupling limit in correlated electron systems. The critical behavior in strong coupling is not by this specification qualitatively affected.

The Bethe-Salpeter equation for the full singlet two-particle vertex $\Gamma_{\sigma\bar{\sigma}}$ with the irreducible vertex $\Lambda^{eh}_{\sigma\bar{\sigma}}$ and $\bar{\sigma}=-\sigma$ in the electron-hole channel is

[TABLE]

where $N$ is the number of lattice sites, $\beta=1/k_{B}T$ and $G_{\sigma}$ is the propagator of the electron with spin $\sigma$.

Analogously the Bethe-Salpeter equation for the same vertex in the electron-electron channel is

[TABLE]

If we introduce a fully two-particle irreducible vertex $I_{\sigma\bar{\sigma}}$ we can use the parquet decomposition of the full vertex

[TABLE]

The parquet decomposition holds if the set of reducible diagrams in the electron-hole channel has no overlap with the set of reducible diagrams in the electron-electron channel Janiš (2009). The two-particle self-consistency is then obtained by replacing the full vertex $\Gamma_{\sigma\bar{\sigma}}$ by the above parquet decomposition in the Bethe-Salpeter equations (1) and (2) to obtain a set of self-consistent equations for the irreducible vertices $\Lambda^{eh}$ and $\Lambda^{ee}$. One standardly chooses the bare Hubbard interaction $U$ as the vertex irreducible in both two-particle scattering channels.

We kept notation $\Gamma$ for the full vertex and used symbol $\Lambda^{\alpha}$ for the irreducible vertex in the channel $\alpha$. Recent publications on the parquet equations as reviewed in Ref. Rohringer et al. (2018) used a slightly different notation. We related these two notations in Table 1 and compared them also with that in a seminal paper on the parquet equations, Ref. Bickers (1991) and with that in an early application of the parquet equations in condensed matter Roulet et al. (1969).

II.1.2 Critical region of the two-particle vertex

One of the Bethe-Salpeter equations approaches a singularity, a divergence in the full vertex at the critical point, when we increase the particle interaction. It is Eq. (1) for the magnetic systems (repulsive interaction) and Eq. (2) for the superconducting systems (attractive interaction). This can be seen from the weak-coupling approximation in the Bethe-Salpeter equation where the irreducible vertex $\Lambda$ is replaced by the bare interaction. We assume that this critical behavior can analytically be continued to strong coupling where the bare interaction must be renormalized and replaced by vertex $\Lambda$. It means that the divergence in the two-particle vertex emerges in the limit $\mathbf{q}\to q_{0}$ and $\nu_{m}\to 0$ and the irreducible vertex in the singular Bethe-Salpeter equation remains bounded and non-singular in the low-energy limit $\mathbf{q}\to\mathbf{q}_{0}$ and $\nu_{m}\to 0$ 111The high-frequency structure of the irreducible vertex may be rather complex with non-analyticities Janiš and Pokorný (2014); Schäfer et al. (2013); Chalupa et al. (2018). They are related to increasing imaginary part of the self-energy and are precursors of the metal-insulator transition. They may affect the low-frequency behavior only beyond the Fermi-liquid regime, since the imaginary part of the self-energy vanishes at the Fermi energy for Fermi liquids.. We will explicitly consider the magnetic case with the repulsive interaction and criticality in Bethe-Salpeter equation (1).

We can single out the relevant fluctuations in the critical region in Eq. (2) in the spirit of the renormalization group. The relevant fluctuations are those that make the dominant contribution to vertex $\Lambda_{\sigma\bar{\sigma}}^{ee}(\mathbf{k},i\omega_{n},\mathbf{k}^{\prime\prime},i\omega_{l};\mathbf{q}+\mathbf{k}^{\prime}-\mathbf{k}^{\prime\prime},i\nu_{m+n^{\prime}-l})$. They are controlled by the transfer momentum $\mathbf{q}+\mathbf{k}^{\prime}-\mathbf{k}^{\prime\prime}\to\mathbf{q}_{0}$ and frequency $i\nu_{m+n^{\prime}-l}\to 0$. For simplicity we assume a homogeneous order and choose $\mathbf{q}_{0}=\mathbf{0}$. The fluctuations in the fermionic variables remain irrelevant in the critical region and can be neglected. Equation (2) in the critical region then reduces to

[TABLE]

where we used the parquet equation (3) with the bare interaction as the completely irreducible vertex, $I_{\sigma\bar{\sigma}}=U$.

II.1.3 Reduction of parquet equations

The assumption about the possibility to follow continuously the critical behavior from weak to strong coupling demands that the irreducible vertex remains free of divergences in the low-energy limit. However, equation (4) may lead to divergences in $\Lambda_{\sigma\bar{\sigma}}^{eh}$ when $\mathbf{k}-\mathbf{k}^{\prime}=\mathbf{0}$ and $\omega_{n}-\omega_{n^{\prime}}=0$. If there is no correction to the bare interaction, the solution of this equation never reaches a critical behavior. The irreducible vertex would become also divergent which is incompatible with the two-particle self-consistency in the parquet equations. The critical behavior would be fully suppressed. To preserve the critical behavior one should extend the two-channel approximation by including the third channel and by replacing the bare interaction by a more complex completely irreducible vertex that would lead to the cancellation of the singular contributions to $\Lambda_{\sigma\bar{\sigma}}^{eh}$ Janiš (2006). This is, however, a tremendous task that would prevent reaching the desired objective of the analytic control of the critical behavior.

Alternatively, one can resort to a “poor-man approach” and keep the critical behavior of the non-renormalized theory by removing the super-divergent term from the convolution of two divergent vertices in Eq. (4). This leads then to a reduction of the parquet equations. This is achieved by the following replacement of the convolution on the right-hand side of Eq. (4)

[TABLE]

where we denoted $K_{\sigma\bar{\sigma}}=\Gamma_{\sigma\bar{\sigma}}-\Lambda^{eh}_{\sigma\bar{\sigma}}$, the reducible vertex in the electron-hole channel with $\Lambda=\Lambda^{eh}$ and used $\mathbf{k}^{\prime\prime}=\mathbf{k}^{\prime}+\mathbf{q}+\mathbf{Q}$, $\mathbf{q}^{\prime}=\mathbf{q}+\mathbf{k}+\mathbf{k}^{\prime}$, $m^{\prime}=m+n+n^{\prime}$. Since the completely irreducible vertex is a static constant, the irreducible vertex $\Lambda_{\sigma\bar{\sigma}}^{eh}(\mathbf{k},i\omega_{n},\mathbf{k}^{\prime},i\omega_{n^{\prime}};\mathbf{q},i\nu_{m})$ depends only on fermionic variables $\mathbf{k}^{\prime}+\mathbf{q}$, $i\omega_{n^{\prime}+m}$ and $\mathbf{k}$, $i\omega_{n}$. We redefine the transfer momentum in the irreducible vertex $\Lambda_{\sigma\bar{\sigma}}(\mathbf{k},\omega_{n};\mathbf{k}^{\prime},i\omega_{n^{\prime}};\mathbf{q},i\nu_{m})\to\Lambda_{\sigma\bar{\sigma}}(\mathbf{k},i\omega_{n},\mathbf{q}+\mathbf{k}^{\prime},i\omega_{m+n^{\prime}})$ to simplify the dependence of the dynamical variables.

The irreducible vertex is then determined from an integral equation

[TABLE]

Its diagrammatic representation is plotted in Fig. 1.

There is no change in the critical Bethe-Salpeter equation (1). The reducible vertex then is

[TABLE]

with its diagrammatic representation in Fig. 2.

Equations (6) and (7) form the set of the reduced parquet equations that introduce a two-particle self-consistency allowing for the extension of the critical behavior from weak-coupling approximations continuously to strong coupling. The reduced parquet equations (6)-(7) are the starting point for the investigation of the critical behavior of the two-particle vertex.

II.1.4 Low momentum and frequency limit

The reduced parquet equations (6)-(7) in their general form do not allow for analytic continuation to real frequencies with contour integrals involving Fermi and Bose distributions. To reach the goal of analytic control, one has to resort to specific limits separating the fermionic and bosonic degrees of freedom. We will consider two cases of the low-momentum and low-frequency limit of external variables of the two-particle vertices, $\mathbf{q},\mathbf{k},\mathbf{k}^{\prime}\to\mathbf{0}$ and $\nu_{m,}\omega_{n},\omega_{n^{\prime}}\to 0$, where we can obtain explicit solutions of the reduced parquet equations. The two cases differ in the ratio of the bosonic and fermionic variables. If we assume $|\mathbf{k}|/|\mathbf{q}|,|\mathbf{k}^{\prime}|/|\mathbf{q}|\to 0$ and simultaneously $\omega_{n}/\nu_{m},\omega_{n^{\prime}}/\nu_{m}\to 0$ we reproduce the static approximation introduced and analyzed in Refs. Janiš et al. (2017a, b). This approximation works well only at zero temperature and in the spin-symmetric state. To cover also non-zero temperatures and the ordered state we have to consider the opposite ratio, $|\mathbf{q}|/|\mathbf{k}|,|\mathbf{q}|/|\mathbf{k}^{\prime}|\to 0$ with $k,k^{\prime}\approx k_{F}$ and $\nu_{m}/\omega_{n},\nu_{m}/\omega_{n^{\prime}}\to 0$. We will now turn the the latter case.

Since the irreducible vertex $\Lambda_{\uparrow\downarrow}(\mathbf{k}^{\prime}+\mathbf{Q},i\omega_{n^{\prime}+l};\mathbf{k}-\mathbf{Q},i\omega_{n-l})$ in Eq. (6)is bounded in the critical region and only small values of momentum $\mathbf{Q}$ and frequency $\nu_{l}$ are relevant, we can neglect its fluctuations. We then obtain an explicit equation for the irreducible verrtex

[TABLE]

The equation for the reducible vertex for small transfer momentum and frequency in the critical region reads as

[TABLE]

We neglected the dependence of the sum on the right-hand side of Eq. (9) on the transfer momentum $\mathbf{q}$ and frequency $\nu_{m}$. It is irrelevant for the critical behavior of the reducible vertex. Equation (9) cannot, however, be solved analytically and further approximations are needed. We do it in the following section III.1.

II.2 Even and odd parts of the self-energy

The reduced parquet equations determine the irreducible and reducible vertices in the scattering channel with a singularity in the respective Bethe-Salpeter equation. The full two-particle vertex is then given as

[TABLE]

We use the convention with $k,k^{\prime}$ being the incoming and outgoing energy-momentum of the electron and $q$ is the transfer energy-momentum between the electron and the hole for the repulsive interaction studied here.

The one-particle propagators in the parquet equations are treated as input. In conserving theories the one- and two-particle Green functions are, however, related. The problem of correlated electron systems is that there are two ways to match consistently the one-particle self-energy with the two-particle vertex. One way is the Ward identity and the other is the dynamical Schwinger-Dyson equation Baym (1962). We recently demonstrated that no approximate solution can obey both relations exactly Janiš et al. (2017a). Neither of the two relations may, however, be disregarded. The former is needed for thermodynamic consistency and making the approximation conserving. The latter comes from the microscopic quantum dynamics. We cannot guarantee both relations with a single vertex and a single self-energy. We must either use a single self-energy and two two-particle vertices or vice versa. Ambiguity in the two-particle vertices leads to ambiguous criticality, thermodynamic inconsistencies and inability to continue the approximations beyond the critical point in the Bethe-Salpeter equation to the ordered phase. We are hence forced to use just a single two-particle vertex and introduce two self-energies to keep the approximate theories free of inconsistencies.

The two self-energies in approximate treatments can be introduced as even and odd parts of a single self-energy. The symmetry is set by the field controlling the critical fluctuations, conjugate to the order parameter. It is the longitudinal magnetic field in the present construction. The odd, or anomalous part of the self-energy will be determined via the Ward identity from the normal part, having the even symmetry, of the irreducible vertex $\Lambda_{\uparrow\downarrow}$ of the singular Bethe-Salpeter equation and the odd/anomalous part of the one-electron propagator. The Ward identity is, however, a functional differential equation that cannot be resolved for the self-energy exactly from the given two-particle vertex. To reach a qualitative thermodynamic consistency, it is sufficient to solve the Ward identity only linearly with respect to the symmetry-breaking field. The odd part or the thermodynamic self-energy is then defined as

[TABLE]

where the normal part of the irreducible vertex with even symmetry with respect to the magnetic field is $\Lambda(\mathbf{k},i\omega_{n},\mathbf{k}^{\prime\prime},i\omega_{l})=\left[\Lambda_{\uparrow\downarrow}(\mathbf{k},i\omega_{n},\mathbf{k}^{\prime\prime},i\omega_{l})+\Lambda_{\downarrow\uparrow}(\mathbf{k},i\omega_{n},\mathbf{k}^{\prime\prime},i\omega_{l})\right]/2$. The odd part of the one-electron propagator is $\Delta G(\mathbf{k}^{\prime\prime},i\omega_{l})=\left[G_{\uparrow}(\mathbf{k}^{\prime\prime},i\omega_{l})-G_{\downarrow}(\mathbf{k}^{\prime\prime},i\omega_{l})\right]/2$. A diagrammatic representation of the Ward identity is presented in Fig. 3. It is evident that the anomalous self-energy is connected with the order parameter and vanishes in the spin-symmetric (paramagnetic) state, unlike the thermodynamic self-energy used in Refs. Janiš et al. (2017a, b)

The normal part of the full self-energy has even symmetry with respect to the magnetic field and will be determined from the Schwinger-Dyson equation for the spin-dependent self-energy

[TABLE]

where the total particle density $n$ and magnetization $m$ read

[TABLE]

The full two-particle vertex is defined in Eq. (10) and the Schwinger-Dyson equation is diagrammatically represented in Fig. 4.

This spin-polarized self-energy contains also odd powers of the magnetic field. The odd part of the self-energy is already determined by the Ward identity and the Schwinger-Dyson equation should relate only the normal part of the self-energy with the two-particle vertex. The spin-dependent self-energy from Eq. (12) must then be symmetrized to acquire the desired even symmetry with respect to the magnetic field. The normal part of the dynamical (spectral) self-energy is

[TABLE]

The full self energy in the one-particle propagator is a sum of the anomalous self-energy, Eq. (11), and the spectral self-energy, Eq. (15). The full one-electron propagator is

[TABLE]

where $\mu$ is the chemical potential, and $\epsilon(\mathbf{k})$ is the dispersion relation.

Equations (11)-(16) define a theory with the full one-particle self-consistency. The only difference from the Baym-Kadanoff approach is the decomposition of the self-energy into its normal and anomalous parts, with even and odd symmetry with respect to the magnetic (symmetry-breaking) field determined from different exact equations. Should the Ward identity be compatible with the Schwinger-Dyson equation, we would recover the solution derived within the standard Baym and Kadanoff approach.

It appears, as we showed in earlier publications and will demonstrate later also in this paper, that the full one-particle self-consistency does not necessarily lead to approximations with the best and most reliable results. It may be convenient to relax the one-particle self-consistency and to replace the normal part of the self-energy in the propagators used in the Bethe-Salpeter and Schwinger-Dyson equation by a simpler one, $\Sigma_{0}(\mathbf{k},\omega)$. It can be selected to optimize the approximate solution. If done so, we call the approximate one-electron Green function a thermodynamic propagator and denote it $G^{T}$. The same propagator is then used also in the equation determining the anomalous self-energy in the Ward identity, Eq. (11). We showed earlier that the simplest approximation best suited for the Kondo asymptotics in SIAM is the Hartree approximation with $\Sigma_{0}(\mathbf{k},\omega)=Un/2$ Janiš et al. (2017a).

III Mean-field approximation for the reducible vertex function

III.1 Mean-field-like decoupling of frequency convolutions

The quantum character of many-body phenomena is manifested via the frequency dependence of the fundamental functions. Pure quantum critical behavior, free of spatial fluctuations, can be observed in the strong-coupling limit of impurity models with correlated electrons. We now resort to a mean-field treatment and suppress the spatial fluctuations in the vertex functions and resort to a local theory. We now apply the general theory to the Kondo behavior of the single-impurity Anderson model to produce an impurity solver for systems with strongly correlated electrons. We keep the full frequency dependence.

The irreducible vertex from Eq. (8) in the local approximation is

[TABLE]

We recall the notation $\bar{\sigma}=-\sigma$. The equation for the reducible vertex reduces to an integral (matrix) equation in Matsubara frequencies

[TABLE]

Although the fermionic variables at low temperatures are relevant only near the Fermi energy, $\omega\approx 0$, the integrals (sums) over the fermionic variables must be appropriately taken into account. We split the sum over Matsubara frequencies in sums over positive and negative frequencies. We use the following notation

[TABLE]

We cannot solve explicitly integral equation (18) but we can resort to a decoupling of the frequency convolutions in the spirit of a mean-value theorem for integrals of products of two positive functions. Fluctuations in the fermionic Matsubara frequencies are not relevant in the critical region of the singularity of the two-particle vertex. The fermionic frequencies are nevertheless important for keeping the approximation reliable also away from the critical region, in particular at non-zero temperatures. We leave only one of the two-particle vertices $X$ and $Y$ in the convolution of type $\langle XGGY\rangle$ dynamic, frequency dependent, while the frequency dependence of the other vertex will be replaced by a single mean value in the decoupling. It is reasonable to assume that the relevant values of the fermionic frequencies at low temperatures are only those from the vicinity of the Fermi energy. We hence choose the mean value to be the limit to the Fermi energy from above for the sum over positive frequencies and from below for the sum over negative ones. The decoupling scheme of the convolutions of fermionic Matsubara frequencies in the critical region of the low-temperature divergency of the reducible vertex then is

[TABLE]

where we denoted $0_{\alpha}=i\eta\ \mathrm{sign}(\alpha)$ and $\eta>0$ is infinitesimally small. This is a natural mean-field-like decoupling neglecting quadratic fluctuations beyond the mean value in averaging of products of operators. The decoupling of convolutions in Eq. (20) holds for functions analytic in the upper and lower complex planes. That is what we assume about the Green and vertex functions in Eq. (18). It is confirmed by explicit analytic continuation of sums in Eq. (18) to contour integrals resulting in an analytic expression for the reducible vertices $\Lambda_{\sigma\bar{\sigma}}(z,z^{\prime})$ and $K_{\sigma\bar{\sigma}}(z,z^{\prime};\zeta)$ for $z,z^{\prime},\zeta\in\mathbb{C}$.

We can solve Eq. (18) explicitly with this decoupling of the frequency convolutions. This approximate solution generates a rich analytic structure of the vertex functions with cuts along the real axis of the fermionic frequencies as well as along positive and negative diagonals in the plane of complex frequencies. The calculations are, however, rather lengthy and we leave the application of the dynamical decoupling from Eq. (20) to a separate publication. Here we resort to a simpler static decoupling of convolutions of fermionic frequencies

[TABLE]

The static decoupling is sufficient to gain a first qualitative picture of the strong-coupling regime. It offers a generalization of the static approximation from Ref. Janiš et al. (2017a) to the spin-polarized state and to non-zero temperatures.

III.2 Static mean-field approximation

We apply the general theory with the static decoupling of convolutions of fermionic Matsubara frequencies, Eq. (21), to the strong-coupling limit of SIAM where the Kondo effect is observed at half filling and zero temperature. The static decoupling with the relevant values of the irreducible vertex only from the Fermi energy simplifies the solution of the reduced parquet equations significantly without missing the strong-coupling limit. We then have $\Lambda_{\uparrow\downarrow}(\omega_{\sigma},\omega^{\prime}_{\tau})\to\Lambda_{\uparrow\downarrow}(0_{\sigma},0_{\tau})$. Consequently we reduce the reducible vertex to $K_{\uparrow\downarrow}(\omega^{\prime}+i\tau\eta^{\prime},\omega+i\sigma\eta;\Omega_{\rho})\to\mathcal{K}_{\uparrow\downarrow}(i\tau 0^{+},i\sigma 0^{+};\Omega_{\rho})=K_{\uparrow\downarrow}(\tau,\sigma)/D_{\uparrow\downarrow}(\Omega_{\rho})$, where $D_{\uparrow\downarrow}(\Omega_{\rho})$ is the determinant of the matrix of the kernel of the equation for the reducible vertex $\mathcal{K}_{\uparrow\downarrow}(i\tau 0^{+},i\sigma 0^{+};\Omega_{\rho})$ with $\Omega\to 0$.

One can solve the equations to determine $\Lambda_{\uparrow\downarrow}(0_{\sigma},0_{\tau})$ and $K_{\uparrow\downarrow}(\tau,\sigma)$. The reduction of the dependence of the vertex functions on the fermionic frequencies to their values at the Fermi energy is justified and works reliably in the quantum critical region of the singularity of the two-particle vertex, that is at zero temperature. Both vertex functions of fermionic frequencies are continuous across the Fermi energy at zero temperature, that is $\Lambda_{\uparrow\downarrow}(\sigma,\tau)=\Lambda_{\uparrow\downarrow}$ and $K_{\uparrow\downarrow}(\tau,\sigma)=K_{\uparrow\downarrow}$. The two values $\Lambda_{\uparrow\downarrow}(+,+)$ and $\Lambda_{\uparrow\downarrow}(+,-)$ split at non-zero temperatures where the values at the Fermi energy play gradually less and less dominant role. If we want to extrapolate consistently the low-temperature approximation also to high temperatures far above the Kondo temperature, one has to return to the full dynamical decoupling of convolutions of fermionic frequencies, Eq. (20) and keep the full dynamics of the irreducible vertex $\Lambda_{\uparrow\downarrow}(\omega_{\sigma},\omega^{\prime}_{\tau})$. A simpler option without leaving the static approximation is to select a single value of the irreducible vertex at the Fermi energy that screens the bare interaction efficiently at all temperatures and is dominating at high temperatures. It is $\Lambda_{\uparrow\downarrow}(0_{+},0_{-})$ that does the job and will be used to extend the zero-temperature static solution continuously to high temperatures Janiš and Klíč (2019).

The single irreducible vertex reduces then to an effective interaction defined from an equation

[TABLE]

where

[TABLE]

We denoted $x_{+}=x+i0^{+}$ and used an equality $f(x)+b(x)=1/\sinh(\beta x)$. We straightforwardly analytically continued the sums over Matsubara frequencies to spectral integrals with the Fermi, $f(x)=1/(e^{\beta x}+1)$ and the Bose, $b(x)=1/(e^{\beta x}-1)$, distributions.

The reduced parquet equations now have an explicit algebraic solution. The full, frequency dependent determinant $D_{\uparrow\downarrow}(\Omega_{+})$ reads as

[TABLE]

where we denoted

[TABLE]

The equation for the reducible vertex with frequencies near the Fermi energy $\mathcal{K}_{\uparrow\downarrow}=\mathcal{K}_{\uparrow\downarrow}(0_{-},0_{+})$ reads as

[TABLE]

We denote $g_{\uparrow\downarrow}(+)=\left\langle\Im\left[G_{\downarrow}(x_{+})G_{\uparrow}(x_{+})\right]\right\rangle_{x}$ and introduce a dimensionless Kondo scale as the zero value of function $D_{\uparrow\downarrow}(\Omega_{+})$ from Eq. (24),

[TABLE]

It measures the distance to the critical point $a=0$ and will affect the low-energy behavior in the strong-coupling limit.

We rewrite Eq. (26) by using Eq. (22) to another, more suitable form

[TABLE]

The right inequality guarantees stability of the solution and integrability and positivity of the integral $X_{\uparrow\downarrow}$ in the strong-coupling regime.

Equations (22)-(26) form a closed set of equations determining self-consistently the values of the vertex functions $\Lambda_{\uparrow\downarrow}$ and $\mathcal{K}_{\uparrow\downarrow}$ at the Fermi energy. They can be solved numerically in a straightforward way via iterations. The bare interaction and the one-electron propagators are input to these equations. The latter contain the normal and anomalous parts of the self-energy. The anomalous self-energy renormalizing the effect of the magnetic field in the static approximation is

[TABLE]

We recall that $\Lambda=\left(\Lambda_{\uparrow\downarrow}+\Lambda_{\downarrow\uparrow}\right)/2$ is the normal part of the irreducible vertex. But due to the electron-hole symmetry $\Lambda_{\uparrow\downarrow}=\Lambda_{\downarrow\uparrow}$.

The spin-polarized self-energy from the Schwinger-Dyson equation is

[TABLE]

where

[TABLE]

is the electron-hole bubble multiplied by the effective interaction.

The normal part of the dynamical self-energy then is $\Sigma(\omega_{+})=\left[\Sigma_{\uparrow}(\omega_{+})+\Sigma_{\downarrow}(\omega_{+})\right]/2$ and the one-electron propagator with the full self-energy for SIAM is

[TABLE]

where $\Delta$ is the width of the local level attached to conducting leads and is set as the energy unit.

The present construction allows for selecting the optimal degree of the one-electron self-consistency. It means that the propagators used in the perturbation expansion, in the Ward identity, the Bethe-Salpeter, and the Schwinger-Dyson equations, need not be the fully renormalized propagator from Eq. (32). We can use a thermodynamic propagator with an appropriately chosen normal self-energy $\Sigma_{0}(\omega)$ with which we control the degree of one-particle self-consistency

[TABLE]

where the Hartree term is calculated with the thermodynamic propagator

[TABLE]

It appears that to reach the best fit for the asymptotic form of the Kondo scale is $\Sigma_{0}(\omega)=0$, which we use in the thermodynamic propagator $G^{T}_{\sigma}(\omega)$ in the following section.

IV Strong-coupling limit at zero temperature of SIAM

The critical behavior of SIAM emerges in the asymptotic limit to strong-coupling of the spin- and charge-symmetric state at zero temperature. We now explicitly solve the equations within the static decoupling of the frequency convolutions at zero temperature where the irreducible vertex is continuous at the Fermi energy. We will analyze the spin polarized solution generally also away from half filling to demonstrate how the critical behavior fades away when moving from the spin and charge symmetric situation and leaving the strong-coupling regime.

IV.1 Analytic solution

The irreducible vertex, the effective interaction at zero temperature is

[TABLE]

with the determinant $D_{\uparrow\downarrow}(x_{+})$ and $K_{\uparrow\downarrow}$ defined in Eqs. (24) and (26).

We separately represent the imaginary and real parts of the dynamical self-energy. The imaginary part has the following representation

[TABLE]

Function $\phi_{\uparrow\downarrow}(x_{+})$ was defined in Eq. (31).

The thermodynamic propagator can be represented via an effective chemical potential and an effective magnetic field

[TABLE]

that are derived from the bare chemical potential and magnetic field together with the thermodynamic particle density and magnetization

[TABLE]

The unrenormalized thermodynamic propagator allows for analytic representations of thermodynamic quantities. The thermodynamic charge density and magnetization are

[TABLE]

The zero-field susceptibility from the thermodynamic propagator is

[TABLE]

The full physical spin-dependent propagator is

[TABLE]

The physical particle density and magnetization are

[TABLE]

respectively and they differ slightly from the thermodynamic ones. The zero-field physical susceptibility from the full propagator is

[TABLE]

Only small frequencies in $D_{\uparrow\downarrow}(\omega_{+})$ are relevant in the strong-coupling limit. We can then replace the full frequency dependence in $D_{\uparrow\downarrow}(\omega_{+})$ by a low-frequency expansion, keeping only the term linear in frequency. That is,

[TABLE]

The linear coefficient reads as

[TABLE]

where $\rho_{\sigma}=-\Im G_{\sigma}(0_{+})/\pi$ is the density of states at the Fermi energy.

The theory with the full one-particle self-consistency and the full propagator in all equations can be solved only numerically at intermediate and not too strong electron interaction. The theory with the thermodynamic propagator, on the other hand, allows for an explicit analytic representation also in the Kondo limit. It is defined when $a=1+\Lambda g_{\uparrow\downarrow}(0)\ll 1$. Using the thermodynamic propagator from Eq. (37) and the low-frequency asymptotics of the determinant from Eq. (47) we can explicitly evaluate the frequency integrals.

The asymptotic form of the effective interaction in strong coupling within the low-frequency asymptotics, Eq. (47) is

[TABLE]

where $G_{\sigma}=G_{\sigma}(0_{+})$ and

[TABLE]

[TABLE]

[TABLE]

The first equality holds always generally for arbitrarily renormalized one-electron propagator while the second one only for the thermodynamic propagator $G^{T}(\omega_{+})$ from Eq. (37). Since the solution is a local Fermi liquid where $\Sigma(0)=0$, the actual degree of the one-particle renormalization does not play a role in the strong-coupling asymptotics. Notice that the imaginary part of the linear coefficient of the determinant $D_{\uparrow\downarrow}(\omega_{+})$ is non-zero only in the spin-polarized state.

The asymptotic form of the imaginary part of the spin-polarized dynamical self-energy in strong coupling is

[TABLE]

where we denoted $\bar{a}=a/\Lambda$.

Finally, we can explicitly evaluate the Kondo scale as a function of the effective chemical potential and the effective magnetic field when we use the critical effective interaction $\Lambda=1/D_{0}$. We then obtain

[TABLE]

This explicit representation can be used to determine the boundary of the strong-coupling region. It is set to the point where the logarithm of the asymptotic Kondo scale goes through zero, that is, $UD_{0}=1$. The Kondo regime is then defined for the interaction strength obeying $UD_{0}\gg 1$.

IV.2 Numerical results

We used thermodynamic and full one-electron propagators to compare the results with the exact ones for SIAM at zero temperature. Since the equation for the irreducible vertex, Eq. (35), slightly differs from the vertex used in Refs. Janiš et al. (2017a, b) we also compare the two versions of the two-particle self-consistency.

The Kondo effect and the Kondo strong-coupling regime occur at half filling of the non-magnetic state. We compare in Fig. 5 the present version of the static approximation with the thermodynamic propagator from Eq. (37) with that of Ref. Janiš et al. (2017a) and with the spectral function from NRG. There is no much difference between the two versions of the static approximation. They both have the same enhancement of the satellite Hubbard bands with a narrower central Kondo-Suhl resonance compared to NRG.

The NRG calculations were performed with the NRG Ljubljana code Žitko (2014). A constant density of states of bandwidth $2D$ with $U/2D>100$ was used. Spectral functions were obtained from the DM-NRG algorithm of Ref. Hofstetter (2000). We opted for not correcting the spectral energies via the so-called self-energy trick. All results were recalculated from the typical NRG units of $D$ into the units of $\Delta$ as used in this paper.

One can improve upon the thermodynamic propagator and take into account the full one-particle self-consistency. We take the same full propagator in the parquet equations as well as in the Ward identity and the Schwinger-Dyson equation. We compared the two solutions for the spectral function of the non-magnetic state at half filling with the NRG result for weak and moderate interactions in Fig. 6. In weak coupling all three approaches give almost the same function. In intermediate coupling the fully self-consistent version delivers a better agreement with the NRG result for high frequencies and in positioning of the Hubbard satellite bands. The width of central quasiparticle peak is, however, missed in the strong-coupling regime of the self-consistent version, Fig. 7. The non-self-consistent version with the thermodynamic propagator from Eq. (37) predicts correctly the linear dependence of the logarithm of the Kondo scale on the interaction strength defined as the half width at half maximum (HWHM) of the central peak while the exponent of the self-consistent solution is one third. We plotted also another definition of the Kondo scale from factor $Z=(1-d\Sigma/d\omega|_{\omega=0})^{-1}$ which shows the same strong-coupling asymptotics. The shift of the curve calculated from our static approximation with the thermodynamic propagator is caused by a difference in the non-universal pre-factor at the logarithm of the Kondo scale. It is $\pi/8$ in the Bethe-ansatz solution with the Lorentzian density of states Hewson (1993) while it comes as $1/\pi$ from our parquet equations.

It is convenient, in particular in more complex models, to have as simple equations in the strong-coupling regime as possible. The asymptotic representation for the self-energy, Eq. (53), can do the job and it is capable to deliver the qualitatively correct three-peak spectral function, see Fig. 8. There is a good agreement in weak coupling with the full solution. The width of the central peak is asymptotically correct and only the satellite Hubbard bands are more pronounced compared to the full solution. The real and imaginary parts of the self-energy are plotted in Fig. 9. Since a one-particle self-consistency is missing, the real part of the self-energy at weak coupling does not reproduce correctly the negative slope at the Fermi energy.

When we move away from half filling the central peak slowly moves away from the Fermi energy as well as the lower Hubbard band (for occupation $n<1$) moves towards the central one and eventually merges with it, Fig. 10. It is the behavior discussed in more details in Ref. Janiš et al. (2017b).

A more interesting situation is when we break the spin-reflection symmetry, the case not discussed in our earlier publications. The magnetic field affects the Kondo strong-coupling behavior more significantly than the shift of the chemical potential from half filling. We show in Fig. 11 the full spectral function for several values of the magnetic field. We can see that already weak magnetic fields of order $h=0.2\Delta$ split the central peak into two separate ones and lower the height of the split peaks. The movement of the central peak in the magnetic field can be better demonstrated in the spectral function of the majority spin (up). We compared the spectral function for the majority spin with NRG in Fig. 12. Qualitative features of the spectral function are well reproduced in the static approximation with the thermodynamic propagator. A dip observed at weak and intermediate field is a consequence of insufficient one-particle self-consistency.

The conclusion that the magnetic field moves the solution from the strong-coupling regime at half filling faster than the normalized chemical potential $\mu-U/2$ can be demonstrated in the dependence of the negative logarithm of the Kondo scale at criticality, Eq. (54), on both variables, Fig. 13. The strong-coupling regime ends where the asymptotic result goes through zero. It happens for much smaller values of the magnetic field than for the normalized chemical potential. The difference between the initial value at $\mu-U/2=0$ of the static solution and the numerically exact one from NRG is due to the fact that $U=12\Delta$ is not yet the true Kondo regime, since the asymptotic limit $UD_{0}\gg 1$ in the asymptotic form in Eq. (54) has not yet been reached. Moreover, the exact value in the limit $U\to\infty$ is $\pi/8$ while the static solution gives $1/\pi$. It is interesting to notice that the bound for the strong-coupling regime from the static solution agrees well with the exact expression for the Kondo scale. There is naturally no boundary between strong and weak coupling in the full solution.

The extent to which quantum fluctuations are relevant can be measured by the strength of the renormalization of the bare interaction. That is, how much the irreducible vertex $\Lambda$, effective interaction, differs from the bare one. We plotted in Fig. 14 the dependence of $\Lambda_{\uparrow\downarrow}$ on the normalized chemical potential and magnetic field. We can again observe that the renormalization decreases faster with the increasing magnetic field than with the chemical potential. The effective interaction approaches the bare one in the weak-coupling limit with decaying quantum fluctuations. It is surprising how well the asymptotic form reproduces the full one in the spin-polarized state.

V Conclusions

The major desired asset of the mean-field theory is its relative simplicity that allows for the analytic control of the critical behavior. It is much easier to achieve this in classical many-body models than in quantum ones. The difference is made by the dynamics brought in by quantum fluctuations in strong coupling. Although the dynamical fluctuations can be studied within a local, dynamical mean-field theory, the only accessible numerical solution does not allow for the analytic control. Further approximations are needed. We presented in this paper a reduction scheme leading to a class of analytic mean-field theories of quantum fluctuations in strongly correlated electron systems. The resulting theories are thermodynamically consistent and conserving. They reconcile the thermodynamic Ward identity and the dynamical Schwinger-Dyson equation and are free of spurious transitions and unphysical behavior of the Hartree weak-coupling theory.

The central objects to be determined form the diagrammatic perturbation theory of the present construction are two-particle vertex functions instead of the one-particle self-energy. The reason for that is to achieve a two-particle self-consistency that guarantees that only integrable singularities and physical phase transitions can exist. The two-particle self-consistency is derived from the parquet equations for two-particle irreducible vertices. The analytic structure of the two-particle vertices is much more complex than that of the self-energy. A set of simplifications must be introduced to reach analytically controlled approximations.

There are three main steps in reaching a tractable analytic mean-field theory with a two-particle self-consistency. First, one has to reduce the full set of parquet equations in order to get rid of super-divergent terms generated by mixing of different scattering channels beyond the random-phase approximation. They are assumed to be canceled by higher-order contributions beyond the two-channel parquet equations with the bare interaction as the completely irreducible vertex. If not suppressed they would prevent reaching the quantum critical behavior. Second, one moves into the critical region of the Bethe-Salpeter equation where the relevant critical fluctuations of bosonic degrees of freedom are decoupled from the non-critical ones of the fermionic degrees of freedom. It is done in analogy with the renormalization group where only small frequencies controlling criticality are explicitly considered. This leads to an explicit form of the non-critical irreducible vertex from the singular Bethe-Salpeter equation. Third, convolutions of non-critical fermionic frequencies in the Bethe-Salpeter equation are decoupled in analogy with with the mean-value theorem. This may be done either fully dynamically or statically. Here we elaborated the static approximation leading to a Hartree-like approximation with a renormalized effective interaction.

The mean-field theory for vertex functions is not complete unless it determines the self-energy renormalizing the one-electron propagators used in the equations for the vertex functions. The consistent and conserving theory must properly relate the one and two-particle functions. We used the symmetry with respect to the external field conjugate to the order parameter of the new phase beyond the critical point of the Bethe-Salpeter equation. The normal part of the self-energy, with even symmetry with respect to the symmetry-breaking field, is determined from the Schwinger-Dyson equation using the two-particle vertex obtained from the parquet equations. The anomalous part, with odd symmetry with respect to the symmetry-breaking field, is determined from the Ward identity and the irreducible vertex from the parquet equations. The Ward identity and the Schwinger-Dyson equations connecting the self-energy and the two-particle vertex in different ways are reconciled in these approximations. Moreover, our approach allows for a flexible handling of the one-particle self-consistency, from non-self-consistent to fully self-consistent one-electron propagators to optimize the quantitative output.

The resulting mean-field theory with a two-particle self-consistency even in its simplest static version gives a qualitatively correct and thermodynamically consistent description of quantum criticality in SIAM. Although it is justified in the critical region of the singularities in the Bethe-Salpeter equations it can be extended to non-critical regions as well as to the ordered phase. It hence offers a qualitative picture of the transition from weak to strong coupling as well as from low to high temperatures in models of correlated electrons. The general construction can also be used beyond the local approximation to describe the low-temperature behavior of the low-dimensional models with no long-range order at non-zero temperatures.

Acknowledgments

VJ, PZ, and AK were supported by Grant No. 19-13525S of the Czech Science Foundation. VP was supported by grant INTER-COST LTC19045.