  Within the tensor network framework, the (positive) thermal density operator
can be approximated by a double layer of infinite Projected Entangled Pair
Operator (iPEPO) coupled via ancilla degrees of freedom. To investigate the
thermal properties of the spin-1/2 Heisenberg model on the square lattice, we
introduce a family of fully spin-$SU(2)$ and lattice-$C_{4v}$ symmetric on-site
tensors (of bond dimensions $D=4$ or $D=7$) and a plaquette-based
Trotter-Suzuki decomposition of the imaginary-time evolution operator. A
variational optimization is performed on the plaquettes, using a full (for
$D=4$) or simple (for $D=7$) environment obtained from the single-site Corner
Transfer Matrix Renormalization Group fixed point. The method is benchmarked by
a comparison to quantum Monte Carlo in the thermodynamic limit. Although the
iPEPO spin correlation length starts to deviate from the exact exponential
growth for inverse-temperature $\beta \gtrsim 2$, the behavior of various
observables turns out to be quite accurate once plotted w.r.t the inverse
correlation length. We also find that a direct $T=0$ variational energy
optimization provides results in full agreement with the
$\beta\rightarrow\infty$ limit of finite-temperature data, hence validating the
imaginary-time evolution procedure. Extension of the method to frustrated
models is described and preliminary results are shown.
