Rotationally invariant slave-boson and density matrix embedding theory: A unified framework and a comparative study on the 1D and 2D Hubbard Model

TL;DR

This paper compares the accuracy of RISB and DMET methods for the 1D and 2D Hubbard models, introduces a unified framework for both, and discusses potential for future generalizations.

Contribution

It provides a unified computational framework for RISB and DMET, enabling direct comparison and new insights into their similarities and differences.

Findings

RISB and DMET have similar accuracy and performance

Unified framework allows consistent implementation of both methods

Facilitates development of generalized approaches

Abstract

We present detailed benchmark ground-state calculations of the one- and two-dimensional Hubbard model utilizing the cluster extensions of the rotationally invariant slave-boson (RISB) mean-field theory and the density matrix embedding theory (DMET). Our analysis shows that the overall accuracy and the performance of these two methods are very similar. Furthermore, we propose a unified computational framework that allows us to implement both of these techniques on the same footing. This provides us with a new line of interpretation and paves the ways for developing systematically new generalizations of these complementary approaches.

Figures

Tables

Table 1. Table 1: Energy E / t 𝐸 𝑡 E/t for DMET and RISB in the AFM phase of the 2D Hubbard model at half-filled n = 1 𝑛 1 n=1 with the nearest neighbor hopping for N c = 2 , 4 , 6 subscript 𝑁 𝑐 2 4 6 N_{c}=2,\ 4,\ 6 at U = 2 ​ t , 4 ​ t , 6 ​ t , 8 ​ t , 12 ​ t 𝑈 2 𝑡 4 𝑡 6 𝑡 8 𝑡 12 𝑡 U=2t,\ 4t,\ 6t,\ 8t,\ 12t . The values in the last two columns are the soltions at N c = 4 subscript 𝑁 𝑐 4 N_{c}=4 and in the TL extracted from Ref. Zheng  et al. , 2017a and LeBlanc  et al. , 2015 .

	$N_c=2$		$N_c=4$		$N_c=6$		$N_c=4$ Ref. Zheng et al.,2017a	TL Ref. LeBlanc et al.,2015
Method	DMET	RISB	DMET	RISB	DMET	RISB	DMET	DMET
$U/t=2$	-1.1804	-1.1673	-1.1790	-1.1693	-1.1790	-1.1704	-1.179	-1.1764
$U/t=4$	-0.8681	-0.8428	-0.8654	-0.8459	-0.8658	-0.8472	-0.863	-0.8604
$U/t=6$	-0.6541	-0.6306	-0.6545	-0.6362	-0.6553	-0.6376	-0.652	-0.6562
$U/t=8$	-0.5115	-0.4942	-0.5155	-0.5023	-0.5157	-0.5100	-	-0.5234
$U/t=12$	-0.3497	-0.3400	-0.3566	-0.3487	-0.3563	-0.3565	-	-0.3685

Table 2. Table 2: Energy E / t 𝐸 𝑡 E/t for DMET and RISB in the PM phase of the 2D Hubbard model at n = 0.8 𝑛 0.8 n=0.8 with the nearest neighbor hopping for N c = 2 , 4 , 6 subscript 𝑁 𝑐 2 4 6 N_{c}=2,\ 4,\ 6 at U = 2 ​ t , 4 ​ t , 6 ​ t , 8 ​ t 𝑈 2 𝑡 4 𝑡 6 𝑡 8 𝑡 U=2t,\ 4t,\ 6t,\ 8t . The values in the last two columns are the solutions at N c = 4 subscript 𝑁 𝑐 4 N_{c}=4 and in the TL extracted from Ref. LeBlanc  et al. , 2015 .

	$N_c=2$		$N_c=4$		$N_c=6$		TL Ref. LeBlanc et al.,2015
Method	DMET	RISB	DMET	RISB	DMET	RISB	DMET
$U/t=2$	-1.312	-1.300	-1.309	-1.302	-1.310	-1.302	-1.306
$U/t=4$	-1.129	-1.083	-1.122	-1.086	-1.120	-1.091	-1.108
$U/t=6$	-1.015	-0.927	-1.002	-0.938	-1.002	-0.942	-0.977
$U/t=8$	-0.950	-0.823	-0.932	-0.838	-0.923	-0.846	-0.880

Table 3. Table 3: Double occupancy ⟨ n ↑ ​ n ↓ ⟩ delimited-⟨⟩ subscript 𝑛 ↑ subscript 𝑛 ↓ \langle n_{\uparrow}n_{\downarrow}\rangle for DMET and RISB in the AFM phase of the half-filled 2D Hubbard model with the nearest neighbor hopping for N c = 2 , 4 , 6 subscript 𝑁 𝑐 2 4 6 N_{c}=2,\ 4,\ 6 at U = 2 ​ t , 4 ​ t , 6 ​ t , 8 ​ t , 12 ​ t 𝑈 2 𝑡 4 𝑡 6 𝑡 8 𝑡 12 𝑡 U=2t,\ 4t,\ 6t,\ 8t,\ 12t . The values in the last column are the solutions in the TL extracted from Ref. LeBlanc  et al. , 2015 .

	$N_c=2$		$N_c=4$		$N_c=6$		TL Ref. LeBlanc et al.,2015
Method	DMET	RISB	DMET	RISB	DMET	RISB	DMET
$U/t=2$	0.1937	0.1942	0.1934	0.1953	0.1935	0.1950	0.1913
$U/t=4$	0.1281	0.1314	0.1274	0.1300	0.1277	0.1300	0.1261
$U/t=6$	0.0819	0.0841	0.0815	0.0829	0.0816	0.0830	0.0810
$U/t=8$	0.0538	0.0548	0.0538	0.0542	0.0539	0.0541	0.0540
$U/t=12$	0.0268	0.0269	0.0272	0.0270	0.0272	0.0270	0.0278

Table 4. Table 4: Staggered magnetic moment m 𝑚 m for DMET and RISB in the AFM phase of the half-filled 2D Hubbard model with the nearest neighbor hopping for N c = 2 , 4 , 6 subscript 𝑁 𝑐 2 4 6 N_{c}=2,\ 4,\ 6 at U = 2 ​ t , 4 ​ t , 6 ​ t 𝑈 2 𝑡 4 𝑡 6 𝑡 U=2t,\ 4t,\ 6t . The values in the last two columns are the solutions at N c = 4 subscript 𝑁 𝑐 4 N_{c}=4 and in the TL extracted from Ref. Zheng  et al. , 2017a .

	$N_c=2$		$N_c=4$		$N_c=6$		$N_c=4$ Ref. Zheng et al.,2017a	TL Ref. Zheng et al.,2017a
Method	DMET	RISB	DMET	RISB	DMET	RISB	DMET	DMET
$U/t=2$	0.161	0.158	0.155	0.147	0.151	0.143	0.152	0.115
$U/t=4$	0.304	0.293	0.298	0.289	0.296	0.288	0.299	0.226
$U/t=6$	0.382	0.376	0.368	0.368	0.367	0.365	0.372	0.275

Equations

$$H=t\sum _{\sigma ,⟨i,j⟩}{c}_{i\sigma }^{†}{c}_{j\sigma }+\sum _{i}U{n}_{i\uparrow }{n}_{i\downarrow },$$

$$H=\sum _{⟨ij⟩,\alpha ,\beta }{\tilde{t}}_{ij}^{\alpha \beta }{c}_{i\alpha }^{†}{c}_{j\beta }+\sum _i{H}_{\text{loc}}[\{c_{i\alpha },c_{i\alpha }^{†}\}],$$

$$\tilde{t}_{ij}^{\alpha\beta}=\Big{\{}\begin{array}[]{c}t_{ij}^{\alpha\beta}\ \ \ \text{if}\ i\neq j\\ \ \ \ \ 0\ \ \ \ \ \text{otherwise}\end{array}.$$

$${\tilde{\epsilon }}_{\mathbf{k}}^{\alpha \beta }=\sum _i{e}^{-i\mathbf{k}\cdot {\mathbf{r}}_i}{\tilde{t}}_{i0}^{\alpha \beta },$$

$$H=\sum _{\mathbf{k}\in \text{RBZ},\alpha ,\beta }{\tilde{\epsilon }}_{\mathbf{k}}^{\alpha \beta }{c}_{\mathbf{k}\alpha }^{†}{c}_{\mathbf{k}\beta }+\sum _{\tilde{i}}{H}_{\text{loc}}[\{c_{i\alpha },c_{i\alpha }^{†}\}],$$

$$H_{eff}=\sum_{\mathbf{k}\in\text{RBZ}}\Big{[}R_{a\alpha}\tilde{\varepsilon}_{\mathbf{k}}^{\alpha\beta}R^{\dagger}_{\beta b}+\lambda_{ab}\Big{]}f_{\mathbf{k}a}^{\dagger}f_{\mathbf{k}b},$$

$$\Sigma (\omega )=-\omega \frac{1-R^{†}R}{R^{†}R}+\frac{1}{R}\lambda \frac{1}{R^{†}}.$$

$$\Sigma (\omega )=\lambda ,$$

$$\displaystyle H_{\text{emb}}=H_{\text{loc}}\big{[}\{c_{\alpha}^{\dagger},c_{\alpha}\}\big{]}$$

$$\displaystyle\qquad+\sum_{\alpha a}\big{(}\mathcal{D}_{a\alpha}c_{\alpha}^{\dagger}f_{a}+\text{H.c.}\big{)}+\sum_{ab}\lambda_{ab}^{c}f_{b}f_{a}^{\dagger},$$

$$\langle O\rangle=\langle\Phi|\hat{O}\big{[}\{c_{\alpha}^{\dagger},c_{\alpha}\}\big{]}|\Phi\rangle\,.$$

$$\mathcal{L}[|\Phi ⟩,R,\lambda ,{\Delta }^p;E^c,\mathcal{D},{\lambda }^c]=$$

$$\displaystyle-\frac{1}{\beta}\frac{N_{c}}{\mathcal{N}}\sum_{\mathbf{k}\in\text{RBZ}}\sum_{i\omega_{n}}\text{Tr log}\big{[}i\omega_{n}\mathbf{1}-R_{a\alpha}\tilde{\varepsilon}_{\mathbf{k}}^{\alpha\beta}R^{\dagger}_{\beta b}-\lambda_{ab}\big{]}e^{i\omega_{n}0^{+}}+\sum_{i}\text{Tr}\Big{[}E^{c}(\langle\Phi|\Phi\rangle-1)+\langle\Phi|H_{\text{emb}}|\Phi\rangle\Big{]}$$

$$\displaystyle-\sum_{iab}\big{(}\lambda_{ab}+\lambda_{ab}^{c}\big{)}\Delta_{ab}^{p}-\sum_{ica\alpha}\big{(}\mathcal{D}_{a\alpha}R_{c\alpha}+c.c\big{)}\big{[}\Delta^{p}(1-\Delta^{p})\big{]}_{ca}^{1/2},$$

$$\displaystyle\Delta_{ab}^{p}=\frac{N_{c}}{\mathcal{N}}\sum_{\mathbf{k}\in\text{RBZ}}\big{[}f_{T}(R\tilde{\varepsilon}_{\mathbf{k}}R^{\dagger}+\lambda)\big{]}_{ba},$$

$$\displaystyle\big{[}\Delta^{p}(1-\Delta^{p})\big{]}^{1/2}_{ac}\mathcal{D}_{c\alpha}=\frac{N_{c}}{\mathcal{N}}\sum_{\mathbf{k}\in\text{RBZ}}\big{[}\tilde{\varepsilon}_{\mathbf{k}}R^{\dagger}f_{T}(R\tilde{\varepsilon}_{\mathbf{k}}R^{\dagger}+\lambda)\big{]}_{\alpha a},$$

$$\displaystyle\sum_{cb\alpha}\frac{\partial}{\partial d^{p}_{s}}\big{[}\Delta^{p}(1-\Delta^{p})\big{]}^{\frac{1}{2}}_{cb}\big{[}\mathcal{D}\big{]}_{b\alpha}\big{[}R\big{]}_{c\alpha}+\text{c.c}+\big{[}l+l^{c}\big{]}_{s}=0,$$

$$H_{\text{emb}}|\Phi ⟩=E^c|\Phi ⟩,$$

$$\displaystyle\Big{[}\mathcal{F}^{(1)}\Big{]}_{ab}\equiv\langle\Phi|f_{b}f_{a}^{\dagger}|\Phi\rangle-\Delta_{ab}^{p}=0,$$

$$\displaystyle\Big{[}\mathcal{F}^{(2)}\Big{]}_{\alpha a}\equiv\langle\Phi|c_{\alpha}^{\dagger}f_{a}|\Phi\rangle-R_{c\alpha}\big{[}\Delta^{p}(1-\Delta^{p})\big{]}_{ca}^{\frac{1}{2}}=0.$$

$${\Delta }^p$$

$${\lambda }^c$$

$$\lambda$$

$$R$$

$$\displaystyle E=\sum_{\mathbf{k}\in\text{RBZ}}\sum_{ab}\big{[}R\tilde{\varepsilon}_{\mathbf{k}}R^{\dagger}f_{T}(R\tilde{\varepsilon}_{\mathbf{k}}R^{\dagger}+\lambda)\big{]}_{ab}$$

$$\displaystyle\qquad\qquad+\sum_{i}\langle\Phi|H_{i,loc}\big{[}c_{i\alpha}^{\dagger},c_{i\alpha}\big{]}|\Phi\rangle,$$

$$E=\sum _i⟨\Phi |\sum _{\alpha a}(D_{\alpha a}{c}_{\alpha }^{†}{f}_a)+H_{i,loc}[\{c_{\alpha }^{†}{c}_{\alpha }\}]|\Phi ⟩.$$

$$\displaystyle\Delta_{ab}^{p}=\frac{N_{c}}{\mathcal{N}}\sum_{\mathbf{k}\in\text{RBZ}}\big{[}f_{T}(\tilde{\varepsilon}_{\mathbf{k}}+\lambda)\big{]}_{ba},$$

$$\displaystyle\big{[}\Delta^{p}(1-\Delta^{p})\big{]}^{1/2}_{ac}\mathcal{D}_{c\alpha}=\frac{N_{c}}{\mathcal{N}}\sum_{\mathbf{k}\in\text{RBZ}}\big{[}\tilde{\varepsilon}_{\mathbf{k}}f_{T}(\tilde{\varepsilon}_{\mathbf{k}}+\lambda)\big{]}_{\alpha a},$$

$$\displaystyle\sum_{cb}\frac{\partial}{\partial d^{p}_{s}}\big{[}\Delta^{p}(1-\Delta^{p})\big{]}^{\frac{1}{2}}_{cb}\big{[}\mathcal{D}\big{]}_{bc}+\text{c.c}+\big{[}l+l^{c}\big{]}_{s}=0,$$

$$H_{\text{emb}}|\Phi ⟩=E^c|\Phi ⟩,$$

$$\displaystyle\Big{[}\mathcal{F}^{(1)}\Big{]}_{ab}\equiv\langle\Phi|f_{b}f_{a}^{\dagger}|\Phi\rangle-\Delta_{ab}^{p},$$

$$\displaystyle\Big{[}\mathcal{F}^{(2)}\Big{]}_{\alpha a}\equiv\langle\Phi|c_{\alpha}^{\dagger}f_{a}|\Phi\rangle-\big{[}\Delta^{p}(1-\Delta^{p})\big{]}_{\alpha a}^{\frac{1}{2}},$$

$$\displaystyle\Big{[}\mathcal{F}^{(3)}\Big{]}_{\alpha\beta}\equiv\langle\Phi|c_{\alpha}^{\dagger}c_{\beta}|\Phi\rangle-\Delta_{\alpha\beta}^{p},$$

$$\displaystyle\lambda_{\text{min}}:=\underset{\lambda}{\text{argmin}}\ \big{(}\|\mathcal{F}^{(1)}\|_{\text{F}}+\|\mathcal{F}^{(2)}\|_{\text{F}}+\|\mathcal{F}^{(3)}\|_{\text{F}}\big{)}\,,$$

Full text

I Introduction

Strongly correlated electron systems are still a most challenging problem in condensed-matter physics. In this area, quantum embedding approaches have proven to be invaluable tools for studying their electronic structure. In particular, dynamical mean-field theory (DMFT) Georges et al. (1996), density matrix embedding theory (DMET) Knizia and Chan (2012) and their respective cluster extensions have been successfully applied to many interacting model Hamiltonians as well as to real materials Georges et al. (1996); Kotliar et al. (2006); Maier et al. (2005); Hettler et al. (1998, 1999); Lichtenstein and Katsnelson (2000); Kotliar et al. (2001); Rohringer et al. (2018); Knizia and Chan (2012, 2013); Wouters et al. (2016); Zheng and Chan (2016); Zheng et al. (2017a, b); LeBlanc et al. (2015); Motta et al. (2017). The common basic idea underlying these schemes is to map the fully interacting lattice to a self-consistently determined impurity problem, for which a fragment of the original lattice, termed cluster, is treated as a correlated impurity coupled to a self-consistently determined non-interacting bath. The accuracy can be systematically improved by increasing the reference cluster size towards the thermodynamic limit (TL) and the size of the Hilbert space representing the non-interacting bath.

Another important theoretical method widely used for studying strongly correlated electron systems is the rotationally-invariant slave-boson theory (RISB) Frésard and Wölfle (1992); Lechermann et al. (2007); Lanatà et al. (2017), which is equivalent to the multi-orbital Gutzwiller approximation at the mean-field level Kotliar and Ruckenstein (1986); Bünemann and Gebhard (2007); Lanatà et al. (2008) and generally provides predictions almost as accurate as DMFT Isidori and Capone (2009); Ferrero et al. (2008, 2009); Mazin et al. (2014); Lanatà et al. (2015, 2017); Piefke and Lechermann (2018); Behrmann and Lechermann (2015) (especially for the ground-state properties) while being much less computationally demanding. Even if the foundation of the RISB mean-field theory is based on seemingly distinct ideas, it turns out that also this framework can be viewed as a quantum-embedding theory. In fact, it has been recently shown Lanatà et al. (2015) that the RISB equations can be cast, similarly to DMET, in terms of ground-state calculations of auxiliary impurity systems named “embedding Hamiltonians”, whose non-interacting bath is determined self-consistently based on the variational principle. Subsequently, it has been also shown Ayral et al. (2017) that DMET can be formally recovered from the RISB equation derived in Ref. Lanatà et al., 2017 by setting to unity the variational parameters encoding the mass renormalization weights.

RISB and DMET are especially essential for the situations where the computational cost of DMFT becomes prohibitively large due to the exponentially growing Hilbert space and/or the sign problem in quantum Monte Carlo. This usually happens for the 5$f$ systems, where the crystal-field effects, spin-orbit-coupling interaction and lattice relaxation have to be taken into account simultaneously, and for the large-scale cluster simulations of the Hubbard model. Many challenging problems, such as the equations of state of elemental actinides and the phase diagram of the high $T_{c}$ superconductors, rely on such approximations to gain a qualitative or even quantitative understanding LeBlanc et al. (2015); Zheng et al. (2017b); Lanatà et al. (2015). Hence, it is of important interest to characterize the respective accuracy and performance of these two approaches.

Here we perform comparative RISB and DMET benchmark calculations on the 1D and 2D Hubbard model against the available exact solution and the DMET values extrapolated to the TL. Zheng et al. (2017a); LeBlanc et al. (2015) Our numerical results indicate that the accuracy and the performance of these two methods are very similar for all the quantities studied, e.g., the total energy and local observables. Small differences between the two methods are found only for small cluster sizes, where RISB provides slightly more accurate predictions for the local observables (such as occupancy, double occupancy and local moments) as well as for the metal-insulator transition in the 2D Hubbard model.

Finally, we derive an alternative numerical implementation of DMET featuring a modified RISB algorithm with mass renormalization weights set to unity Ayral et al. (2017), which provides us with a new line of interpretation and paves the way for developing new generalizations and synergistic combination of these approaches (e.g., to systems at finite temperature and/or with inter-site electron-electron interactions or electron-phonon interactions Lanatà et al. (2015); Wang et al. (2010); Sandri et al. (2013); Motta et al. (2017); Sandhoefer and Chan (2016); Reinhard et al. ). This implementation makes it also possible to pattern an interface between density functional theory (DFT) and DMET after previous DFT+RISB and DFT+DMFT works Lanatà et al. (2015); Kotliar et al. (2006).

The paper is organized as follows: The Hubbard model is introduced in Sec. II. The RISB and DMET formalism and algorithmic structure are outlined in Sec. III. In Section IV are presented our benchmark simulation of the Hubbard model in 1D and 2D. Finally, Sec. V is devoted to concluding remarks.

II Model

Let us consider the 1D and 2D Hubbard model with the nearest neighbor hopping,

[TABLE]

where $t$ is the hopping amplitude, $i$ and $j$ are the indices for the lattice sites, and the $\sigma$ is the spin label, and $U$ is the local Coulomb interaction. $c_{i\sigma}^{(\dagger)}$ is the annihilation (creation) operator for the electron at site $i$ and spin $\sigma$.

The cluster extensions of RISB and DMET are both implemented by tiling the original lattice with clusters of increasing size Maier et al. (2005). Thus, the degrees of freedom of the single-band Hubbard model belonging to each cluster are treated as a single impurity, i.e., as if they were elementary (orbital) degrees of freedom of a multi-orbital Hubbard Hamiltonian represented as follows:

[TABLE]

where the indices $i,\ j=1,...,\mathcal{N}/N_{c}$ denote the enlarged unit cell, $\mathcal{N}$ is the total number of atoms and $N_{c}$ is the number of atoms within each cluster and the labels $\alpha,\ \beta=1,...,2N_{c}$ indicate the cluster spin and atom degrees of freedom.

In order to utilize the RISB and DMET theory, it is useful to define the inter-cluster hopping matrix as follows:

[TABLE]

The terms corresponding to the intra-cluster hopping parameters $t_{i\alpha,i\beta}$ are included within the operator $H_{\text{loc}}[\{c_{i\alpha},c_{i\alpha}^{\dagger}\}]$, along with the chemical potential and the local Coulomb interaction.

In our calculations, the translational invariance is exploited only partially, i.e., we represent the hopping matrix defined as:

[TABLE]

where the momentum $\mathbf{k}$ belongs to the reduced Brillouin zone (RBZ) of the enlarged unit cell containing the cluster. The resulting Hamiltonian in the momentum space is represented as follows:

[TABLE]

where $H_{\text{loc}}[\{c_{i\alpha,}c_{i\alpha}^{\dagger}\}]$ contains all the local one- and two-body terms.

III Methods

As shown in Refs. Lanatà et al., 2015; Ayral et al., 2017; Knizia and Chan, 2012, the RISB and DMET ground-state solution of the Hubbard Hamiltonian [Eq. (5)] is obtained by solving recursively two auxiliary systems: (i) a non-interacting system termed “effective-medium” or “quasiparticle Hamiltonian” and (ii) an interacting embedding impurity problem called “embedding Hamiltonian.”

The structure of the effective-medium Hamiltonian is the following:

[TABLE]

where $\tilde{\varepsilon}_{\mathbf{k}}$ was defined in Eq. (4), $R$ and $\lambda$ are $2N_{c}\times 2N_{c}$ complex matrices (the factor 2 arises from the spin degrees of freedom) and $\lambda$ is Hermitian. As we are going to show in Sec. III.1, in RISB both $R$ and $\lambda$ are determined self-consistently Lanatà et al. (2017) and their converged entries are connected to the self-energy $\Sigma(\omega)$ as follows: Lechermann et al. (2007); Lanatà et al. (2017)

[TABLE]

On the other hand, in DMET only the entries of $\lambda$ (called $u$ in the DMET literature) can vary while $R=\mathbf{1}$, i.e., the self-energy consist exclusively of the part representing the on-site energy shifts: Knizia and Chan (2012)

[TABLE]

see Sec. III.1.

The embedding Hamiltonian describes a multi-orbital dimer molecule containing a correlated impurity $c_{\alpha}^{(\dagger)}$ and a non-correlated bath $f_{a}^{(\dagger)}$. It reads:

[TABLE]

where $H_{\text{loc}}$ is defined in Eq. (2), $\mathcal{D}$ and $\lambda^{c}$ are $2N_{c}\times 2N_{c}$ complex matrices and $\lambda^{c}$ is Hermitian. The entries of both matrices are determined self-consistently Knizia and Chan (2012); Lanatà et al. (2015, 2017); Ayral et al. (2017), see Secs. III.1 and III.2. After convergence, the reduced density matrix of the impurity degrees of freedom (which is formally obtained by tracing out the bath degrees of freedom) provides the local reduced density matrix of the original physical system. In other words, the expectation value of any local operator \hat{O}\big{[}\{c_{\alpha}^{\dagger},c_{\alpha}\}\big{]}, such as the double occupancy or the local stagger magnetic moment, can be calculated from the ground state wavefunction $|\Phi\rangle$ of $H_{\text{emb}}$ as follows:Lanatà et al. (2015)

[TABLE]

III.1 Rotationally invariant slave-boson mean-field theory

The RISB theory is, in principle, an exact reformulation of the Hubbard system constructed by introducing auxiliary “slave” bosons coupled to “quasiparticle” fermionic degrees of freedom. Lechermann et al. (2007); Lanatà et al. (2015, 2017) As shown in Ref. Lanatà et al., 2015, the RISB mean-field theory is entirely encoded in the following Lagrange function:

[TABLE]

where: $R$ and $\lambda$ are the renormalization coefficients of the quasiparticle Hamiltonian introduced in Eq. (6), $H_{\text{emb}}$, $\mathcal{D}$ and $\lambda^{c}$ are the parameters of the embedding Hamiltonian introduced in Eq. (9), $|\Phi\rangle$ is the ground state wavefunction of $H_{\text{emb}}$, $E^{c}$ is a Lagrange multiplier enforcing the normalization of $|\Phi\rangle$ and $\Delta^{p}$ is the local density matrix of $H_{\text{eff}}$ (see Eq. (12)).

The self-consistency conditions determining the parameters of $H_{\text{emb}}$ and $H_{\text{eff}}$, see Eqs. (6) and  (9), are obtained by extremizing the mean-field Lagrange function with respect to $|\Phi\rangle,\ R,\ \lambda,\ \Delta^{p},\ E_{c},\ \mathcal{D},\ \text{and\ }\lambda^{c}$, which leads to the following equations:

[TABLE]

where the symbol $f_{T}$ stands for the Fermi function of a single-particle matrix at temperature $T$ and we utilized the following matrix parameterizations:

[TABLE]

where the set of matrices $h_{s}$ are an orthonormal basis of the space of Hermitian matrices (with respect to the canonical trace inner product). The parameters $d^{p}_{s}$, $l^{c}_{s}$ and $l_{s}$ are real, while $r_{s}$ is complex. The RISB saddle-point equations can be solved as follows:

Starting with an initial guess of $R$ and $\lambda$, compute $\Delta^{p}$ from Eq. (12). 2. 2.

From $\Delta^{p}$, calculate $\mathcal{D}$ from Eq. (13). 3. 3.

With $\mathcal{D}$ and $\Delta^{p}$, compute $\lambda^{c}$ from Eq. (14). 4. 4.

From $\mathcal{D}$ and $\lambda^{c}$, construct $H_{\text{emb}}$ from Eq. (9) and calculate its ground state $|\Phi\rangle$. 5. 5.

From $|\Phi\rangle$ and $\Delta^{p}$, calculate Eqs. (16) and (17) and utilize quasi-Newton methods to estimate the new $R$ and $\lambda$. 6. 6.

The convergence is achieved if Eqs. (16) and  (17) are satisfied. Otherwise, continue the root searching with the new $R$ and $\lambda$.

This structure is summarized schematically in Fig. 1.

Note that the Lagrange function [Eq. 11] evaluated for the converged parameters reduces to:

[TABLE]

which is the total energy of the system. Lanatà et al. (2017) It can be straightforwardly verified that, as long as Eqs. (12)-(17) are satisfied, the total energy can be equivalently expressed also as follows:

[TABLE]

III.2 Density matrix embedding theory

The self-consistency conditions determining the parameters of $H_{\text{emb}}$ and $H_{\text{eff}}$ in DMET can be formulated as follows:Ayral et al. (2017)

[TABLE]

where the symbol $\|...\|_{\text{F}}$ in Eq. 31 indicates the Frobenius norm. Note that Eqs. (24)-(29) are equivalent to Eqs. (12)-(17) with $R=\mathbf{1}$ and the constraint Eq. (30) was originally considered also in the Gutzwiller approximation (equivalent to RISB), but later was found to be unnecessaryFabrizio (2007).

The DMET equations can be solved as follows, see Fig. 1:

Starting with an initial guess of $\lambda$, calculate $\Delta^{p}$ using Eq. (24). 2. 2.

Compute $\mathcal{D}$ and $\lambda_{c}$ from Eq. (25) and Eq. (26) and construct the $H_{\text{emb}}$. 3. 3.

Compute the ground state $|\Phi\rangle$ and the corresponding single-particle density matrix, i.e.: $\langle\Phi|f_{b}f_{a}^{\dagger}|\Phi\rangle$, $\langle\Phi|c_{\alpha}^{\dagger}f_{a}|\Phi\rangle$ and $\langle\Phi|c_{\alpha}^{\dagger}c_{\beta}|\Phi\rangle$. 4. 4.

From $\langle\Phi|f_{b}f_{a}^{\dagger}|\Phi\rangle$, $\langle\Phi|c_{\alpha}^{\dagger}f_{a}|\Phi\rangle$ and $\langle\Phi|c_{\alpha}^{\dagger}c_{\beta}|\Phi\rangle$, determine the entries of $\lambda_{\text{min}}$ that minimize Eq. 31Zheng (note that such a minimum is generally larger than zero in interacting systems Knizia and Chan (2012); Ayral et al. (2017)). 5. 5.

Iterate until $\lambda_{\text{min}}$ is converged.

A quasi-Newton method Pulay (1980) is usually utilized to accelerate the convergence of DMET iteration. Once convergence is reached, the DMET total energy is computed from Eq. (23).Knizia and Chan (2012)

IV Results

Here, we benchmark RISB and DMET with cluster sizes $N_{c}=1,\ 2,\ 4,\ 6$ on the Hubbard model with the nearest neighbor hopping in 1D and 2D (on a square lattice). The DMET calculations below are all performed utilizing the implementation outlined in Sec. III.2, featuring a modified RISB algorithm with mass renormalization weights set to unity. Our results are compared to the DMET data obtained in Refs. Zheng et al., 2017a and LeBlanc et al., 2015.

IV.1 1D Hubbard model

In Fig. 2 the DMET and RISB behaviors of the energies as a function of the occupation $n$ for $U\ =\ 1t,\ 4t,\ 8t$ with $N_{c}\ =\ 1,\ 2,\ 4$ are shown in comparison with the exact Bethe Ansatz (BA)Lieb and Wu (1968) solutions. Overall, the DMET and RISB approximations to the total energies are very similar for all cluster sizes, and both techniques reproduce the BA results with less than 2% error already for $N_{c}=4$. The only difference observed is that the DMET energies are slightly more accurate at half-filling, while the RISB energies are more accurate away from half-filling.

In Figure 3 are shown the behaviors of the DMET and RISB occupancies $n$ as a function of the chemical potential $\mu$ for $U=\ 1t,\ 4t,\ 8t$ with $N_{c}=\ 1,\ 2,\ 4$, in comparison with the BA. The Mott insulating phase is characterized by a constant $n$ with compressibility $\frac{dn}{d\mu}=0$. At the Mott insulator-metal transition point $\mu_{c}$ the compressibility $\frac{dn}{d\mu}$ divergesCapone et al. (2004). In the metallic phase, $n$ decreases monotonically by decreasing $\mu$. We observe that both DMET and RISB capture the correct behavior for $N_{c}\geq 2$. Moreover, RISB yields more accurate $n$ and $\mu_{c}$ at $N_{c}=2$. However, at $N_{c}=4$ both DMET and RISB predicts very precise occupancy and $\mu_{c}$ with less than 5% error.

In Fig. 4 are shown the behaviors of the DMET and RISB double occupancies $\langle n_{\uparrow}n_{\downarrow}\rangle$ with $N_{c}\ =\ 1,\ 2,\ 4$, in comparison with the BA. At $N_{c}=1$ the DMET solutions are always metallic for every $U$; consequently, the double occupancy deviates from the BA results at large $U$. On the other hand in RISB, the double occupancy vanishes at the critical point $U_{c}\sim 10t$, i.e., the charge fluctuations are not captured in the Mott phaseBrinkman and Rice (1970). For $N_{c}=2$ both methods predict behaviors of $\langle n_{\uparrow}n_{\downarrow}\rangle$ that closely follow the BA values, although RISB is slightly more accurate. At $N_{c}=4$, both methods are very accurate with less than 7% error compared to BA.

We also analyze the convergence of the energy as a function of cluster size at filling $n=1$ and $n=0.75$ with $U=4t$ and $U=8t$ for DMET and RISB as shown in Fig. 5. DMET gives a better estimation for the ground-state energy at half-filling, while RISB yields more accurate energies at $n=0.75$. However, as the cluster size grows, both methods converge to the BA value rapidly. Our results are consistent with the data extracted from Ref. Zheng et al., 2017a, where an antiferromagnetic ground state was assumed (in 1D the ground state is non-magnetic).

IV.2 2D Hubbard model

Here we investigate the behaviors of the RISB and DMET solutions of the 2D Hubbard model on a square lattice with cluster sizes $N_{c}=1,\ 2,\ 4,\ 6$, see Fig. 6. These geometries are chosen so that the antiferromagnetic (AFM) ground state can be reproduced for $N_{c}\geq 2$ and that the paramagnetic (PM) and the AFM energetics can be compared on the same footing.

In Fig. 7 are shown the behaviors of the DMET and RISB total energy $E$ as a function of the Hubbard interaction $U$ at half-filling $n=1$ in the PM metal, PM insulating and AFM insulating phase, with cluster sizes $N_{c}=1,\ 2,\ 4$.

At $N_{c}=1$, DMET does not capture the Mott metal-insulator transition (MIT), i.e., it predicts a metallic solution for every value of U. On the other hand, RISB predicts a MIT at $U_{c}=12.6t$, where the total energy vanishes Brinkman and Rice (1970). For $N_{c}\geq 2$, both methods capture a MIT, as indicated by the crossing of the PM metal and PM insulator energies. Moreover, the energies of the AFM solutions are lower than the PM solutions, consistently with previous studies Knizia and Chan (2012).

It is also interesting to see how $U_{c}$ varies with the cluster size. We observe that in DMET $U_{c}$ is almost independent of the cluster size, e.g., $U_{c}=8.95t$ for $N_{c}=2$ and $U_{c}=9.65t$ for $N_{c}=4$. On the other hand, in RISB $U_{c}$ decreases from $12.6t$ for $N_{c}=1$ to $6.4t$ for $N_{c}=4$ (which is very close to the CDMFT value $U_{c}=6.05t$ for the same cluster size Park et al. (2008)).

Figure 8 shows the DMET and RISB occupancy $n$ as a function of chemical potential $\mu$ at $U=12t$ with $N_{c}\ =\ 2,\ 4$. We observe that in DMET the difference in the occupancy and the $\mu_{c}$ between $N_{c}=2$ and $N_{c}=4$ is large, while in RISB, the discrepancy between the two cluster sizes is small (less than 3% error). We conclude that RISB provides a slightly better description of the PM solutions.

The ground-state energy predicted from DMET and RISB are shown in Tabs. 1 and 2 for $n=1$ AFM phase and $n=0.8$ PM phase, respectively, with various $U$ and $N_{c}$. Our numerical values are compared to the DMET results at $N_{c}=4$ and in the TL in Refs. LeBlanc et al., 2015 and Zheng et al., 2017a, which are also shown as black solid dots in Fig. 7 at $n=1$.

We observe that at half-filling $n=1$ DMET gives overall more accurate predictions to the ground-state energies in the AFM phase compared to the TL energies LeBlanc et al. (2015) (see Tab. 1 and Fig. 7). However, the discrepancy between the two methods is already small at $N_{c}=4$ (less than 3% error). Away from half-filling ($n=0.8$), the ground-state energies predicted by RISB and DMET are equally accurate compared to the energies in the TL LeBlanc et al. (2015). Our DMET results are consistent with previous studies Zheng et al. (2017a); LeBlanc et al. (2015).

The double occupancies $\langle n_{\uparrow}n_{\downarrow}\rangle$ at $n=1$ in the AFM phase with different $N_{c}$ and $U$ are shown in Tab. 3. DMET yields slightly more precise double occupancy at $N_{c}=2$ for smaller $U$ compared to the TL results LeBlanc et al. (2015). However, for $N_{c}=4$, both methods obtained very accurate double occupancy close to the TL (less than 3% error).

In Tab. 4 we present the prediction of the AFM magnetic moment $m$ for both methods with different cluster sizes $N_{c}$ and $U$. Overall, we found the DMET and RISB magnetic moment are very similar, with RISB slightly closer to the TLLeBlanc et al. (2015).

V Conclusions

We have performed comparative benchmark calculations of RISB and DMET on the 1D and 2D (square lattice) Hubbard model with cluster sizes ranging from $N_{c}=1$ to $6$. We found that the overall performances of the two methods are very similar. Small differences are observed only for small cluster sizes, where RISB generally predicts slightly more accurate Mott MIT critical points, magnetic moments, occupancies and double occupancies. The DMET ground-state energy is usually more accurate around half-filling, while the RISB ground-state energy is more precise away from half-filling.

Furthermore, we proposed an alternative implementation of DMET featuring a modified RISB algorithm with a unity mass renormalization matrix. This formalism paves the ways for many generalizations. For example, the DFT+RISB derived in Ref. Lanatà et al., 2015 can now be readily transposed to DFT+DMET. The non-equilibrium extensions of both methods are also available Schiró and Fabrizio (2010, 2011); Mazza and Georges (2017); Kretchmer and Chan (2018). A systematic way of improving the accuracy of RISB without breaking translational symmetry has been recently proposed by introducing auxiliary “ghost” degrees of freedom Lanatà et al. (2017), and similar ideas have been applied also within the DMET framework Fertitta and Booth (2018). Other possible directions may be to generalize DMET to finite-temperature Sandri et al. (2013); Lanatà et al. (2015); Mazza and Georges (2017) or extending RISB to systems with electron-phonon interactions or inter-site electron-electron interactions Sandhoefer and Chan (2016); Reinhard et al. ; Motta et al. (2017).

VI Acknowledgements

T.-H. L. thanks G. Booth and Q. Chen for useful discussions on the DMET algorithm. Y. Y. thanks for the supports from BNL CMS center. T.-H. L, T. A., and G. K. were supported by the Department of Energy under Grant No. DE-FG02-99ER45761. N. L. was supported by the VILLUM FONDEN via the Centre of Excellence for Dirac Materials (Grant No. 11744). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) funded by NSF under Grants No. TG-DMR170121.