  Progress in describing thermodynamic phase transitions in quantum systems is
obtained by noticing that the Gibbs operator $e^{-\beta H}$ for a
two-dimensional (2D) lattice system with a Hamiltonian $H$ can be represented
by a three-dimensional tensor network, 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 $D$. The coarse-graining is performed by a tree tensor network of
isometries. The isometries are optimized variationally --- taking into account
full tensor environment --- to maximize the accuracy of the PEPO. The algorithm
is applied to the isotropic quantum compass model on an infinite square lattice
near a symmetry-breaking phase transition at finite temperature. From the
linear susceptibility in the symmetric phase and the order parameter in the
symmetry-broken phase the critical temperature is estimated at ${\cal
T}_c=0.0606(4)J$, where $J$ is the isotropic coupling constant between $S=1/2$
pseudospins.
