  Tree tensor network states (TTNS) decompose the system wavefunction to the
product of low-rank tensors based on the tree topology, serving as the
foundation of the multi-layer multi-configuration time-dependent Hartree
(ML-MCTDH) method. In this work, we present an algorithm that automatically
constructs the optimal and exact tree tensor network operators (TTNO) for any
sum-of-product symbolic quantum operator.The construction is based on the
minimum vertex cover of a bipartite graph. With the optimal TTNO, we simulate
open quantum systems such as spin relaxation dynamics in the spin-boson model
and charge transport in molecular junctions. In these simulations, the
environment is treated as discrete modes and its wavefunction is evolved on
equal footing with the system. We employ the Cole-Davidson spectral density to
model the glassy phonon environment, and incorporate temperature effects via
thermo field dynamics. Our results show that the computational cost scales
linearly with the number of discretized modes, demonstrating the efficiency of
our approach.
