Solving nonequilibrium dynamical mean-field theory using matrix product states

We solve nonequilibrium dynamical mean-field theory (DMFT) using matrix product states (MPS). This allows us to treat much larger bath sizes and by that reach substantially longer times (factor $\sim 2–3$) than with exact diagonalization. We show that the star geometry of the underlying impurity problem can have substantially better entanglement properties than the previously favored chain geometry. This has immense consequences for the efficiency of an MPS-based description of general impurity problems: in the case of equilibrium DMFT, it leads to an orders-of-magnitude speedup. We introduce an approximation for the two-time hybridization function that uses time-translational invariance, which can be observed after a certain relaxation time after a quench to a time-independent Hamiltonian.

Figure 1

(a) Semielliptic bath spectral function $-\frac{1}{\pi }\mathrm{Im}\Lambda \left(\omega \right)$ defined in Eq. (20). (b) Corresponding couplings in the chain (${\tilde{V}}_l$) and the star ($V_l$) geometry. (c) Corresponding potentials in the chain (${\tilde{\epsilon }}_l$) and the star (${\epsilon }_l$) geometry. Parameters are given in units of the hopping $v$.

Figure 2

Sketch of the three setups studied. The star geometry can be mapped with the unitary transform $U$ to the chain geometry (i). It can also be mapped to an auxiliary chain by sorting the indices ascendingly to their potential energy. If one places the impurity in the center of this chain, one obtains the layout (ii), if one places it on the left edge of the chain, layout (iii) is obtained. Layouts (ii) and (iii) differ by the range over which the couplings $V_l$ couple different lattice sites.

Figure 3

(a)–(c) Properties of the MPS approximation $|E_0\rangle$ of the ground state of a SIAM in the chain (full blue lines) in the three setups sketched in Fig. 2: (i) chain geometry (full blue lines), (ii) and (iii) star geometries with central and edge impurity (dashed green lines and red dotted lines). This is for the semielliptic bath spectral function (20) shown in Fig. 1 and for $U/v=4$. (d), (e) Properties of the initial state $c_{0\sigma }^{†}|E_0\rangle$ for the time evolution needed to compute the retarded Green's function. (a) Density distribution $n_l$. (b), (d) Bond dimension $m$. (c), (e) Bond entanglement entropy $S_b$. $L_b=39$ sites are used to approximate the bath. The total chain length is $L=L_b+1=40$. Ground states have been computed with a maximum bond dimension of $m=500$. In the case of the chain geometry (i) this sufficed to reach a variance of $\langle {[(H^{\mathrm{chain}}-E_0)/v]}^2\rangle \sim {10}^{-4}$, whereas in the case of the star geometry, (ii) and (iii), one could reach $\langle {[(H^{\mathrm{star}}-E_0)/v]}^2\rangle \sim {10}^{-6}$. Here, $E_0$ denotes the numerical value of the ground-state energy.

Figure 4

Time evolution of a SIAM with semielliptic bath spectral function (20) at $U/v=4$ for the three setups (i)–(iii) shown in Fig. 2. The properties of the initial states for the evolution are shown in Fig. 3. Panels (a), (d), and (g) refer to the chain geometry (i), panels (b), (e), and (h) refer to the star geometry with the impurity located at the center (ii), and panels (c) and (f) to the star geometry with the impurity located at the left edge (iii). Panels (a), (b), and (c) show the local bond dimension $m$ plotted versus bond and time, panels (d), (e), and (f) show the bond entanglement entropy $S_b$, and panels (g) and (h) show the density distribution substracted from its initial value $n_l\left(t\right)-n_l\left(0\right)$ plotted versus site $l$ and time.

Figure 5

Panels (a)–(c) refer to the computation with semiellipitic bath spectral function [Fig. 1]. Panels (d)–(h) refer to a computation with the self-consistently determined bath spectral function (d) for $U/v=4$. Panels (a) and (g) show computation time versus physical time. Therein, the blue solid line refers to the chain geometry (i) whereas the green [red] dashed [dotted] line refers to the star geometry with the impurity located at the center (ii) [at the edge (iii)]. Panels (b) and (h) show the time evolution of the greater Green's function and panels (c) and (i) show the Fourier transform of the Green's function, which is the same for all three setups (i)–(iii), and therefore only one curve is shown. The oscillations in the resolution of $-\mathrm{Im}G\left(\omega \right)$ in panel (c) can be removed by convolution with a Gaussian or a Lorentzian of small width $\eta$. On the (real-) time domain, this would correspond to a slight damping (windowing) of $G^{>}\left(t\right)$ with a Gaussian or Lorentzian of large width $1/\eta$ and maximum at $t=0$. This suppresses contributions for times $t\gtrsim 1/\eta$. Alternatively, one can compute the real-time evolution of the Green's function up to higher times, until it has converged to zero, or use an extrapolation technique such as linear prediction [32, 39, 40, 41].

Figure 6

Time evolution for SIAMs in both geometries for a SIAM with self-consistently determined bath spectral function as shown in Fig. 5. The definition of panels is analogous to the one of Fig. 4.

Figure 7

Time evolution of double occupancy $\langle d\left(t\right)\rangle$ (a) and kinetic and total energy (b) for $U=10$ and different bath sizes. For the largest bath shown $L_b=24$ we could reach a time $t_{\max }\sim 7/v_0$ in a controlled way, meaning that the total energy is conserved. We did not limit the maximal allowed bond dimension $m$, but we imposed an upper error bound for the time evolution in single time step of $\Delta t=0.05/v_0$ of ${\epsilon }_{\mathrm{err}}={10}^{-6}$, as defined in Eq. (23).

Figure 8

Solid lines refer to the $U=10$ case, dashed lines to the $U=4$ case. Panel (a) shows the error needed to compute one time step in the DMFT scheme of $\Delta t=0.05/v_0$. Panel (b) shows the computation time $max|\Lambda (t,t^{\prime })-{\Lambda }_{\mathrm{Cholesky}}(t,t^{\prime })|$ of the hybridization function. Panel (c) shows the maximal bond dimension that occurs in set of Krylov states needed to expand $exp[-iH(t\left)\Delta t\right]$. Panel (d) shows the average bond dimension in these states. The average dimension is much lower than the maximal dimension, which can be understood when looking at the spatially resolved bond dimensions that are shown in Figs. 4 and 6 for the equilibrium case, but are typical also for the nonequilibrium case.

Figure 9

Time evolution of double occupancy $\langle d\left(t\right)\rangle$ (a) and kinetic and total energy (b) for $U=4$ and different bath sizes. For the largest bath shown $L_b=18$ we could reach a time $t_{\max }\sim 5.5/v_0$ in a controlled way, meaning that the total energy is conserved.

Figure 10

(a) Hybridization function $\Lambda (t,t^{\prime })$ obtained from a calculation in which the self-consistency has been computed on the “full” $t-t^{\prime }$ grid. (b) Difference of hybridization functions ${\Lambda }_{\mathrm{relax}}(t,t^{\prime })-\Lambda (t,t^{\prime })$, where for ${\Lambda }_{\mathrm{relax}}(t,t^{\prime })$, the self-consistency has only been solved for times $t^{\prime }<t_{\mathrm{relax}}=1/v_0$.

Figure 11

Comparison between calculation based on the approximated $\Lambda (t,t^{\prime })$, for which self-consistency has only been computed for $min(t,t^{\prime })<t_{\mathrm{relax}}$, and the exact calculation.

Figure 12

Bond dimension $m$ versus time and bond for different interactions strengths $U$ as given in the figure. We compare the chain geometry (i) (panels (a), (c), (e), (g)) and the star geometry with the impurity located at the center (ii) (panels (b), (d), (f), (h)). The geometries are defined in Fig. 2.

Figure 13

Results for the solution of the impurity model with the hybridization function defined in Eq. (B1) for two interactions $U=10$ (a) and $U=4$ (b). This does not involve the solution of a full self-consistent computation.

Figure 14

Time evolution of couplings for the case with homogeneous potential ${\epsilon }_{l\sigma }=0$ (full lines) and the case with ${\epsilon }_{l\sigma }=-2v_0+l4v_0/L_b$, where $l$ runs over $L_b$ bath sites.