  Tensor network methods have become a powerful class of tools to capture
strongly correlated matter, but methods to capture the experimentally
ubiquitous family of models at finite temperature beyond one spatial dimension
are largely lacking. We introduce a tensor network algorithm able to simulate
thermal states of two-dimensional quantum lattice systems in the thermodynamic
limit. The method develops instances of projected entangled pair states and
projected entangled pair operators for this purpose. It is the key feature of
this algorithm to resemble the cooling down of the system from an infinite
temperature state until it reaches the desired finite-temperature regime. As a
benchmark we study the finite-temperature phase transition of the Ising model
on an infinite square lattice, for which we obtain remarkable agreement with
the exact solution. We then turn to study the finite-temperature Bose-Hubbard
model in the limits of two (hard-core) and three bosonic modes per site. Our
technique can be used to support the experimental study of actual effectively
two-dimensional materials in the laboratory, as well as to benchmark optical
lattice quantum simulators with ultra-cold atoms.
