Strong-coupling formula for momentum-dependent susceptibilities in dynamical mean-field theory

TL;DR

This paper introduces a simplified, computationally efficient formula for calculating momentum-dependent susceptibilities in DMFT, applicable across different coupling regimes, enabling easier evaluation of spin and orbital susceptibilities.

Contribution

The authors derive a strong-coupling formula for susceptibilities in DMFT that reduces computational cost and is valid in both strong and weak coupling regimes.

Findings

The formula accurately reproduces intersite interactions like kinetic exchange and RKKY.

Numerical tests show wide applicability including weak-coupling regions.

Practical approximation levels facilitate easier susceptibility calculations.

Abstract

Computing momentum-dependent susceptibilities in the dynamical mean-field theory (DMFT) requires solving the Bethe-Salpeter equation, which demands large computational cost. Exploiting the strong-coupling feature of local fluctuations, we derive a simplified formula that can be solved at a considerably lower cost. The validity and the physical meaning of the formula are confirmed by deriving the effective intersite interactions in the strong-coupling limit, such as the kinetic exchange and RKKY interactions. Furthermore, numerical calculations for single-orbital and multiorbital models demonstrate surprisingly wider applicability including weak-coupling region. Based on this formula, we propose three levels of practical approximations that can be chosen depending on complexity of problems. Simpler evaluations of spin and orbital susceptibilities in multiorbital systems thus become possible within DMFT.

Tables

Table 1. Table 1: Classification of 16 operators σ ξ ​ τ η superscript 𝜎 𝜉 superscript 𝜏 𝜂 \sigma^{\xi}\tau^{\eta} in the two-orbital system. The column ‘irrep’ shows irreducible representation in cubic point group in Mulliken and Bethe notation. The index g and u (or + + and − - ) represent time-reversal even and odd, respectively.

label	operator(s)	degeneracy	irrep in Oh
charge	${\sigma }^0{\tau }^0$	1	A${}_{\text{1g}}$	${\Gamma }_1^{+}$
spin (AFM)	${\sigma }^x{\tau }^0$, ${\sigma }^y{\tau }^0$ ${\sigma }^z{\tau }^0$	3	T${}_{\text{1u}}$	${\Gamma }_4^{-}$
orbz (HS/LS)	${\sigma }^0{\tau }^z$	1	$\}$ E${}_{\text{g}}$	${\Gamma }_3^{+}$
orbx	${\sigma }^0{\tau }^x$	1
orb${}^{\text{3}}$	${\sigma }^x{\tau }^y$, ${\sigma }^y{\tau }^y$, ${\sigma }^z{\tau }^y$	3	T${}_{\text{2g}}$	${\Gamma }_5^{+}$
sp-orb1	${\sigma }^0{\tau }^y$	1	A${}_{\text{2u}}$	${\Gamma }_2^{-}$
sp-orb3x (EI)	${\sigma }^x{\tau }^x$, ${\sigma }^y{\tau }^x$, ${\sigma }^z{\tau }^x$	3	T${}_{\text{1u}}$	${\Gamma }_4^{-}$
sp-orb3z	${\sigma }^x{\tau }^z$, ${\sigma }^y{\tau }^z$, ${\sigma }^z{\tau }^z$	3	T${}_{\text{2u}}$	${\Gamma }_5^{-}$

Equations

$$O_{m\sigma ,m^{\prime }{\sigma }^{\prime }}(i)=c_{im\sigma }^{†}{c}_{i{m}^{\prime }{\sigma }^{\prime }},$$

$$[\mathbf{\chi }(\mathbf{q},i\Omega ){]}_{12,34}={\int }_0^{\beta }d\tau ⟨O_{12}(\mathbf{q},\tau )O_{43}(-\mathbf{q})⟩e^{i\Omega \tau },$$

$${\chi }_{12,34}(\mathbf{q},i\Omega )=T\sum _{\omega {\omega }^{\prime }}{X}_{12,34}(i\omega ,i{\omega }^{\prime };\mathbf{q},i\Omega ).$$

$$X_{12,34}(i\omega ,i{\omega }^{\prime };\mathbf{q},i\Omega )\equiv [\mathbf{X}(\mathbf{q},i\Omega ){]}_{(12,\omega ),(34,{\omega }^{\prime })}.$$

$${\mathbf{X}}^{-1}(\mathbf{q},i\Omega )$$

$$[X_0(\mathbf{q},i\Omega ){]}_{(12,\omega ),(34,{\omega }^{\prime })}$$

$$=-\frac{{\delta }_{\omega {\omega }^{\prime }}}{N}\sum _{\mathbf{k}}{G}_{31}(\mathbf{k},i\omega )G_{24}(\mathbf{k}+\mathbf{q},i\omega +i\Omega ).$$

$${\mathbf{X}}_{\text{loc}}^{-1}(i\Omega )$$

$${\mathbf{X}}^{-1}(\mathbf{q},i\Omega )$$

$$\mathcal{H}=\sum _{\mathbf{k}\sigma }{ϵ}_{\mathbf{k}}{c}_{\mathbf{k}\sigma }^{†}{c}_{\mathbf{k}\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$$({\chi }_{\mathbf{q}}^{\text{RPA}}{)}^{-1}=({\chi }_{0,\mathbf{q}}^{U=0}{)}^{-1}-U.$$

$$({\chi }_{\mathbf{q}}^{\text{RRPA}}{)}^{-1}={\chi }_{0,\mathbf{q}}^{-1}-{\chi }_{\text{0,loc}}^{-1}+{\chi }_{\text{loc}}^{-1}.$$

$$X_{\text{loc}}(i\omega ,i{\omega }^{\prime })\simeq \Phi (i\omega )\Phi (i{\omega }^{\prime }).$$

$$X_{\text{loc}}(i\omega ,i{\omega }^{\prime })$$

$$\Phi (i\omega )=\sqrt{s_0}{u}_0(i\omega )=\sqrt{s_0}{v}_0^{\ast }(i\omega ).$$

$$T{\chi }_{\text{loc}}\simeq {\overline{\Phi }}^2=s_0{\overline{u}}_0^2,$$

$${\mathbf{X}}_{\mathbf{q}}^{-1}={\mathbf{X}}_{\text{loc}}^{-1}-{\mathbf{\Lambda }}_{\mathbf{q}},$$

$${\mathbf{\Lambda }}_{\mathbf{q}}\equiv {\mathbf{X}}_{\text{0,loc}}^{-1}-{\mathbf{X}}_{0,\mathbf{q}}^{-1}.$$

$$\sum _{{\omega }^{\prime \prime }}{X}_{\text{loc}}(i\omega ,i{\omega }^{\prime \prime }){\Lambda }_{\mathbf{q}}(i{\omega }^{\prime \prime })X_{\text{loc}}(i{\omega }^{\prime \prime },i{\omega }^{\prime })\simeq {\lambda }_{\mathbf{q}}{X}_{\text{loc}}(i\omega ,i{\omega }^{\prime }),$$

$${\lambda }_{\mathbf{q}}\equiv \sum _{{\omega }^{\prime \prime }}{\Lambda }_{\mathbf{q}}(i{\omega }^{\prime \prime })\Phi (i{\omega }^{\prime \prime }{)}^2.$$

$${\chi }_{\mathbf{q}}^{\text{SCL}}=(1-{\lambda }_{\mathbf{q}}{)}^{-1}{\chi }_{\text{loc}}.$$

$$I_{\mathbf{q}}\equiv {\lambda }_{\mathbf{q}}/{\chi }_{\text{loc}},$$

$${\chi }_{\mathbf{q}}^{\text{SCL}}=({\chi }_{\text{loc}}^{-1}-I_{\mathbf{q}}{)}^{-1}.$$

$$I_{\mathbf{q}}\simeq T\sum _{\omega }{\Lambda }_{\mathbf{q}}(i\omega )\varphi (i\omega {)}^2,$$

$$u_0^{\text{atom}}(i\omega )=\frac{A}{2}\left(\frac{1}{i\omega +\mu }-\frac{1}{i\omega +\mu -U}\right),$$

$$A^2\simeq 8UTr(1-r),$$

$$s_0\simeq \frac{1}{4UTr(1-r)}.$$

$${\varphi }^{\text{atom}}(i\omega )=\frac{1}{i\omega +\mu }-\frac{1}{i\omega +\mu -U}.$$

$${\Lambda }_{\mathbf{q}}(i\omega )\simeq -2t^2{\gamma }_{\mathbf{q}}.$$

$$I_{\mathbf{q}}\simeq -\frac{4t^2{\gamma }_{\mathbf{q}}}{U},$$

$$\mathcal{H}$$

Full text

I Introduction

Strongly correlated materials including $d$ and/or $f$ electrons show rich quantum phases that attract academic interest and promise novel functionality. To find further interesting or useful properties, new compounds are searched for continuously. On the theoretical side, a qualitative or even quantitative prediction of phase transitions is desired to promote experimental researches for finding new materials.

Dynamical mean-field theory (DMFT) is a method which is extensively employed in investigations of strongly correlated systems Georges et al. (1996). The DMFT takes full account of local correlations, which is responsible for formation of a local moment, and therefore offers a good starting point in addressing long-range ordering. Applications of the DMFT have expanded from models to materials by combining first-principles calculations within the density functional theory (DFT+DMFT)Kotliar et al. (2006).

The DMFT is exploited not only for single-particle excitations but also for two-particle responses by means of the momentum-dependent susceptibilities $\chi(\bm{q},\omega)$. The dynamical part gives information of two-particle excitations including collective modes. On the other hand, the static component $\chi(\bm{q})$ allows systematic investigations of phase transitions. This method has been applied to various correlated models such as the Hubbard model Jarrell (1992); Lin et al. (2012); Kuneš and Augustinský (2014); Geffroy et al. (2019), the periodic Anderson model Jarrell (1995), and variants of the Kondo lattice model Otsuki et al. (2009a, b); Hoshino et al. (2010, 2013). It has also been applied to material calculations within the DFT+DMFT framework Matsumoto et al. (2009); Boehnke and Lechermann (2012); Yin et al. (2014).

An actual calculation of $\chi(\bm{q})$ requires much more effort compared to the ordinary self-consistency calculations: In addition to the single-particle Green function (or the self-energy), one needs to compute the two-particle Green function (or the vertex part) which depends on three frequencies. This is plugged into the Bethe-Salpeter (BS) equation to compute $\bm{q}$-dependent $\chi(\bm{q})$. This additional process of computing $\chi(\bm{q})$ becomes particularly harder in multiorbital systems, since the number of susceptibility components increases in proportion to $N_{\textrm{orb}}^{4}$, where $N_{\textrm{orb}}$ is the number of orbitals. Toward realistic susceptibility calculations including full orbitals, efforts have been made to develop an efficient treatment of the BS equation Tagliavini et al. (2018) and the vertex part Kuneš (2011); Rohringer et al. (2012); Kaufmann et al. (2017); Wentzell et al. ; Shinaoka et al. (2018); Krien . Related developments are addressed, for example, in solving the parquet equation Li et al. (2016) and diagrammatic extensions of the DMFT Rohringer et al. (2018); Galler et al. (2018).

In this paper, we establish a practical approximation that enables quick computations of $\chi(\bm{q})$, rather than pursuing the rigorous solution of the original BS equation. Using a strong-coupling feature of the vertex part, we first reduce the BS equation into a simple form which is similar to the random phase approximation (RPA). The difference with the RPA, however, is that the starting point of our formula is the local susceptibility, which fully takes local correlations into account. The effective interaction becomes non-local, describing, for example, the kinetic exchange interaction in the Hubbard model and the RKKY interaction in the periodic Anderson model. Thus, our simple formula captures correct physics of well-defined local moments in the strong-coupling limit (SCL). For this reason, we refer to this formula as SCL formula. It will be shown by numerical calculations that the SCL formula provides accurate results over wide range from weak- to strong-coupling region.

This paper is organized as follows. In Sec. II, we first review necessary equations for computing $\chi(\bm{q})$, and present an explicit example that shows the issue to be addressed. In Sec. III, we derive the SCL formula of $\chi(\bm{q})$ by reducing the BS equation. In Sec. IV, we prove that the formula contains proper physics around the atomic limit such as the kinetic exchange interactions and the RKKY interactions. The equation is extended to multiorbital systems in Sec. V. Based on the SCL formula for $\chi(\bm{q})$, we propose three levels of practical approximations in Sec. VI. An application to a two-orbital model is presented in Sec. VII to compare those approximation schemes. The paper is summarized in Sec. VIII.

II Fundamentals

II.1 Susceptibility in the DMFT

We first review how to evaluate the momentum dependent susceptibilities within the DMFT. The basic idea is presented in Ref. Jarrell, 1992, which showed that the antiferromagnetic (AFM) transition temperature covers the paramagnetic Mott transition in the Hubbard model. For early and recent reviews of the method, see Refs. Georges et al., 1996; Kuneš et al., 2017, respectively.

We define a local one-particle operator $O_{m\sigma,m^{\prime}\sigma^{\prime}}$ at site $i$ by

[TABLE]

where $m$ and $\sigma$ denote orbital and spin components, respectively. Any orbital or multipole operators that are invariant under the point-group symmetry are represented by a linear combinations of $O_{m\sigma,m^{\prime}\sigma^{\prime}}$. In this paper, we consider momentum- and frequency-dependent susceptibility corresponding to $O_{m\sigma,m^{\prime}\sigma^{\prime}}$. Regarding $O_{m\sigma,m^{\prime}\sigma^{\prime}}$ as a vector, the susceptibility is defined as a matrix $\bm{\chi}$ as

[TABLE]

where $\Omega$ denotes the bosonic Matsubara frequency, and abbreviations such as $1\equiv(m,\sigma)$ were introduced. Eigenvalues of $\bm{\chi}$ gives physical susceptibilities $\chi^{(\gamma)}(\bm{q},i\Omega)$. A phase transition can be detected by divergence of the static component $\chi^{(\gamma)}(\bm{q},0)$.

In the DMFT calculations, we treat the two-particle Green function, $X_{12,34}(i\omega,i\omega^{\prime};i\Omega)$, defined with two additional fermionic Matsubara frequencies, $i\omega$ and $i\omega^{\prime}$. The physical susceptibility $\chi$ is obtained by summing over the fermionic frequencies as

[TABLE]

For convenience, we introduce a matrix notation including $i\omega$ and $i\omega^{\prime}$ as

[TABLE]

Then, the BS equation which $\bm{X}(\bm{q},i\Omega)$ satisfies is represented as a matrix equation of the form

[TABLE]

where $\bm{X}_{0}(\bm{q},i\Omega)$ is defined by

[TABLE]

Here, $N$ denotes the number of lattice sites, and $G_{31}(\bm{k},i\omega)$ is the single-particle Green function within DMFT. The vertex part $\bm{\Gamma}_{\rm loc}(i\Omega)$ does not depend on momentum in accordance with the local self-energy in the DMFT. Hence, $\bm{\Gamma}_{\rm loc}(i\Omega)$ can be evaluated in the effective impurity problem through a similar BS equation

[TABLE]

where $\bm{X}_{\textrm{0,loc}}$ is defined by $\bm{X}_{\textrm{0,loc}}(i\Omega)=N^{-1}\sum_{\bm{q}}\bm{X}_{0}(\bm{q},i\Omega)$. Eliminating $\bm{\Gamma}_{\rm loc}$ in Eqs. (5) and (7), we obtain

[TABLE]

This is the equation we need to solve to compute $\bm{\chi}({\bm{q}})$. The size of the matrices in Eq. (8) is, in general, $4N_{\textrm{orb}}^{2}N_{\omega}\times 4N_{\textrm{orb}}^{2}N_{\omega}$, where 4 comes from spin, $N_{\textrm{orb}}$ denotes the number of orbitals, and $N_{\omega}$ is the number of the fermionic Matsubara frequencies. Solving the matrix equation as well as computing $\bm{X}_{\textrm{loc}}$ in the effective impurity problems becomes hard for multiorbital models at low temperatures.

Hereafter, we consider the static susceptibility, $\Omega=0$, and omit the argument $i\Omega$.

II.2 Exemplary results and aim of this paper

To illustrate the formalism so far, and to make it clear the purpose of this paper, we show exemplary results for the Hubbard model on a square lattice. The Hamiltonian reads

[TABLE]

where $\epsilon_{\bm{k}}=-2t(\cos k_{x}+\cos k_{y})$. The chemical potential is fixed at $\mu=U/2$ to consider a half-filled case. The band width is given by $W=8t$, and we use $t=1$ as the unit of energy. All computations were performed in a system of size $N=32\times 32$ with the periodic boundary condition.

Regarding the impurity solver, we used hybridization expansion algorithm of the continuous-time quantum Monte Carlo (CTQMC) method Werner et al. (2006); Gull et al. (2011) in the region $U/W\lesssim 1.5$. In the strong-coupling regime $U/W\gtrsim 1.0$, we used the Hubbard-I approximation to obtain a stable solution, avoiding statistical errors (see Appendix A for some remarks). In solving the matrix equation in Eq. (8), we introduced a cutoff for the fermionic frequencies. See Appendix B for the technical details.

Figure 1 shows the momentum dependence of the static spin susceptibility $\chi_{\bm{q}}^{\textrm{sp}}$. Here, $\chi_{\bm{q}}^{\rm sp}$ is defined as the correlation between the operator $O^{\textrm{sp}}=(O_{\uparrow\uparrow}-O_{\downarrow\downarrow})/\sqrt{2}$. From this figure, we confirm that fluctuations at $\bm{q}=(\pi,\pi)\equiv\bm{Q}$ is the strongest. Then, we fix $\bm{q}$ at $\bm{Q}$ and look into the $T$ dependence of $\chi_{\bm{Q}}^{\rm sp}$. Figure 2 shows the inverse of $\chi_{\bm{Q}}^{\rm sp}$ as a function of $T$ for $U/W=0.5$, 1, and 1.5 (red squares). In the DMFT, $\chi_{\bm{Q}}^{\rm sp}$ basically follows the Curie-Weiss law, $\chi_{\bm{Q}}^{\rm sp}\simeq C/(T-T_{\textrm{N}})$, due to the mean-field nature of the intersite correlations111In other words, the two-particle sum-rule $N^{-1}\sum_{\bm{q}}\chi_{\bm{q}}=\chi_{\textrm{loc}}$ is not fulfilled in the DMFT. To recover it, spatial correlations need to be included. See Ref. Rohringer et al., 2018 for a review of such extended theories.. The Curie constant $C$ agrees with that estimated from the localized spin, $C=1/2$, in the strong-coupling limit, and is reduced in the weak-coupling regime. The transition temperature determined from the divergence of $\chi_{\bm{Q}}^{\rm sp}$ is plotted as a function of $U$ in Fig. 3. Convergence to the Heisenberg limit $T_{\textrm{N}}\sim 4t^{2}/U$ is confirmed.

Two reference results are shown in Figs. 2 and 3 for comparison. One is the RPA computed by

[TABLE]

Here, $\chi_{0,\bm{q}}^{U=0}$ denotes the susceptibility for $U=0$ (no self-energy correction). The RPA yields $T_{\textrm{N}}$ that increases monotonically with increasing $U/W$ , and hence $T_{\textrm{N}}$ is considerably overestimated in the whole parameter region, as expected. Improvement to the RPA susceptibility is made by including correlations from DMFT calculations. We consider the following quantity, which we refer to as renormalized RPA (RRPA):

[TABLE]

This expression corresponds to replacing all $X$ with $\chi$ in Eq. (8) (namely, summation over the fermionic Matsubara frequencies are taken independently on each $X$). The RRPA formula takes two corrections into account compared with that in RPA: Firstly, $G$ in $\chi_{0}$ includes the self-energy correction, and secondly, the bare interaction $U$ is replaced with the effective one, $U_{\textrm{eff}}=\chi_{0,\rm loc}^{-1}-\chi_{\textrm{loc}}^{-1}$. Figure 3 demonstrates that the RRPA yields a reasonable value between the exact and RPA for $U\lesssim W$. However, it underestimates $T_{\textrm{N}}$ for $U\gtrsim W$ and does not gives the correct asymptotics for large $U$. Consequently, the RRPA results in a wrong phase diagram in which some phase is missing (orbz) in multiorbital cases (see Sec. VII for an explicit result). In the rest of this paper, we propose an alternative approach which is valid in the strong-coupling regime. The formula is as simple as the RPA formula and is executable with a low numerical cost comparable to the RRPA.

III Strong-coupling formula

In this section, we derive a simple susceptibility formula which gives the correct asymptotics in the large-$U$ limit. To avoid complexity and illustrate our main idea clearly, we restrict ourselves to a single-band model in this section. Extension to multiorbital systems will be presented in Sec. V.

III.1 Decoupling

In Eq (8), the frequency dependence in $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$ is particularly important in the strong coupling regime. The RRPA formula filters out the frequency dependence over $(i\omega,i\omega^{\prime})$ space and hence fails in reproducing the correct results for $U\gtrsim W$. Let us begin with observing $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$ to understand its structure in the $\omega$-$\omega^{\prime}$ plane. Figure 4(a) shows numerical results for $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$ in the spin channel (we represent it as $X_{\textrm{loc}}^{\textrm{sp}}$). It consists of sharp contribution at the diagonal, $\omega=\omega^{\prime}$, and a broad peak around the origin, $\omega=\omega^{\prime}=0$.

Our finding is that the broad contribution describes fluctuations of local degrees of freedom, which are relevant for $\chi_{\bm{q}}$, and the diagonal components are irrelevant. Furthermore, the broad peak in the $\omega$-$\omega^{\prime}$ plane can be decoupled as

[TABLE]

Mathematically, this approximation is carried out by considering $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$ as a matrix with respect to $\omega$ and $\omega^{\prime}$, and performing the singular value decomposition (SVD) to retain only the leading term:

[TABLE]

Here, $\{u_{i}(i\omega)\}$ and $\{v_{i}(i\omega)\}$ each is an orthogonal set and normalized by $\sum_{\omega}|u_{i}(i\omega)|^{2}=\sum_{\omega}|v_{i}(i\omega)|^{2}=1$. Provided that the phase factor of $u_{0}(i\omega)$ and $v_{0}^{*}(i\omega)$ are properly chosen, those functions satisfy $u_{0}(i\omega)=v_{0}^{*}(i\omega)$ (for details, see Appendix C). Thus, $\Phi(i\omega)$ can be obtained by

[TABLE]

The left panel of Fig. 5(a) shows the singular values $s_{i}$ of $X_{\textrm{loc}}^{\textrm{sp}}(i\omega,i\omega^{\prime})$ in the Hubbard model. It turns out that the first component $s_{0}$ is larger than the second component $s_{1}$ by 1–2 orders of magnitude. We found that $s_{0}$ increases in proportion to $1/T$. The ratio $s_{1}/s_{0}$ plotted as a function of $U$ in Fig. 5(b) demonstrates a rapid reduction of the higher-order contribution as $U$ increases and $T$ decreases. We conclude from these results that the decoupling in Eq. (13) is valid in the strong-coupling and low-temperature regime.

After the decoupling, $X_{\textrm{loc}}^{\textrm{sp}}(i\omega,i\omega^{\prime})$ in Fig. 4(a) is approximated into Fig. 4(b), and Fig. 4(c) shows the neglected part. It is clear from these figures that the broad peak is extracted by $s_{0}$, and the diagonal components are neglected, demonstrating that our central idea is supported in terms SVD.

For comparison, the charge channel, $X_{\textrm{loc}}^{\textrm{ch}}(i\omega,i\omega^{\prime})$, is plotted in Fig. 4(d). Here, the charge operator is defined by $O^{\textrm{ch}}=(O_{\uparrow\uparrow}+O_{\downarrow\downarrow})/\sqrt{2}$. In contrast to the spin channel, the dominant structure appears to be the diagonal component, and a broad peak around $\omega=\omega^{\prime}=0$ is missing. In this case, the singular values $s_{i}$ vary continuously as shown in the right panel of Fig. 5(a). The decoupling approximation, therefore, does not work for the charge channel.

The contrasting behaviors of $s_{i}$ in the spin and charge channels indicate that the singular peak at $s_{0}$ is indicative of residual degrees of freedom in the ground state. As the local spin fluctuations get stronger according to the Curie law, the singularity in $s_{0}$ becomes more conspicuous. Since long-range order takes place in such situations, the decoupling approximation is expected to be suitable for descriptions of phase transitions.

Let us finally derive an explicit relation between $s_{0}$ and $\chi_{\textrm{loc}}$. For this purpose, we use Eq. (3) as a “sum-rule” for $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$. Substituting Eq. (12) into Eq. (3), we obtain

[TABLE]

where the overline signifies an average defined by $\overline{\Phi}\equiv T\sum_{\omega}\Phi(i\omega)e^{i\omega 0^{+}}$. At low-$T$, $\overline{u}_{0}^{2}$ is proportional to $T$ (see Sec. IV.1), and we obtain $s_{0}\propto\chi_{\textrm{loc}}$.

III.2 The SCL formula

Now, we apply the decoupling in Eq. (12) to the BS equation in Eq. (8). To this end, we first rewrite Eq. (8) in the form which is suitable in the strong-coupling regime:

[TABLE]

where

[TABLE]

Equation (16) differs from the weak-coupling formula like the RRPA in the sense that the starting point is the local susceptibility and the non-local correction is taken into account by $\bm{\Lambda}_{\bm{q}}$. Expanding the right-hand side of Eq. (16) requires evaluation of the product $\bm{X}_{\textrm{loc}}\bm{\Lambda}_{\bm{q}}\bm{X}_{\textrm{loc}}$, which can be evaluated under the decoupling approximation as

[TABLE]

where we defined $\lambda_{\bm{q}}$ by

[TABLE]

Using this result, Eq. (16) can be expressed as $\bm{X}_{\bm{q}}\simeq(1-\lambda_{\bm{q}})^{-1}\bm{X}_{\textrm{loc}}$. The summations over $\omega$ and $\omega^{\prime}$ can be now taken explicitly to yield

[TABLE]

Here, the superscript SCL stands for strong-coupling limit. Equation (20) is our susceptibility formula derived in the strong-coupling limit.

It is convenient to express Eq. (20) in an RPA-like form. Defining a quantity $I_{\bm{q}}$ by

[TABLE]

Eq. (20) is rewritten as

[TABLE]

One can see that $I_{\bm{q}}$ represents effective intersite interactions, which give rise to the correction to the local susceptibility, $\chi_{\textrm{loc}}$. In the next section (Sec. IV), we will show that, indeed, $I_{\bm{q}}$ correctly reproduces the kinetic exchange and RKKY interactions in the Hubbard and periodic Anderson models, respectively.

Finally, it is important to confirm, in Eq. (21), that $\chi_{\textrm{loc}}$ is canceled with $\lambda_{\bm{q}}$ and $I_{\bm{q}}$ is independent of $\chi_{\textrm{loc}}$ (otherwise, $I_{\bm{q}}$ does not make physical sense). To see this, we use the relation in Eq. (15), which holds in the strong-coupling limit, and express $I_{\bm{q}}$ as 222Actually, this expression overestimates the interactions because Eq. (15) holds only in the strong-coupling limit and generally we have $T\chi_{\textrm{loc}}>s_{0}\overline{u}_{0}^{2}$. Nevertheless, we will use this expression in the SCL3 scheme (Sec. VI.2) because no estimation of $s_{0}$ is necessary.

[TABLE]

where we introduced a normalized function $\phi(i\omega)=\Phi(i\omega)/\overline{\Phi}$. It is now clear that, since $\phi(i\omega)$ is normalized, $I_{\bm{q}}$ does not depend on the magnitude of local fluctuations.

III.3 Numerical verification of the SCL formula

Let us apply the SCL formula, Eq. (22), to the Hubbard model to check the validity of the approximation. Figure 6 shows a comparison between $\chi_{\bm{q}}^{\textrm{SCL}}$ and $\chi_{\bm{q}}$ computed with the original BS equation. For the strong-coupling parameter, $U/W=1.5$, the two results agree quite well as expected. For weaker coupling, on the other hand, $\chi_{\bm{q}}^{\textrm{SCL}}$ is larger than $\chi_{\bm{q}}$, namely, overestimate fluctuations at high $T$. However, the deviation gets smaller as $T$ decreases, and the two results coincide at the transition temperature. This result indicates that the decoupling approximation does not cause a loss in quality of estimation of the transition.

The Néel temperature $T_{\textrm{N}}$ is shown in Fig. 7 as a function of $U/t$. Surprisingly, our formula, which was derived using the strong-coupling features of $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$, shows complete agreement with the correct $T_{\textrm{N}}$ even in the weak-coupling regime down to $U/W=2/8$. From this result, we conclude that our SCL formula, Eq. (22), provides the correct estimation of the transition temperature, though the fluctuations are overestimated for $T>T_{\textrm{N}}$. Multiorbital models will be examined in Sec. VII.

We close this section by discussing the effective nonlocal interactions, $I_{\bm{q}}$, which can be derived using the SCL formula by virtue of the simplified treatment of the BS equation. Figure 8 shows $I_{\bm{q}}^{\textrm{sp}}$ in the spin channel computed by Eq. (21). The positive (negative) values indicate enhancement (suppression) of the fluctuations, and the peak at M-point ($\bm{q}=\bm{Q}$) means antiferromagnetic interactions. The result for $U/W=1.5$ perfectly agrees with the analytic formula in the atomic limit, $I_{\bm{q}}^{\textrm{atom}}=-(4t^{2}/U)(\cos k_{x}+\cos k_{y})$, demonstrating the validity of the SCL formula. This will be proved analytically in the next section. For smaller $U$, the result shows deviation from the simple cosine function (especially at X-point), which indicates that the second-neighbor or even longer-range interactions are effectively active. Figure 9 shows the $T$-dependence of $I_{\bm{Q}}^{\textrm{sp}}$. It turns out that the value only weakly depends on $T$. Thus, it has been confirmed that $I_{\bm{q}}$ can be regarded as a physical effective interaction between the localized spins.

IV Atomic limit of the SCL formula

Our SCL formula in Eq. (22) is expected to capture correct physics in the strong-coupling regime. In this section, we verify analytically that the formula yields the well-known effective interactions emerging in the strong-coupling limit of two fundamental lattice models: the kinetic exchange interaction in the Hubbard model and the RKKY interaction in the periodic Anderson model.

IV.1 The function $u_{0}(i\omega)$ in the atomic limit

The expression for $I_{\bm{q}}$ in Eq. (23) includes two kinds of contribution: $\Lambda_{\bm{q}}(i\omega)$ takes account of the energy-band structure, and $\phi(i\omega)$ or $u_{0}(i\omega)$ takes local correlations into account. We first clarify the features of the local quantities in this subsection, and then derive the effective interactions in two lattice models in the succeeding subsections.

Figure 10 shows the numerical result for $u_{0}(i\omega)$ in a single-orbital atom with the Hubbard interaction. Here, the phase of $u_{0}(i\omega)$ is fixed so that $\sum_{\omega}u_{0}(i\omega)$ is positive real. The parameters $U$ and $\mu$ are chosen so that the ground state is the singly-occupied states without the particle-hole symmetry. We found numerically that $u_{0}(i\omega)$ at low temperatures can be perfectly fitted by the two-pole function of the form

[TABLE]

where $A$ is a real normalization factor determined by the condition $\sum_{\omega}|u_{0}^{\textrm{atom}}(i\omega)|^{2}=1$. The form of Eq. (24) can be interpreted in terms of single-particle excitations: The first contribution is due to the excitation to empty state, and the second due to the excitation to doubly occupied state. In the case with the particle-hole symmetry, namely $\mu=U/2$, the function $u_{0}^{\textrm{atom}}(i\omega)$ is reduced to a Lorentzian of the form $u_{0}^{\textrm{atom}}(\omega)=A(U/2)[\omega^{2}+(U/2)^{2}]$.

For the sake of later use, we derive explicit expressions of $A$ and $s_{0}$. For simplicity, we consider the low-temperature region, i.e., $T\ll\mu$ and $T\ll U-\mu$, in which the average $\overline{u}_{0}$, for example, is simplified into $\overline{u}_{0}\simeq A/2$. From the normalization condition of $u_{0}^{\textrm{atom}}(i\omega)$, $A$ is evaluated as

[TABLE]

where $r\equiv\mu/U$. Furthermore, using the Curie law $\chi_{\textrm{loc}}=1/2T$ ($1/2Tn_{\textrm{orb}}$ in the case with orbital degeneracy) in Eq. (15), we obtain

[TABLE]

This result is consistent with Fig. 5(a), which exhibits an increase of $s_{0}$ with decreasing $T$.

Finally, the function $\phi(i\omega)$ defined in connection with Eq. (23) is evaluated by $\phi(i\omega)=u_{0}(i\omega)/\overline{u}_{0}$, which becomes

[TABLE]

IV.2 Kinetic exchange interaction in the Hubbard model

The Hubbard model is reduced to the Heisenberg model in the strong-coupling limit. The exchange coupling constant is proportional to $t^{2}/U$. Applying Eq. (22) to the Hubbard model, we demonstrate below that $\chi_{\bm{q}}^{\textrm{SCL}}$ corresponds to the mean-field susceptibility of the effective Heisenberg model.

A key observation is that the quantity $\Lambda_{\bm{q}}(i\omega)$ defined in Eq. (17) becomes independent of $\omega$ in the strong-coupling limit. Given that the energy dispersion $\epsilon_{\bm{k}}$ takes the form $\epsilon_{\bm{k}}=2t\gamma_{\bm{k}}$ with $\gamma_{\bm{k}}\equiv\cos k_{x}+\cos k_{y}$, the quantity $\Lambda_{\bm{q}}(i\omega)$ becomes

[TABLE]

The detailed derivation is presented in Appendix D. The numerical result in Fig. 11(a) shows that $\Lambda_{\bm{q}}(i\omega)$ is well represented by Eq. (28) for $U/W=1.5$. We note that the quantities $X_{0,\bm{q}}(i\omega)$ and $X_{\textrm{0,loc}}(i\omega)$ [see Eq. (17)] do depend on $\omega$ [Fig. 11(a)], however, the $\omega$-dependence cancels out in $\Lambda_{\bm{q}}(i\omega)$.

Substitution of Eqs. (28) and (27) into Eq. (23) yields the coupling constant $I_{\bm{q}}$ of the form

[TABLE]

where we used $A^{2}=2UT$, assuming the half-filling, $\mu=U/2$. This expression corresponds to the kinetic exchange interaction between neighboring sites. The transition temperature is determined by the condition $\chi_{\textrm{loc}}I_{\bm{q}}=1$. Using $I_{\bm{q}}=8t^{2}/U$ and the Curie law $\chi_{\textrm{loc}}=1/2T$, we obtain the Néel temperature $T_{\textrm{N}}=4t^{2}/U$, which is consistent with the Heisenberg limit of the Hubbard model.

IV.3 RKKY interaction in the periodic Anderson model

We next consider the periodic Anderson model and derive the RKKY interaction between localized electrons. The Hamiltonian consists of band electrons written with the $c$ operator and localized (correlated) electrons written with $f$ operators:

[TABLE]

where $n_{fi\sigma}=f_{i\sigma}^{{\dagger}}f_{i\sigma}$. The ordinary derivation of the RKKY interaction is given as the fourth-order perturbation with respect to $V$ (see for example Ref. Shiba et al., 1999). We will derive a similar result based on the DMFT, starting from the SCL formula.

We consider the susceptibility $\chi(\bm{q})$ of $f$ electrons. Although there are two bands ($c$ and $f$) in this model, the BS equation is closed only with $f$ quantities because no two-body interaction is supposed between $c$ electrons Jarrell (1995). Hence, we simply regard all equations in Sec. III as those for $f$-components.

In the atomic limit of $I_{\bm{q}}$ in Eq. (23), $\phi(i\omega)$ is common between the Hubbard model and the periodic Anderson model. Hence, the difference between the two models arises from $\Lambda_{\bm{q}}(i\omega)$. In the strong-coupling limit, the analytical expression of $\Lambda_{\bm{q}}(i\omega)$ reads

[TABLE]

where

[TABLE]

The detailed derivation is presented in Appendix E. Figure 11(b) shows $\Lambda_{\bm{q}}(i\omega)$ in the periodic Anderson model, and demonstrates a striking difference with that in the Hubbard model. We note that no qualitative difference is seen in $X_{\textrm{0},\bm{q}}(i\omega)$, confirming that the quantity $\Lambda_{\bm{q}}(i\omega)$ is the relevant quantity which characterizes each model. In Fig. 11(b), we remark that $\Lambda_{\bm{q}}(i\omega)$ decays rapidly above $\omega\sim t$ [because $t$ is the only energy scale that enters in $g_{c,\bm{k}}(i\omega)$]. On the other hand, $\phi(i\omega)$ shows the power-law decay with the length $\omega\sim U$ [see Eq. (27) and Fig. 10]. Therefore, we can replace the $\omega$-dependence in $\Lambda_{\bm{q}}(i\omega)$ with $\delta_{\omega,0}$ if $U\gg t$:

[TABLE]

where $\chi_{\textrm{c},\bm{q}}$ is the static susceptibility of the conduction electrons defined by $\chi_{\textrm{c},\bm{q}}=T\sum_{\omega}[X_{\textrm{c},\bm{q}}(i\omega)-\bar{X}_{\textrm{c}}(i\omega)]$. Replacing $\Lambda_{\bm{q}}(i\omega)$ in Eq. (23) and using $\phi(0)=4/U$ [from Eq. (27)], we obtain

[TABLE]

where the Kondo exchange coupling $J_{\textrm{K}}$ is defined by $J_{\textrm{K}}=4V^{2}/U$. The above expression for $I_{\bm{q}}$ is consistent with the RKKY interaction.

As demonstrated above, our SCL formula correctly reproduces the RKKY interaction in $f$ electron systems and the kinetic exchange interaction in $d$ electron systems, depending on the form of $\Lambda_{\bm{q}}(i\omega)$ (delta function and constant, respectively). It provides a solid ground for applications of the SCL formula.

V Extension to Multiorbital systems

We now recall the orbital indices in $\bm{\chi}_{\textrm{loc}}$ and $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$, and consider multiorbital systems (or spin-orbital coupled systems if the spin-orbit coupling exists). In this case, we diagonalize $\bm{\chi}_{\textrm{loc}}$ with respect to the indices $\alpha\equiv(1,2)$ and $\alpha^{\prime}\equiv(3,4)$, and consider “eigen-modes” $\xi$ of local fluctuations as

[TABLE]

Using the unitary matrix $\bm{W}$, we transform $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ into the eigen-mode basis

[TABLE]

We neglect off-diagonal components ($\xi\neq\xi^{\prime}$), since they disappear when summations over $\omega$ and $\omega^{\prime}$ are taken. Thus, the discussion up to here for the single-orbital model can be applied to each eigen-mode $X_{\textrm{loc}}^{\xi}(i\omega,i\omega^{\prime})$.

As in Eq. (13), we apply SVD to $X_{\textrm{loc}}^{\xi}(i\omega,i\omega^{\prime})$ and retain the leading term as

[TABLE]

We show, in Fig. 12, an exemplary numerical result for the singular values $s_{i}^{\xi}$ in a two-orbital atom. The explicit form of the interactions will be presented in Sec. VII. Here, we only emphasize that the eight independent fluctuation modes are classified into two types. Figures 12(b), 12(c), 12(e), and 12(g) exhibit a singular behavior at $s_{0}$, as in the spin channel in the Hubbard model [left panel in Fig. 5(a)]. The singular enhancement of $s_{0}$ indicates the existence of local degrees of freedom, which gives rise to the Curie law in the local susceptibility. For such fluctuation modes, the decoupling in Eq. (38) is justified. On the other hand, the rest fluctuation modes [Figs. 12(a), 12(d), 12(f), and 12(h)] are continuous around $i=0$ as in the charge channel in the Hubbard model [right panel in Fig. 5(a)]. The decoupling is not valid for those modes. However, since such fluctuations are absolutely small, the decoupling does not affect principal results in $\bm{\chi}_{\bm{q}}$. Thus, the decoupling approximation in Eq. (38) is expected to provide a reasonable description in multiorbital models.

Using the decoupling in Eq. (38), we can represent Eq. (37) as follows:

[TABLE]

where

[TABLE]

The phases of $u_{0}^{\xi}(i\omega)$ and $v_{0}^{\xi}(i\omega)$ have been fixed so that $u_{0}^{\xi}(i\omega)={v_{0}^{\xi}}^{*}(i\omega)$ (for details, see Appendix C). The frequency dependence has thus been decoupled as in Eq. (12), while the spin-orbital indices, $\alpha$ and $\alpha^{\prime}$, are retained to be evaluated by the matrix product.

Once the decoupling in Eq. (39) is applied, the BS equation in Eq. (16) is evaluated to yield

[TABLE]

where we defined a multiorbital extension of $\lambda_{\bm{q}}$ in Eq. (19) by

[TABLE]

We can now take the summations over $\omega$ and $\omega^{\prime}$ independently in Eq. (41), and obtain the SCL formula for the spin and orbital susceptibility

[TABLE]

This is the multiorbital extension of Eq. (20). Here, we used the relation $\overline{\bm{\Phi}}^{2}\simeq T\bm{\chi}_{\textrm{loc}}$ with $\overline{\bm{\Phi}}\equiv T\sum_{\omega}\bm{\Phi}(i\omega)e^{i\omega 0^{+}}$ [multiorbital version of Eq. (15)].

As in the case of the single-orbital formula, we can represent $\bm{\chi}_{\bm{q}}^{\textrm{SCL}}$ in an RPA-like expression. We define

[TABLE]

and then rewrite Eq. (43) as

[TABLE]

$\bm{I}_{\bm{q}}$ represents the intersite spin and orbital interactions. In the strong-coupling regime, we can eliminate $\bm{\chi}_{\textrm{loc}}$ in the expression of $\bm{I}_{\bm{q}}$. Following the procedure leading to Eq. (23), we can show

[TABLE]

where $\bm{\phi}(i\omega)\equiv\bm{\Phi}(i\omega)\overline{\bm{\Phi}}^{-1}$. Since $\bm{\phi}(i\omega)$ is normalized, $\bm{I}_{\bm{q}}$ has no explicit dependence on local fluctuations, and reflects lattice properties predominantly. We note again that $\bm{I}_{\bm{q}}$ describes both the kinetic exchange process and the RKKY interactions depending on the functional form of $\bm{\Lambda}_{\bm{q}}(i\omega)$ as demonstrated in Sec. IV.

If the $\xi$-dependence in $\phi^{\xi}(i\omega)$ is neglected as a rough approximation (see Sec. VI.2), $\bm{I}_{\bm{q}}$ is reduced to a simpler expression

[TABLE]

which is similar to Eq. (23).

VI Approximation scheme

Our susceptibility formulas [Eq. (22) for single-orbital cases and Eq. (45) for multiorbital cases] simplify the process of computing $\bm{\chi}_{\bm{q}}$ in the DMFT. In particular, it has been shown, in the single-orbital Hubbard model, that $T_{\textrm{N}}$ computed using the formula shows a complete agreement with the rigorous solution of the BS equation. This result indicates that the decoupling of $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ essentially does not worsen the estimation of transition temperatures. In a practical point of view, however, the computational cost is still not low enough, since the decoupling in Eq. (13) requires as much information as in solving the BS equations. Simpler calculations of $\bm{\chi}_{\bm{q}}$ are achieved by evaluating the decoupled functions $\bm{\Phi}(i\omega)$ or $\bm{\phi}(i\omega)$ without SVD of $X_{\textrm{loc}}^{\xi}(i\omega,i\omega^{\prime})$.

In order to reduce the computational cost for $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$, we propose three levels of approximation in evaluating $\bm{\Phi}(i\omega)$:

SCL1

Using SVD [Eq. (13)]. Full information of $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ is required within a certain cutoff frequency. Results are already shown in Sec. III.3.

SCL2

Estimating $u_{0}^{\xi}(i\omega)$ from one-dimensional data of $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$. For example, $i\omega$ is varied with fixing $i\omega^{\prime}=i\omega_{0}$.

SCL3

Assuming an analytical form of $u^{\xi}_{0}(i\omega)$ which is valid in the atomic limit. No calculation of $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ is required.

The difference in the necessary information of $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ is schematically depicted in Fig. 13. The cost of computing $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ is $\mathcal{O}(N_{\omega}^{2})$, $\mathcal{O}(N_{\omega}^{1})$, and $\mathcal{O}(1)$ for SCL1, SCL2, and SCL3, respectively. Here, $N_{\omega}$ is the number of fermionic Matsubara frequencies in each axis.

VI.1 SCL2: A crude estimation of $u_{0}(i\omega)$

In the second scheme, SCL2, only a part of data in $\omega$-$\omega^{\prime}$ plane is used. This approximation is based on the fact that, if the decoupling in Eq. (12) is exact, one-dimensional data (cut of the plane) is sufficient for evaluation of $u_{0}(i\omega)$. Actually, since the result depends on the location of the cut, one can improve accuracy by considering more than one cut and taking the average over different estimations. Technical details are presented in Appendix F for a practical implementation.

The phase diagram obtained by the SCL2 scheme is shown in Fig. 14(a). The value of $T_{\textrm{N}}$ turns out to be underestimated, in particular, in the weak-coupling regime, $U/W\lesssim 1.0$. The deviation becomes smaller as $U$ increases, and the SCL2 reproduces $T_{\textrm{N}}$ almost perfectly in the strong-coupling regime, $U/W\gtrsim 1.5$.

The fact that the SCL2 scheme underestimates $T_{\textrm{N}}$ can be explained in terms of a variational principle as follows. We first remark that the square of the singular value, $s_{i}^{2}$, corresponds to the eigenvalue of $\bm{X}_{\textrm{loc}}\bm{X}_{\textrm{loc}}^{{\dagger}}$. Hence, the decoupling approximation, Eq. (13), is represented by $s_{0}^{2}\simeq\langle u_{0}|\bm{X}_{\textrm{loc}}\bm{X}_{\textrm{loc}}^{{\dagger}}|u_{0}\rangle$. The right-hand side becomes maximum when $u_{0}(i\omega)$ is the true eigenvector of $\bm{X}_{\textrm{loc}}\bm{X}_{\textrm{loc}}^{{\dagger}}$. In other words, an approximate evaluation of $u_{0}(i\omega)$ results in a smaller estimation of $s_{0}$, namely, the underestimate of $\chi_{\bm{q}}$.

This fact gives us “controlability” in applications of the SCL2 scheme. Since the SCL2 scheme underestimates all fluctuation modes, a competition among different modes would still be captured correctly. Indeed, the result for the two-orbital model in Sec VII supports the validity of the SCL2 approximation.

VI.2 SCL3: Two-pole approximation

In the single-orbital model, the function $u_{0}(i\omega)$ is reduced to a simple two-pole function (Lorentzian) in the atomic limit, as discussed in Sec. IV.1. This results from the fact that the intermediate empty and doubly occupied states are both singlet. In multiorbital models, on the other hand, $u_{0}(i\omega)$ is not of Lorentzian, since intermediate states have multiplet. Nevertheless, the two-pole function gives a reasonable approximation for relevant fluctuation modes that follow the Curie law. Detailed analysis in a two-orbital model is presented in Appendix G. For this reason, we adopt the two-pole function also in multiorbital cases and construct an approximation that does not require explicit computations of $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$.

We generalize the two-pole function in Eq. (27) into multiorbital cases as333This function fulfills the normalization condition $\overline{\phi^{\xi}}=1$ only at low-$T$ of $T\ll E_{\pm}^{\xi}$. However, it does not affect an estimation of transition temperatures. Furthermore, the weights of the two terms can actually be different from each other (keeping the sum at 2), but we do not introduce such a parameter for simplicity.

[TABLE]

where $E_{-}^{\xi}=E_{n-1}-E_{n}>0$ and $E_{+}^{\xi}=E_{n+1}-E_{n}>0$. Here, $E_{n}$ denotes a “typical” value of $n$-particle states, e.g., the average or the lowest energy. In general, $E_{\pm}^{\xi}$ depends on $\xi$, because the dominant contribution to $X_{\textrm{loc}}^{\xi}(i\omega,i\omega^{\prime})$ occurs from different $(n\pm 1)$-particle states depending on $\xi$. Once the parameters $E_{\pm}^{\xi}$ are fixed, we can readily evaluate $\bm{\chi}_{\bm{q}}^{\textrm{SCL}}$ using Eqs. (45) and (46), which we refer to as SCL3 approximation. In practical applications, it is convenient to replace $E_{\pm}^{\xi}$ with the lowest excitation energy $E_{\pm,\textrm{min}}$ estimated in the local Hamiltonian. For details, see the example in Appendix G and Sec. VII.

The accuracy of the SCL3 scheme is examined in Fig. 14(b) for the single-orbital model, in which $E_{\pm}^{\xi}$ can be uniquely fixed at $E_{\pm}^{\xi}=U/2$ for the half-filling. The transition temperature turns out to be slightly overestimated for $U/W\gtrsim 0.5$ and converges to the strong-coupling tail. The deviation is large in the weak-coupling regime of $U/W\lesssim 0.5$, since the analytic expression in Eq. (48) is justified only for large-$U$.

Though the SCL3 approximation gives reasonable results for the single-orbital model, the situation is rather complicated in multiorbital models. The SCL3 scheme corresponds to an approximation which neglects splitting of energy levels of excited multiplet, and may result in a misjudgment of competing orders. See the next section for a test in a multiorbital model.

VII Application: Two-orbital model

VII.1 Model and background

For a demonstration of the multiorbital version of the strong-coupling formula, Eq. (45), we consider a two-orbital model with an energy splitting $\Delta$. Letting the spin index $\sigma=\uparrow$, $\downarrow$, and the orbital index $\tau=\textrm{a}$, b, the Hamiltonian reads

[TABLE]

The kinetic energy term is assumed to be orbital-diagonal. We consider the tight-binding band on the square lattice in common with the previous example in Eq. (9). The second term expresses the energy splitting $\Delta$ between the two bands. The rest terms are the Slater-Kanamori interactions. We consider the half-filled case, $n=2$, and vary $\Delta$ keeping other parameters, $U$, $U^{\prime}$, and $J$. We fix $U$ at $U/W=1.5$, and use the Hubbard-I solver, which is expected to be reasonable at this parameter as demonstrated in the case with the single-band model (see for example Fig. 3). The Hund’s coupling $J$ is fixed at $J=U/4$, and $U^{\prime}$ at $U^{\prime}=U-2J$ to keep the rotational symmetry of the orbitals when $\Delta=0$.

This model has been discussed in the context of the spin-state transition Werner and Millis (2007); Suzuki et al. (2009); Kuneš and Křápek (2011) and an excitonic insulator (EI) Kaneko et al. (2012); Kuneš and Augustinský (2014); Kuneš (2014); Nasu et al. (2016); Hoshino and Werner (2016). Let us first review the basic properties of this model. When $\Delta\ll W$, the two orbitals are almost degenerate, and the spin-triplet state is favored according to the Hund’s coupling. The two orbitals are equally occupied as shown in Fig. 15. When $\Delta\gg W$, on the other hand, only the lower band is occupied and the spin state becomes singlet. It is known that the change between the two spin-states becomes more abrupt as $T$ is lowered Werner and Millis (2007) as seen around $\Delta\simeq W$ in Fig. 15. Hence, large orbital fluctuations are expected. As a results of competing spin and orbital fluctuations, a spin-orbital coupled state called EI state emerges in this region Kaneko et al. (2012); Kuneš and Augustinský (2014). We examine, in the following, whether this non-trivial ordered state can be reproduced by our simplified calculation scheme.

VII.2 Susceptibilities

There are 16 fluctuation modes represented by the operator $O_{\tau\sigma,\tau^{\prime}\sigma^{\prime}}$ defined in Eq. (1). It is convenient to represent the spin and orbital indices in terms of the Pauli matrix $\sigma^{\xi}$ ($\xi=0,x,y,z$) as

[TABLE]

These operators are classified into eight classes as summarized in Table 1. The number of degeneracy is either 1 or 3, where threefold degeneracy is due to the spin part, ($\sigma^{x}$, $\sigma^{y}$, $\sigma^{z}$), while no extra degeneracy takes place because the orbital indices are all independent. For reference, correspondence to the irreducible representations in the cubic point-group symmetry is also listed J. Ohkawa (1983); Shiina et al. (1997). The indices g (u) and $+$ ($-$) stand for time-reversal even (odd). Here, the odd time-reversal symmetry comes from $\sigma^{x}$, $\sigma^{y}$, $\sigma^{z}$, and $\tau^{y}$. The third to fifth classes (labeled with orb) are kinds of orbital orders, and the 6th to 8th (labeled with sp-orb) are spin-orbital coupled orders with broken time-reversal symmetry. In particular, the 3rd operator (orbz) is the same as the level splitting, and its staggard order corresponds to high-spin–low-spin (HS/LS) order. The 7th operator (sp-orb3x) corresponds to the EI ordered state.

The left panel of Fig. 16(a) shows the $\bm{q}$-dependence of $\chi_{\bm{q}}$ obtained by solving the original BS equation in Eq. (8). We chose the parameter $\Delta/W=1.1$, where multiple fluctuation modes are competing. It turns out that $\bm{\chi}_{\bm{q}}$ exhibits enhancement at M-point [$\bm{q}=(\pi,\pi)\equiv\bm{Q}$]. The leading fluctuations are the EI order, HS/LS order, AFM order, and threefold orbital fluctuations (orb3). Other fluctuation modes exhibit no distinguished enhancement (not shown in the figure). The transition temperatures are determined from the divergence of $\chi_{\bm{Q}}$ [the center panel of Fig. 16(a)]. The obtained phase diagram [the right panel of Fig. 16(a)] includes AFM phase for $\Delta/W\lesssim 1.0$, and EI and HS/LS phases for $\Delta/W\gtrsim 1.0$. The HS/LS order appears only in the limited region at finite-$T$ because this state has degeneracy without additional symmetry breaking. The orb3 order is not realized since other transition takes place at higher $T$.

Let us now examine the accuracy of other approximation schemes. Figure 16(b) shows the SCL1 result. It is quite similar to Fig. 16(a) with only a little enhancement of the transition temperatures, justifying the validity of our SCL formula even for the multi-orbital systems.

The SCL2 result in Fig. 16(c) also reproduces the phase diagram very well. The qualitative feature such as the severe competition between EI and HS/LS orders are correctly reproduced. Quantitatively, the transition temperatures are slightly underestimated, in consistent with the conclusion in Sec. VI.1.

The result in the SCL3 scheme is shown in Fig. 16(d). In this calculation, we neglected the $\xi$-dependence of $\phi^{\xi}(i\omega)$ for simplicity, and employed Eq. (47). The parameter $E_{\pm}^{\xi}$ in Eq. (48) were replaced with the lowest excitation energy $E_{\pm,\textrm{min}}$ determined from the eigenvalues of the local Hamiltonian. Despite these approximations, we see the diverging susceptibility for the four enhanced modes. However, the estimate of the transition temperature is less accurate than those by SCL1 and SCL2. The phase diagram shows that the AFM fluctuations are largely overestimated, while the fluctuations for the HS/LS order is underestimated so that the HS/LS order is missing (covered by the EI order). The source of these artifacts is discussed in the next subsection in terms of the effective nonlocal interactions.

Finally, Fig. 16(e) shows the RRPA result. Although the RRPA reproduces the divergence of susceptibility in AFM, EI, and orb3 orders, the HS/LS susceptibility are largely underestimated and does not diverge at finite temperature. This artifact suggests that, as in the single-orbital case, the usage of the RRPA in the strong-coupling regime is not justified.

VII.3 Effective nonlocal interactions

In our SCL formula, it is possible to derive the effective nonlocal interactions $I_{\bm{q}}$, which enters in the RPA-like susceptibility formula in Eq. (45). Figure 17 shows $I_{\bm{q}}$ computed in the three approximate schemes. For all schemes, the momentum dependence exhibits the simple cosine form, $\cos k_{x}+\cos k_{y}$, as in the strong-coupling regime of the single-band Hubbard model ($U/W=1.5$ in Fig. 8). It indicates that the effective interaction works only between nearest-neighbor sites.

With the results of $I_{\bm{q}}$ in Fig. 17, we can understand the results of each SCL scheme in Fig. 16. For SCL1 and SCL2, the HS/LS mode has the strongest interaction, which helps the realization of the HS/LS order around $\Delta/W=1.1$. On the other hand, the SCL3 scheme exhibits no mode-dependence in $I_{\bm{q}}$, and hence underestimates the HS/LS transition temperature in Fig. 16(d).

The artificial degeneracy of $I_{\bm{q}}$ in the SCL3 scheme originates from the approximation that we neglected the $\xi$-dependence of $E_{\pm}^{\xi}$ (instead $E_{\pm,\textrm{min}}$ is used for all modes) in Eq. (48). According to the SVD analysis in Appendix G, $E_{\pm}^{\textrm{HS/LS}}$ is actually smaller than the others, so replacing $E_{\pm}^{\textrm{HS/LS}}$ with $E_{\pm,\textrm{min}}$ reduces the HS/LS fluctuation. On the contrary, the actual value of $E_{\pm}^{\textrm{spin}}$ is about two times larger than $E_{\pm,\textrm{min}}$, meaning that the use of $E_{\pm,\textrm{min}}$ results in an overestimate of the spin exchange interactions. The above analysis of $I_{\bm{q}}$ indicates that a proper estimation of $E_{\pm}^{\xi}$ is needed when the SCL3 scheme is applied to competing orders.

VIII Summary

The BS equation for the calculations of $\bm{\chi}_{\bm{q}}$ in the DMFT was reduced to an RPA-like form called the SCL formula [Eq. (22) for single-orbital cases and Eq. (45) for multiorbital cases]. Its derivation is based on the fact that the two-particle Green function $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ can be decoupled [Eqs. (12) and (39)] if local moments are well-defined and there exist local fluctuations which give rise to the Curie law. Since a long-range order often occurs in such situations, the SCL formula is suitable for the use of investigations of phase transitions. Indeed, the SCL formula is found to yield an accurate transition temperature over a surprisingly wide range of parameters, though $\bm{\chi}_{\bm{q}}$ is overestimated at higher temperatures in the weak-coupling region.

So far, the bottleneck in computing $\bm{\chi}_{\bm{q}}$ was twofold: (i) calculations of $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ in the effective local problem and (ii) solving the BS equation with large matrices. The SCL formula provides the correct phase diagram in accuracy comparable to the BS equation [see SCL1 results in Figs. 7 and 16(b)], and therefore settles the latter issue. Furthermore, approximate estimations of the function $\bm{\Phi}(i\omega)$ or $\bm{\phi}(i\omega)$ (SCL2 and SCL3 schemes) relieve the former issue. The cost of computing $\bm{X}_{\textrm{loc}}(i\omega,i\omega^{\prime})$ is reduced from $\mathcal{O}(N_{\omega}^{2})$ to $\mathcal{O}(N_{\omega}^{1})$ in the SCL2 scheme and $\mathcal{O}(1)$ in the SCL3 scheme, where $N_{\omega}$ is the number of fermionic Matsubara frequencies. We note that, in the SCL3 scheme, the mode-dependence of the effective excitation energy, $E_{\pm}^{\xi}$, is relevant when different types of fluctuation modes are competing.

By combining the three levels of approximation in the SCL formula and the RRPA formula in the weak-coupling side, one can evaluate $\bm{\chi}_{\bm{q}}$ in the DMFT without solving the BS equation. It will enable systematic investigations of complicated spin-orbital orders in multiorbital systems, in particular, realistic $d$ and $f$-electron materials within the DFT+DMFT framework.

Acknowledgements.

We acknowledge useful discussions with S. Ishihara, M. Naka, J. Kuneš, H. Ikeda, Y. Ōno, T. Yamada, Y. Iizuka. The numerical calculations were performed using the following open source libraries and softwares: TRIQSParcollet et al. (2015), DFTToolsAichhorn et al. (2016), pomerolAntipov et al. (2017), pomerol2triqs, and pomerol2triqs1.4. We were supported by JSPS KAKENHI Grant No. 18H01158. J.O. was supported by JSPS KAKENHI Grant No. 18H04301 (J-Physics). K.Y. was supported by Building of Consortia for the Development of Human Resources in Science and Technology, MEXT, Japan. H.S. was supported by JSPS KAKENHI Grant No. 16K17735. Y.N. was supported by JSPS KAKENHI Grant Nos. 16H06345 and 17K14336.

Appendix A Hubbard-I approximation

In practical computations for multiorbital models with complicated interactions, it is difficult to solve the effective impurity problem at arbitrary temperature and interaction strength without approximation. Furthermore, the self-consistency loop of the DMFT may become unstable in the strong-coupling regime. In this situation, it is useful to consider the atomic limit expressed by

[TABLE]

and plug only the one-body levels into the impurity problem. This is called the Hubbard-I approximation in the context of Mott gap in the single-particle excitation spectra.

With the approximation in Eq. (51), we can compute the local self-energy $\Sigma_{\rm atom}(i\omega)$ and the two-particle Green function $X_{\textrm{atom}}(i\omega,i\omega^{\prime};i\Omega)$ by diagonalizing the atomic Hamiltonian explicitly. We used pomerol implementationAntipov et al. (2017) in our calculations.

Some remarks are made on the atomic-limit calculations for two-particle quantities. Because of the absence of self-consistent determination of $\Delta(i\omega)$, the self-consistency condition of DMFT

[TABLE]

is not fulfilled. Here, $G_{\textrm{imp}}(i\omega)$ denotes the Green function computed in the impurity problem. It leads that a similar equality for the bare two-particle Green function

[TABLE]

is also violated, though the self-consistently determined solution of the DMFT does have this equality. Here, $\bm{X}_{0,\textrm{imp}}(i\Omega)$ is given by $[\bm{X}_{0,\textrm{imp}}]_{\omega,\omega^{\prime}}(i\Omega)=-\delta_{\omega\omega^{\prime}}G_{\textrm{imp}}(i\omega)G_{\textrm{imp}}(i\omega+i\Omega)$. Thus, there is a choice between $\bm{X}_{0,\textrm{imp}}(i\Omega)$ and $\bm{X}_{0,\textrm{loc}}(i\Omega)$ in solving the BS equation in Eq. (8). We found that using $\bm{X}_{\rm 0,loc}(i\Omega)$ gives reasonable results, because the difference between $\bm{X}_{0}(\bm{q},i\Omega)$ and $\bm{X}_{\rm 0,loc}(i\Omega)$ plays an essential role in the $\bm{q}$-dependence of the susceptibility.

Appendix B Frequency cutoff

In practical computations of the BS equation, one needs to introduce cutoff in the fermionic frequencies. Here, we present a practical expression that is suitable for introducing cutoff. We begin with the BS equation in the form of Eq. (16), which is rewritten as follows:

[TABLE]

Here, $\Delta\bm{X}_{\bm{q}}$ describes the deviation with respect to $\bm{X}_{\textrm{loc}}$. Taking the summations over the fermionic Matsubara frequencies, we obtain

[TABLE]

The infinite summations in the first term have been taken exactly. We now introduce a frequency cutoff for the summations in the second term to evaluate the deviation from $\bm{\chi}_{\textrm{loc}}$.

Appendix C Relations between $u_{0}(i\omega)$ and $v_{0}(i\omega)$

In this section, we derive relations between $u_{0}(i\omega)$ and $v_{0}(i\omega)$ from symmetry properties of two-particle Green functions. We first use the relation

[TABLE]

which implies that $u_{0}(i\omega)$ and $v_{0}(i\omega)$ are basically equivalent. Replacing $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$ with the decoupling formula in Eq. (13), we can prove $u_{0}(i\omega)=Cv_{0}^{*}(i\omega)$, where $C$ is a constant. Furthermore, the normalization condition for $u_{0}(i\omega)$ and $v_{0}(i\omega)$ leads to $|C|=1$, and thus the relation is expressed as

[TABLE]

where $\theta$ is a real number.

Next, we use another relation that connects $i\omega$ and $-i\omega$

[TABLE]

The same procedure as above leads to a similar relation between $u_{0}(i\omega)$ and $u_{0}^{*}(-i\omega)$, namely, $u_{0}(i\omega)=e^{i\theta^{\prime}}u_{0}^{*}(-i\omega)$. Furthermore, combining with Eq. (58), we can prove $\theta=\theta^{\prime}$, and thus obtain

[TABLE]

It is also possible to eliminate the phase factor by combining Eqs. (58) and (60) to yield

[TABLE]

which does not depend on the choice of $\theta$.

In practical calculations, it is convenient to fix the phase so that $\theta=0$. To this end, we define $\widetilde{u}_{0}(i\omega)\equiv u_{0}(i\omega)e^{-\theta/2}$ and $\widetilde{v}_{0}(i\omega)\equiv v_{0}(i\omega)e^{-\theta/2}$, and thus $\widetilde{u}_{0}(i\omega)$ and $\widetilde{v}_{0}(i\omega)$ follow the above equations with $\theta=0$, namely,

[TABLE]

Appendix D Derivation of $\Lambda_{\bm{q}}(i\omega)$ in the Hubbard model

We derive the analytical expression for $\Lambda_{\bm{q}}(i\omega)$ in Eq. (28). The single-particle Green function in the Hubbard model is given by

[TABLE]

In the atomic limit, the local Green function $G_{\textrm{loc}}(i\omega)$ is

[TABLE]

We expand $G_{\bm{k}}(i\omega)$ with respect to $\epsilon_{\bm{k}}$ and retain the first correction

[TABLE]

Using this expression, $X_{\bm{q}}^{0}(i\omega;i\Omega)$ is evaluated as

[TABLE]

Here, we used $\langle\epsilon_{\bm{k}}\rangle_{\bm{k}}=0$. Inserting this approximated result into $\Lambda_{\bm{q}}(i\omega;i\Omega)$ in Eq. (17), we obtain

[TABLE]

Taking the average over $\bm{k}$ and using $\langle\gamma_{\bm{k}}\gamma_{\bm{k}+\bm{q}}\rangle_{\bm{k}}=\gamma_{\bm{k}}/2$ for a nearest-neighbor hopping model, we obtain Eq. (28).

Appendix E Derivation of $\Lambda_{\bm{q}}(i\omega)$ in the periodic Anderson model

The analytical expression for $\Lambda_{\bm{q}}(i\omega)$ in Eq. (31) is derived in this appendix. The single-particle Green function in the periodic Anderson model is given by

[TABLE]

where $g_{\textrm{c},\bm{k}}(i\omega)$ is the conduction electron Green function defined by $g_{\textrm{c},\bm{k}}(i\omega)=1/(i\omega+\mu-\epsilon_{\bm{k}})$. The local Green function $G_{\textrm{loc}}$ is given by

[TABLE]

In the atomic limit, $\Delta(i\omega)$ is approximated into $\Delta(i\omega)\simeq\langle V_{\bm{k}}^{2}g_{\textrm{c},\bm{k}}(i\omega)\rangle_{\bm{k}}$, and we expand $G_{\bm{k}}^{-1}(i\omega)=G_{\textrm{loc}}^{-1}(i\omega)-[V_{\bm{k}}^{2}\epsilon_{\bm{k}}-\Delta(i\omega)]$ around $G_{\textrm{loc}}(i\omega)$ as

[TABLE]

Then, $X_{\bm{q}}^{0}(i\omega;i\Omega)$ is evaluated as

[TABLE]

Inserting this expression into $\Lambda_{\bm{q}}(i\omega;i\Omega)$ in Eq. (17) and retaining the term of order $V^{4}$, we obtain Eq. (31).

Appendix F Technical details of SCL2 scheme

Let us assume that the decoupling in Eq. (12) is exactly satisfied, namely,

[TABLE]

Then, $u_{0}(i\omega)$ can be determined from a one-dimensional data of $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$ as shown in Fig. 13: $u_{0}(i\omega)\propto X_{\textrm{loc}}(i\omega,-i\omega_{0})$, where $\omega_{0}$ is, e.g., the lowest Matsubara frequency. In practice, however, a good estimation is not achieved in this way, because actual data of $X_{\textrm{loc}}(i\omega,i\omega^{\prime})$ has a structure around the diagonal $\omega=\omega^{\prime}$, which is not included in $u_{0}(i\omega)$. For better estimation of $u_{0}(i\omega)$, we need to extract only the broad structure away from the diagonal. For this reason, we use only the half data on $\omega>0$ side as depicted in Fig. 18, and the negative side is recovered using the symmetry relation in Eq. (60). In the following, we describe how to construct $u_{0}(i\omega)$ in this way.

We first determine $\theta$, which is necessary to use Eq. (60). To this end, we input two data points, P and $\textrm{P}^{\prime}$, in Fig. 18 and define

[TABLE]

Here, $\omega_{1}\neq\omega_{0}$ is, e.g., the second lowest Matsubara point. Then, we shift the phase of $X_{\textrm{loc}}(i\omega,-i\omega_{0})$ by $-\Theta$ to define

[TABLE]

where $C=1/s_{0}|v_{0}^{*}(-i\omega)|$ is a real number. The negative side, $\omega<0$, can be evaluated using the symmetry $\widetilde{u}_{0}(-i\omega)=\widetilde{u}_{0}^{*}(i\omega)$, which can be derived from Eq. (60). The factor $C$ is fixed from the normalization condition. Finally, the singular value $s_{0}$ is determined using the relation $X_{\textrm{loc}}(i\omega_{0},-i\omega_{0})=s_{0}|u_{0}(i\omega_{0})|^{2}$.

Appendix G The function $u_{0}(i\omega)$ in a two-orbital model

In this Appendix, we present explicit results for the function $u_{0}^{\xi}(i\omega)$ in the two-orbital model, and provide a grounding of the two-pole approximation introduced in Sec. VI.2.

We first show the mode-dependence of the ratio $s^{\xi}_{1}/s^{\xi}_{0}$ in Fig. 19(a). The data are taken from Fig. 12 (the lowest temperature, $T=0.1$). The decoupling approximation is justified if $s^{\xi}_{1}/s^{\xi}_{0}$ is small, and four modes, i.e., spin, orbz, orb3, and sp-orb3x, apply to this case as expected.

The function $u_{0}^{\xi}(i\omega)$ is shown for those four modes in Fig. 20. The lines show fitting results by the two-pole function in Eq. (48), demonstrating that reasonable fitting was achieved. The fitting parameter $E_{+}=E_{-}$ are depicted in Fig. 19(b). The horizontal line indicates the lowest excitation energy $E_{\pm,\textrm{min}}\approx 3.40$ from $n=2$ to $n=2\pm 1$. The fitting results excepting the spin mode turn out to agree with $E_{\pm,\textrm{min}}$, indicating that replacing $E_{\pm}$ with $E_{\pm,\textrm{min}}$ provides a good approximation. Regarding the spin mode, this replacement is expected to result in an overestimate of fluctuations.

Finally, we verify Eq. (15), which relates $s_{0}^{\xi}$ to $\chi_{\textrm{loc}}^{\xi}$. Figure 19(c) compares the left- and right-hand sides of Eq. (15), demonstrating that Eq. (15) is satisfied in reasonable accuracy. All the results above support the SCL3 approximation scheme as a simple and reasonable approximation for solving the BS equation.