  The variational tensor network renormalization approach to two-dimensional
(2D) quantum systems at finite temperature is applied for the first time to a
model suffering the notorious quantum Monte Carlo sign problem --- the orbital
$e_g$ model with spatially highly anisotropic orbital interactions.
Coarse-graining of the tensor network along the inverse temperature $\beta$
yields a numerically tractable 2D tensor network representing the Gibbs state.
Its bond dimension $D$ --- limiting the amount of entanglement --- is a natural
refinement parameter. Increasing $D$ we obtain a converged order parameter and
its linear susceptibility close to the critical point. They confirm the
existence of finite order parameter below the critical temperature $T_c$,
provide a numerically exact estimate of~$T_c$, and give the critical exponents
within $1\%$ of the 2D Ising universality class.
