  We develop a hierarchical equations of motion (HEOM) approach to compute
real-time single-particle correlation functions and thermodynamic properties of
the Holstein model at finite temperature. We exploit the conservation of the
total momentum of the system to formulate the momentum-space HEOM whose
dynamical variables explicitly keep track of momentum exchanges between the
electron and phonons. Our symmetry-adapted HEOM enable us to overcome the
numerical instabilities inherent to the commonly used real-space HEOM. The HEOM
method is then used to study the spectral function and thermodynamic quantities
of chains containing up to ten sites. The HEOM results compare favorably to
existing literature. To provide an independent assessment of the HEOM approach
and to gain insight into the importance of finite-size effects, we devise a
quantum Monte Carlo (QMC) procedure to evaluate finite-temperature
single-particle correlation functions in imaginary time and apply it to chains
containing up to twenty sites. QMC results reveal that finite-size effects are
quite weak, so that the results on 5 to 10-site chains, depending on the
parameter regime, are representative of larger systems. A detailed comparison
between the HEOM and QMC data place our HEOM method among reliable methods to
compute real-time finite-temperature correlation functions in parameter regimes
ranging from low- to high-temperature, and weak- to strong-coupling regime.
