A low-temperature quantum Fokker-Planck equation that improves the numerical stability of the hierarchical equations of motion for the Brownian oscillator spectral density

We investigate the numerical stability of the hierarchical equations of motion (HEOM) method applied to systems with the Brownian oscillator (BO) and multimode BO (MBO) spectral densities. It is shown that, with the increase in the system-bath coupling strength, the standard HEOM may become unstable, and a simple increase in the truncation depth of the HEOM cannot remove the instability at long times. To solve this problem, we first show that the high-temperature approximation of the HEOM with the BO spectral density is equivalent to the celebrated quantum Fokker-Planck equation (QFPE). By starting from the HEOM, we then derive a new multidimensional phase space differential equation that generalizes the QFPE to arbitrary temperature. It is further shown that the numerical instability can be removed if the new low-temperature QFPE is expanded in a basis set different than the one that leads to the conventional HEOM. The matrix product state method is also employed to propagate the new equation based on the low-temperature QFPE and to resolve the numerical instability problem for an electron transfer model with the MBO spectral density presented in the recent literature.