The partition function (PF) plays a key role in the calculation of quantum thermodynamic properties of a system that interacts with a heat bath. The imaginary-time hierarchical equations of motion (imHEOM) approach was developed to evaluate in a rigorous manner the PF of a system strongly coupled to a non-Markovian bath. In this paper, we present a numerically efficient scheme to evaluate the imHEOM utilizing the β-differentiated imHEOM (BD-imHEOM) that are obtained by differentiating the elements of the imHEOM with respect to the inverse temperature. This approach allows us to evaluate the system, system-bath interaction, and heat-bath parts of the PF efficiently. Moreover, we employ a polyharmonic decomposition method to construct a concise hierarchical structure with better convergence, thus reducing the cost of numerical integrations. We demonstrate the proposed approach by compute thermodynamic quantities of a spin-boson system and a 2 × 2 antiferromagnetic triangular spin lattice system with an Ohmic spectral distribution.
The partition function (PF) plays a key role in the calculation of quantum thermodynamic properties of a system that interacts with a heat bath. The imaginary-time hierarchical equations of motion (imHEOM) approach was developed to evaluate in a rigorous manner the PF of a system strongly coupled to a non-Markovian bath. In this paper, we present a numerically efficient scheme to evaluate the imHEOM utilizing the β-differentiated imHEOM (BD-imHEOM) that are obtained by differentiating the elements of the imHEOM with respect to the inverse temperature. This approach allows us to evaluate the system, system-bath interaction, and heat-bath parts of the PF efficiently. Moreover, we employ a polyharmonic decomposition method to construct a concise hierarchical structure with better convergence, thus reducing the cost of numerical integrations. We demonstrate the proposed approach by compute thermodynamic quantities of a spin-boson system and a 2 × 2 antiferromagnetic triangular spin lattice system with an Ohmic spectral distribution.
For a system strongly coupled to a heat bath, the quantum coherence of the system and the heat bath plays an important role in the system dynamics. This is particularly true in the case of non-Markovian noise. We rigorously investigate the influence of system-bath coherence by deriving the reduced hierarchal equations of motion (HEOM), not only in real time, but also in imaginary time, which represents an inverse temperature. It is shown that the HEOM in real time obtained when we include the system-bath coherence of the initial thermal equilibrium state possess the same form as those obtained from a factorized initial state. We find that the difference in behavior of systems treated in these two manners results from the difference in initial conditions of the HEOM elements, which are defined in path integral form. We also derive HEOM along the imaginary time path to obtain the thermal equilibrium state of a system strongly coupled to a non-Markovian bath. Then, we show that the steady state hierarchy elements calculated from the real-time HEOM can be expressed in terms of the hierarchy elements calculated from the imaginary-time HEOM. Moreover, we find that the imaginary-time HEOM allow us to evaluate a number of thermodynamic variables, including the free energy, entropy, internal energy, heat capacity, and susceptibility. The expectation values of the system energy and system-bath interaction energy in the thermal equilibrium state are also evaluated.
We consider a quantum mechanical system represented in phase space (referred to hereafter as "Wigner space"), coupled to a harmonic oscillator bath. We derive quantum hierarchal Fokker-Planck (QHFP) equations not only in real time but also in imaginary time, which represents an inverse temperature. This is an extension of a previous work, in which we studied a spin-boson system, to a Brownian system. It is shown that the QHFP in real time obtained from a correlated thermal equilibrium state of the total system possesses the same form as those obtained from a factorized initial state. A modified terminator for the hierarchal equations of motion is introduced to treat the non-Markovian case more efficiently. Using the imaginary-time QHFP, numerous thermodynamic quantities, including the free energy, entropy, internal energy, heat capacity, and susceptibility, can be evaluated for any potential. These equations allow us to treat non-Markovian, non-perturbative system-bath interactions at finite temperature. Through numerical integration of the real-time QHFP for a harmonic system, we obtain the equilibrium distributions, the auto-correlation function, and the first- and second-order response functions. These results are compared with analytically exact results for the same quantities. This provides a critical test of the formalism for a non-factorized thermal state and elucidates the roles of fluctuation, dissipation, non-Markovian effects, and system-bath coherence. Employing numerical solutions of the imaginary-time QHFP, we demonstrate the capability of this method to obtain thermodynamic quantities for any potential surface. It is shown that both types of QHFP equations can produce numerical results of any desired accuracy. The FORTRAN source codes that we developed, which allow for the treatment of Wigner space dynamics with any potential form (TanimuranFP15 and ImTanimuranFP15), are provided as the supplementary material.
We investigate the application of the imaginary time hierarchical equations of motion method to calculate real time quantum correlation functions. By starting from the path integral expression for the correlated system-bath equilibrium state, we first derive a new set of equations that decouple the imaginary time propagation and the calculation of auxiliary density operators. The new equations, thus, greatly simplify the calculation of the equilibrium correlated initial state that is subsequently used in the real time propagation to obtain the quantum correlation functions. It is also shown that a periodic decomposition of the bath imaginary time correlation function is no longer necessary in the new equations such that different decomposition schemes can be explored. The applicability of the new method is demonstrated in several numerical examples, including the spin-Boson model, the Holstein model, and the double-well model for proton transfer reaction.
We have developed a computer code for the thermodynamic hierarchical equations of motion derived from a spin subsystem coupled to multiple Drude baths at different temperatures, which are connected to or disconnected from the subsystem as a function of time. The code can simulate the reduced dynamics of the subsystem under isothermal, isentropic, thermostatic, and entropic conditions. The extensive and intensive thermodynamic variables are calculated as physical observables, and Gibbs and Helmholtz energies are evaluated as intensive and extensive work. The energy contribution of the system-bath interaction is evaluated separately from the subsystem using the hierarchical elements of the hierarchical equations of motion. The accuracy of the calculated results for the equilibrium distribution and the two-body correlation functions is assessed by contrasting the results with those obtained from the time-convolution-less Redfield equation. It is shown that the Lindblad master equation is inappropriate for the thermodynamic description of a spin-boson system. Non-Markovian effects in thermostatic processes are investigated by sequentially turning on and off the baths at different temperatures with different switching times and system-bath coupling. In addition, the Carnot cycle is simulated under quasi-static conditions. To analyze the work performed for the subsystem in the cycle, thermodynamic work diagrams are plotted as functions of intensive and extensive variables. The C++ source codes are provided as supplementary material.
We have developed a computer code for the thermodynamic hierarchical equations of motion derived from a spin subsystem coupled to multiple Drude baths at different temperatures, which are connected to or disconnected from the subsystem as a function of time. The code can simulate the reduced dynamics of the subsystem under isothermal, isentropic, thermostatic, and entropic conditions. The extensive and intensive thermodynamic variables are calculated as physical observables, and Gibbs and Helmholtz energies are evaluated as intensive and extensive work. The energy contribution of the system--bath interaction is evaluated separately from the subsystem using the hierarchical elements of the HEOM. The accuracy of the calculated results for the equilibrium distribution and the two-body correlation functions are assessed by contrasting the results with those obtained from the time-convolution-less Redfield equation. It is shown that the Lindblad master equation is inappropriate for thermodynamic description of a spin--boson system. Non-Markovian effects in thermostatic processes are investigated by sequentially turning on and off the baths at different temperatures with different switching times and system--bath coupling. In addition, the Carnot cycle is simulated under quasi-static conditions. To analyze the work done for the subsystem in the cycle, thermodynamic work diagrams are plotted as functions of intensive and extensive variables. The C++ source codes are provided as supplementary material.
We develop a quantum-classical hierarchical equations of motion (QC-HEOM) approach for simulating non-Markovian open quantum systems. The method combines the ensemble-averaged classical path reference of the quantum-classical path integral formalism with a hierarchy of auxiliary quantum influence functionals. By incorporating thermal fluctuations through an ensemble average over reference trajectories, the hierarchy is required to represent only the residual quantum memory associated with the imaginary part of the bath response function. Consequently, unlike conventional hierarchical equations of motion, QC-HEOM does not require Matsubara or Padé expansions of the thermal kernel and exhibits only weak temperature dependence of the hierarchy size. Furthermore, because thermal fluctuations are supplied through reference classical trajectories, the framework naturally extends beyond harmonic baths and enables the incorporation of anharmonic and molecular environments through externally generated trajectories. We derive the formalism and demonstrate its exactness for a harmonic bath. Applications to an asymmetric spin-boson model and the seven-site Fenna--Matthews--Olson complex illustrate the accuracy of QC-HEOM. It reproduces benchmark quasi-adiabatic path integral and hierarchical equations of motion results while requiring substantially fewer auxiliary objects, particularly at low temperatures. These results establish QC-HEOM as an efficient framework for treating residual quantum memory in quantum-classical descriptions of open-system dynamics. The separation of thermal fluctuations from residual quantum memory through the use of Wigner trajectories provides an approximate route toward hierarchical treatments of complex anharmonic environments that are inaccessible to conventional HEOM approaches.