  The hierarchical equations of motion (HEOM) method is a numerically exact
open quantum system dynamics approach. The method is rooted in an exponential
expansion of the bath correlation function, which in essence strategically
reshapes a continuous environment into a set of effective bath modes that allow
for more efficient cutoff at finite temperatures. Based on this understanding,
one can map the HEOM method into a Schr\"odinger-like equation with a
non-Hermitian super Hamiltonian for an extended wavefunction being the tensor
product of the central system wave function and the Fock state of these
effective bath modes. Recognizing that the system and these effective bath
modes form a star-shaped entanglement structure, in this work, we explore the
possibility of representing the extended wave function as an efficient tree
tensor network state (TTNS), the super Hamiltonian as a tree tensor network
operator of the same structure, as well as the application of a time
propagation algorithm using the time-dependent variational principle. Our
benchmark calculations based on the spin-boson model with a slow-relaxing bath
show that, the proposed HEOM+TTNS approach yields consistent results with that
of the conventional HEOM method, while the computation is considerably sped up
by a factor of a few orders of magnitude. Besides, the simulation with a
genuine TTNS is four times faster than a one-dimensional matrix product state
decomposition scheme.
