  We provide an efficient approximation for the exponential of a local operator
in quantum spin systems using tensor-network representations of a cluster
expansion. We benchmark this cluster tensor network operator (cluster TNO) for
one-dimensional systems, and show that the approximation works well for large
real- or imaginary-time steps. We use this formalism for representing the
thermal density operator of a two-dimensional quantum spin system at a certain
temperature as a single cluster TNO, which we can then contract by standard
contraction methods for two-dimensional tensor networks. We apply this approach
to the thermal phase transition of the transverse-field Ising model on the
square lattice, and we find through a scaling analysis that the cluster-TNO
approximation gives rise to a continuous phase transition in the correct
universality class; by increasing the order of the cluster expansion we find
good values of the critical point up to surprisingly low temperatures.
