Variational tensor network renormalization in imaginary time : benchmark results in the Hubbard model at finite temperature

A Gibbs operator $e^{-\beta H}$ for a two-dimensional (2D) lattice system with a Hamiltonian H can be represented by a 3D tensor network, with the third dimension being the imaginary time (inverse temperature) $\beta$. Coarse graining the network along \beta results in a 2D projected entangled-pair operator (PEPO) with a finite bond dimension. The coarse graining is performed by a tree tensor network of isometries. They are optimized variationally to maximize the accuracy of the PEPO as a representation of the 2D thermal state $e^{-\beta H}$. The algorithm is applied to the two-dimensional Hubbard model on an infinite square lattice. Benchmark results at finite temperature are obtained that are consistent with the best cluster dynamical mean-field theory and power-series expansion in the regime of parameters where they yield mutually consistent results.