Nonequilibrium diagrammatic many-body simulations with quantics tensor trains

The nonequilibrium Green's function formalism provides a versatile and powerful framework for numerical studies of nonequilibrium phenomena in correlated many-body systems. For calculations starting from an equilibrium initial state, a standard approach consists of discretizing the Kadanoff-Baym contour and implementing a causal time-stepping scheme in which the self-energy of the system plays the role of a memory kernel. This approach becomes computationally expensive at long times, because of the convolution integrals and the large amount of computer memory needed to store the Green's functions. A recent idea for the compression of nonequilibrium Green's functions is the quantics tensor train representation. Here, we explore this approach by implementing equilibrium and nonequilibrium simulations of the two-dimensional Hubbard model with a second-order weak-coupling approximation to the self-energy. We show that calculations with compressed two-time functions are possible without any loss of accuracy, and that the quantics tensor train implementation shows a much improved scaling of the computational effort and memory demand with the length of the time contour.

Figure 1

Second-order contribution to the self-energy in the real-space representation. $z$ and $z^{\prime }$ are time points on the KB contour $\mathcal{C}$.

Figure 2

Discretization of the three-legged KB contour $\mathcal{C}$ and weight factors for the trapezoidal rule integration.

Figure 3

Equilibrium unfolded Green's function $G_k(z,z^{\prime })$ for $U=2$, $\beta =5$, $N_t=N_{\tau }=800$, and $k=(\pi ,\pi )$. The top (bottom) panel shows the real (imaginary) part. In the top panel we also indicate the different Green's function components and symmetry relations with respect to different axes.

Figure 4

Illustration of the quantics tensor train representation of the function $f\left(z\right)$. First, the argument $z$ is expressed in a binary representation, which introduces the bits (or spins) $z_i$, $i=1,...,R$. Then, the function defined on the $2^R$ dimensional space is decomposed into a tensor train.

Figure 5

Bond dimensions of the Green's function (top) and self-energy (bottom) on a logarithmic scale, as well as the results for the individual components, for cutoff ${\epsilon }_{\text{cutoff}}={10}^{-15}$. The parameters are $k=(\pi ,\pi )$, $U=2$, $\beta =5$, $N_t=N_{\tau }=800$. The QTT has $R^{\prime }=2R=24$ bits ($R=12$ for each time variable) in case of the full function, 20 bits for the retarded, lesser, and left-mixing components, and 10 bits in case of the Matsubara component. The dashed line represents the worst case scenario without scale separation.

Figure 6

Maximum bond dimension as a function of $k$ in a quarter of the first BZ for the full Green's function $G_k$ (top left) and its lesser (top right), retarded (bottom left), and left-mixing (bottom right) components. $U=2$, $\beta =5$, $N_t=N_{\tau }=800$, ${\epsilon }_{\text{cutoff}}={10}^{-15}$.

Figure 7

Maximum bond dimension as a function of $k$ in a quarter of the first BZ for the full self-energy ${\mathrm{\Sigma }}_k$ (top left) and its lesser (top right), retarded (bottom left), and left-mixing (bottom right) components. $U=2$, $\beta =5$, $N_t=N_{\tau }=800$, ${\epsilon }_{\text{cutoff}}={10}^{-15}$.

Figure 8

Exponential convergence of the QTT representation with $R$. The figure plots SMAPE for the reference state $G_{\text{matrix}}$ (circles). An offset $C=1.6$ is added, because the exponential convergence is towards the infinite-resolution function $G_{\infty }$.

Figure 9

(Top) Maximum norm error for $G_{\vec{k}}^{\left(l\right)}-G_{\vec{k}}^{(l+1)}$ as a function of iterations $l$ for $U=2$ and 4, $\beta =2$, $(k_x,k_y)=(3.14,3.14)$ and $(1.57,1.26)$ (near the Fermi surface). The circles show the results of the QTT implementation and the lines indicate the reference data from the matrix implementation. (Bottom) SMAPE for $G_{\text{QTT}}$ and $G_{\text{matrix}}$ as a function of iterations $l$, for the same parameters. We use ${\epsilon }_{\text{cutoff}}={10}^{-15}$ for the left panels and ${\epsilon }_{\text{cutoff}}={10}^{-8}$ for the right ones.

Figure 10

Real (left) and imaginary (right) noninteracting (top) and interacting (middle: $U=2$, bottom: $U=4$) Green's functions for $\beta =2$ and $k=(\pi ,\pi )$. We restrict the color bar to the interval $[-0.01,0.01]$ to emphasize small structures.

Figure 11

CPU and memory demand as a function of increasing mesh size for fixed $\beta =2$, $t_{\text{max}}=2$ (iteration $l=7$ for the CPU measurement). On the horizontal axis, we report the number of digits $R$ in the binary representation, which corresponds to $2^R$ time discretization points on the contour $\mathcal{C}$. Fitting the results for the matrix implementation to ${\left(2^R\right)}^b$ yields the exponent $b=1.70\pm 0.03(1.57\pm 0.04)$ for the CPU (RAM) scaling. In the QTT calculation, we set the maximum bond dimension to $D=100$ and the cutoff to ${\epsilon }_{\text{cutoff}}={10}^{-8}$.

Figure 12

(Top) Maximum bond dimension of the converged Green's function ($l=7$) with $k=(\pi ,\pi )$ as a function of increasingly fine grid, corresponding to the setup of Fig. 11. (Bottom) Maximum bond dimension of the converged Green's function ($l=7$) with $k=(\pi ,\pi )$ as a function of maximum time $t_{\text{max}}$, corresponding to the setup of Fig. 13.

Figure 13

CPU and memory demand as a function of increasing number of digits $R$ in the binary representation. The corresponding values of $N_t$ are indicated in the brackets on the horizontal axis ($t_{\text{max}}=N_{t}dt$). Here, $\beta =1$, $d\tau =0.009$, and $N_{\tau }=108$ are fixed. Fitting to ${\left(2^R\right)}^b$ yields the exponent $b=2.098\pm 0.003(1.78\pm 0.06)$ for the CPU (RAM) scaling. In the QTT calculation, we set the maximum bond dimension to $D=100$ and the cutoff to ${\epsilon }_{\text{cutoff}}={10}^{-8}$.

Figure 14

(Top) Maximum norm error for $G_k^{\left(l\right)}-G_k^{(l+1)}$ as a function of iterations $l$ for interaction ramps $U\left(t\right)$ to the indicated values of $U_{\text{final}}$, $\beta =2$ in the initial state, $(k_x,k_y)=(3.14,3.14)$ and $(1.57,1.26)$ (near the Fermi surface). The circles are the result of the QTT implementation and the lines indicate the reference data from the matrix implementation. (Bottom) SMAPE for $G_{\text{QTT}}$ and $G_{\text{matrix}}$ as a function of iterations $l$, for the same parameters.

Figure 15

Real (top) and imaginary (bottom) part of the converged Green's function for the interaction ramp $U\left(t\right)$ to $U_{\text{final}}=4$, $\beta =2$ in the initial state, and $k=(\pi ,\pi )$.

Figure 16

Kinetic energy of the lattice system subject to an interaction ramp $U\left(t\right)$ to $U_{\text{final}}=2$ and 4. The lines (circles) show the results from the QTT (matrix) implementation. The initial inverse temperature is $\beta =2$, $\kappa =0.5$ and $t_{\text{max}}=2$.