  The hierarchical equations of motion theory for Drude dissipation is
optimized, with a convenient convergence criterion proposed in advance of
numerical propagations. The theoretical construction is on basis of a Pad\'{e}
spectrum decomposition that has been qualified to be the best sum-over-poles
scheme for quantum distribution function. The resulting hierarchical dynamics
under the {\em apriori} convergence criterion are exemplified with a benchmark
spin-boson system, and also the transient absorption and two-dimensional
spectroscopy of a model exciton dimer system.
