Imaginary-Time Matrix Product State Impurity Solver for Dynamical Mean-Field Theory

We present a new impurity solver for dynamical mean-field theory based on imaginary-time evolution of matrix product states. This converges the self-consistency loop on the imaginary-frequency axis and obtains real-frequency information in a final real-time evolution. Relative to computations on the real-frequency axis, required bath sizes are much smaller and no entanglement is generated, so much larger systems can be studied. The power of the method is demonstrated by solutions of a three-band model in the single- and two-site dynamical mean-field approximation. Technical issues are discussed, including details of the method, efficiency as compared to other matrix-product-state-based impurity solvers, bath construction and its relation to real-frequency computations and the analytic continuation problem of quantum Monte Carlo methods, the choice of basis in dynamical cluster approximation, and perspectives for off-diagonal hybridization functions.

Popular Summary

Dynamical mean-field theory obtains an approximation to the electron self-energy of a full (typically infinite) system in terms of the solution of a simpler “quantum impurity” problem that consists of a set of localized, interacting levels coupled to a bath. While the approach has been very powerful, its applicability to interesting electronic phenomena occurring in broad classes of complex oxides has been limited because methods do not exist to solve large (many-orbital) impurity models with realistic interactions. Here, we significantly relax this limit, and we present a method for solving much more general classes of multiorbital impurity models.

Our key advance is a representation of the impurity model Green’s function via matrix product state-based computations formulated on the imaginary-time axis. Imaginary-time evolution does not create entanglement and hence dramatically increases the range of situations that can be studied. This development establishes, for the first time, the density matrix renormalization group, which is traditionally a computational approach to one-dimensional systems (and indeed to ground states and selected excited states), as a flexible, low-cost approach to a much wider class of problems. These problems include a general and previously inaccessible class of quantum impurity models relevant to the dynamical mean-field solution of realistic materials problems. We use our method to solve the lowest and next-to-lowest level of cluster dynamical mean-field approximation to study intersite physics of vanadate Ruddlesden-Popper materials with realistic rotationally invariant interactions.

We expect that with modest further optimization, many problems will be within reach, including ab initio studies of superconductivity in iron pnictide materials.

Figure 1

(a) Imaginary-time correlation functions ${\tilde{G}}^{\gtrless }(\tau )$ defined in Eq. (4a) for an impurity model arising in the context of the two-site dynamical cluster approximation to the single-band Hubbard model on the square lattice with next-nearest neighbor hopping $t^{\prime }/t=0.3$, half-bandwidth $D=4t$, interaction $U=2.5D$, and band filling $n=0.96$ in the paramagnetic phase. See Ref. [52] for definition of the model and the meaning of the orbital (patch) quantum number $K=\pm$. The dashed line is obtained using linear prediction for times $\tau D\ge 200$. (b) Maximal bond dimension $m$ of time-evolved states. The MPS computation uses a Hamiltonian representation of the discrete (approximate) impurity model with $L_c=2$ correlated sites and $L_b=14$ bath sites. The hybridization function of the impurity model ${\mathrm{\Lambda }}^{\text{discr}}$ is fitted using Eq. (13) for ${\beta }_{\mathrm{eff}}=315/D$ and $\alpha =0$. For the ground state optimization, we enforce a maximal bond dimension of $m=300$. The Krylov time evolution uses a time step of $\mathrm{\Delta }t=0.1/D$ and allows for a maximal global truncation error of $10^{-4}$ at each time step, adjusting bond dimensions automatically. This leads to an immediate decay of $m$ at $\tau \simeq 0$ from $m=300$ down to $m\simeq 110$, as seen in (b). We use the global SU(2) symmetry of the Hamiltonian to reach these low values of the bond dimension.

Figure 2

Fit of the hybridization function in the two-site DCA problem studied in Figs. 1 and 3, but here for the case $U=0$. The minimization [Eq. (13)] is done using $\alpha =0$ and a frequency grid defined by ${\beta }_{\mathrm{eff}}=100/D$ and a cutoff frequency of ${\omega }_c=6D$. Evidently, the quality of the fit does not improve any more for $L_b/L_c\gtrsim 9$.

Figure 3

Real- and imaginary-frequency Green’s functions computed by converging the DMFT self-consistency equation [Eq. (11)] for the two-site dynamical cluster approximation to the single-band Hubbard model on the square lattice with next-nearest neighbor hopping $t^{\prime }/t=0.3$, half-bandwidth $D=4t$, interaction $U=2.5D$, and band filling $n=0.96$ in the paramagnetic phase (as in Fig. 1). See Ref. [52] for definition of the model and the meaning of the orbital (patch) quantum number $K=\pm$. (a)–(c) Electron spectral function $A_{+}(\omega )$ (a) and $A_{-}(\omega )$ (b) obtained by converging on imaginary-frequency axis (ITMPS) using number of bath sites and different broadenings as specified in the figure, and compared to unbroadened ($\eta =0$) real-frequency axis computation using $L_b/L_c=39$ bath sites per correlated site of Ref. [39]. (d) Converged Matsubara Green’s function for number of bath sites shown, compared to numerically exact quantum Monte Carlo result of Ref. [52], computed at $\beta =200/D$.

Figure 4

Density per orbital as function of chemical potential for three-band Hubbard-Kanamori model Eq. (14) using the semielliptic density of states Eq. (15) and $U=12t$, obtained from single-site DMFT approximation evaluated using imaginary-time MPS (crosses) and CTQMC data (circles, Fig. 1 of Ref. [61], inverse temperature $\beta =50/t$). In the DMRG computations the bath fitting is performed using ${\beta }_{\mathrm{eff}}=100/t$, ${\omega }_c=6t$, and $\alpha =1$, with three bath sites per correlated site ($L_b/L_c=3$). The maximal matrix dimensions is $m=300$ for the ground state calculation, exploiting the SU(2) symmetry, which leads to the high precision $⟨(H-E{)}^2⟩\simeq 10^{-14}$. For the time evolution, we compute ${\tilde{G}}_a^{\gtrless }(\tau )$ in Eq. (4a) in steps of $\mathrm{\Delta }\tau =0.1/t$ allowing for a global truncation error of $5\times 10^{-4}$ per step, up to imaginary time ${\tau }_{\max }=100/t$, and use linear prediction for higher times.

Figure 5

(a) Imaginary part of Matsubara axis self-energy $\mathrm{\Sigma }$ and (b) imaginary part of hybridization function $\mathrm{\Lambda }$ for densities shown obtained from converged ITMPS solution of single-site DMFT for three-band Hubbard-Kanamori model [Eq. (16)] for $U=8t$ and $J=U/6$. Crosses represent ITMPS data and black circles depict CTQMC data from Fig. 3 of Ref. [62], computed at inverse temperature $\beta =100/t$. We choose all parameters as described in the caption of Fig. 4; in particular, for the bath fitting [Eq. (13)], we use ${\beta }_{\mathrm{eff}}=100/t$, ${\omega }_c=6t$, and $\alpha =1$. Choosing $\alpha =1$ enforces agreement for low frequencies at the price of disagreement at high frequencies, which is observed in both (a) and (b). In (b), $\mathrm{\Lambda }$ denotes the hybridization function that is fitted with the hybridization function ${\mathrm{\Lambda }}^{\text{discr}}$ of the discrete impurity model.

Figure 6

Comparison of results obtained using imaginary-time MPS with $L_b/L_c=3$ for single-site ($1s$) and two-site ($2s$) DMFT approximations to the Hubbard-Kanamori model [Eq. (14)] on the two-dimensional square lattice with half-bandwidth $D=4t$, ${ϵ}_k^{ab}=-2t(\cos k_x+\cos k_y){\delta }^{ab}$, $U^{\prime }=U-2J$, $J=U/4$, and $n=3$ ($\mu =5U/2-5J$), that is, in the particle-hole symmetric case. (a) Spectral functions for broadening $\eta =0.2D$; (b) broadening $\eta =0.05D$. In (c), we show the imaginary-time evolution of ${\tilde{G}}^{>}(\tau )$ as defined in Eq. (4a), confirming by comparison to a calculation for a smaller bath $L_b/L_c=2$ that this quantity has been converged with respect to the bath size. The maximal bond dimension for the ground state search is $m=1000$.

Figure 7

Sketch of a typical cluster-bath Hamiltonian ($L_c=3$, $L_b=6$) when it is mapped to a one-dimensional chain. Dashed lines depict couplings that do not mix symmetry sectors, and solid lines depict couplings that mix symmetry sectors.

Figure 8

Single-site DMFT for three-band Hubbard-Kanamori model as studied in Fig. 5. Here for the case $n=1.77$ ($\mu =5.0$) and $L_b/L_c=3$. To obtain the solution for $n=1.79$ as shown in Fig. 5, we choose $\mu =5.1$ and start from the $n=1.77$ solution. (a) Convergence of Matsubara Green’s function in the DMFT loop, starting from the noninteracting solution. (b) Convergence of the density in the DMFT loop. (c) Convergence of the ground state energy per particle in the DMFT loop. (d) Computation time. An iteration on the Matsubara axis takes about 30 min. The final real-axis computation (iteration 31) is considerably more expensive, but can still be optimized.

Figure 9

Single-particle excitation $|{\psi }_0⟩$ of the ground state $|E_0⟩$ and the subspace $\mathcal{N}=\{|\psi ⟩|⟨\psi |H|{\psi }_0⟩\ne 0\}$ of the Hilbert space that is relevant for the computation of a single-particle spectral function of the form $⟨{\psi }_0|\delta (\omega -H)|{\psi }_0⟩$. The single-particle excitation is very lowly entangled, the subspace is more strongly entangled, but still in general more lowly entangled than the rest of the Hilbert space.