Efficient implementation of the Gutzwiller variational method

We present a self-consistent numerical approach to solve the Gutzwiller variational problem for general multiband models with arbitrary on-site interaction. The proposed method generalizes and improves the procedure derived by Deng et al. [Phys. Rev. B 79, 075114 (2009)], overcoming the restriction to density-density interaction without increasing the complexity of the computational algorithm. Our approach drastically reduces the problem of the high-dimensional Gutzwiller minimization by mapping it to a minimization only in the variational density matrix, in the spirit of the Levy and Lieb formulation of density functional theory (DFT). For fixed density the Gutzwiller renormalization matrix is determined as a fixpoint of a proper functional, whose evaluation requires only ground-state calculations of matrices defined in the Gutzwiller variational space. Furthermore, the proposed method is able to account for the symmetries of the variational function in a controlled way, reducing the number of variational parameters. After a detailed description of the method we present calculations for multiband Hubbard models with full (rotationally invariant) Hund’s rule on-site interaction. Our analysis shows that the numerical algorithm is very efficient, stable, and easy to implement. For these reasons this method is particularly suitable for first-principles studies (e.g., in combination with DFT) of many complex real materials, where the full intra-atomic interaction is important to obtain correct results.

Figure 1

Flow chart representing the numerical calculation of the functional ${\mathcal{T}}_{n^0}$.

Figure 2

Comparison of results for two degenerate bands on the 3D cubic lattice with rotational-invariant Hunds interaction, from Gutzwiller (solid lines) and slave-boson[43] (dotted lines), showing the quasiparticle weight $\mathcal{Z}$ as a function of $U$ and (from right to left), $J/U=0$, $0.01$, $0.02$, $0.05$, $0.10$, $0.20$, $0.45$. Additionally, the higher energy fixpoint solutions of Eq. (88) at finite $J/U$ are reported.

Figure 3

Filling per spin of orbital 1 for $\Delta =0.2$ and different values of $J/U$. From bottom to top, $J/U=0,0.01,0.02,0.05,0.1,0.15,0.25$. The continue lines correspond to our Gutzwiller results for the metallic phase at half-filling. Open (full) symbols correspond to metallic (insulating) solutions obtained from DMFT[52] at inverse temperature $\beta =25$.

Figure 4

Gutzwiller expectation value of ${\widehat{S}}^2$ for different values of $J/U$ and crystal field splittings $\Delta$, the total filling is six electrons per site and $U=1$. Comparison with DMFT data at inverse temperature $\beta =25$ from Ref. 53.

Figure 5

Gutzwiller renormalization matrix $\mathcal{Z}$ and filling of the bonding-antibonding bands for the two-bands bilayer Hubbard model with equal bandwidths, local hybridization $V=0.25$, $U^{\prime }=J^{\prime }=J=0$, and filling $N=1.88$ (solid lines). Comparison with slave-boson results from Ref. 43 (dotted lines).

Figure 6

Convergence of the forward recursion scheme (circles) and of the Newton method (squares) for five bands at $N=5$, $J=0$, and $\Delta =0$. The quasiparticle weight $\mathcal{Z}\equiv {\mathcal{R}}_{\alpha \alpha }^2$ goes to zero as $U$ goes to the critical coupling $U_c\approx 10$. Simultaneously, the leading eigenvalue $\lambda$ of the Jacobian of the recursion function goes to 1, and the number of forward iterations required to reach a fixed relative precision (here $|{\mathcal{R}}_i-{\mathcal{R}}_{i+1}|\le {10}^{-6}$ is used) diverges. On the contrary, the number of Newton iterations is almost independent of $U$. The matrix ${\left[{\mathcal{R}}_0\right]}_{\alpha \beta }={\delta }_{\alpha \beta }$ is used as initial condition of both the forward recursion series and the Newton method $\forall U$.

Figure 7

Sweep in the filling of the antibonding orbital $n_{-}$ for the two-bands bilayer Hubbard model with equal bandwidths, local hybridization $V=0.25$, filling $N=1.88$, $U^{\prime }=J^{\prime }=J=0$, and $U=2.5$. Renormalization matrix $\mathcal{Z}$ (top panel) and total energy (bottom panel).