  Numerical implementation of an explicit phonon bath requires a large number
of oscillator modes in order to maintain oscillators at the initial temperature
when modeling energy relaxation processes. An additional thermalization
algorithm may be useful in controlling the local temperature. In this paper we
extend our previously proposed thermalization algorithm [M. Jaku\v{c}ionis and
D. Abramavi\v{c}ius, Phys. Rev. A 103, 032202 (2021) ] to be used with the
numerically exact multiple-Davydov-D2 trial wave function for simulation of
relaxation dynamics and spectroscopic signals of open quantum systems using the
time-dependent Dirac-Frenkel variational principle. By applying it to the
molecular aggregate model, we demonstrate how the thermalization approach
significantly reduces the numerical cost of simulations by decreasing the
number of oscillators needed to explicitly simulate the aggregate's environment
fluctuations while maintaining correspondence to the exact population
relaxation dynamics. Additionally, we show how the thermalization can be used
to find the equilibrium state of the excited molecular aggregate, which is
necessary for simulation of the fluorescence and other spectroscopic signals.
The thermalization algorithm we present offers the possibility to investigate
larger system-bath models than was previously possible using the
multiple-Davydov-D2 trial wave function and local heating effects in molecular
complexes.
