A statistical characteristics preserving watermarking scheme for time series databases

1.   Introduction

Over the past decades, the rapid developments of science and technology have allowed people to explore the quantum world. Research on quantum impurity systems (QIS) has become a significant subject and has received extensive attention for its widespread applications in physics^{[1, 2]}, chemistry^[3–6], nanoelectronics^{[7, 8]}, quantum information processing^{[9, 10]}, spintronics^{[11, 12]}, and quantum computation^{[13, 14]}. QIS are typical open systems involving few degrees of freedom, which differs from the surrounding material environment. The material environment contains a huge (practically infinite) number of near-degenerate degrees of freedom, which usually serves as an electron reservoir and thermal bath. The strongly correlated states of QIS mentioned in this review originate from the strong coupling between the system and environment. A typical and famous example is the Kondo state that emerges in the low-temperature regime^{[15, 16]}. Physically, the appearance of Kondo states originates from the local spin moment of the system being screened by the spins of itinerant electrons in the surrounding environment^[17]. Because of the presence of the Kondo state, the strongly correlated QIS exhibit diversified novel quantum phenomena. These include the Kondo effect^{[18, 19]}, local magnetic anisotropy^[16], quantum memristive effect^[20–22], unconventional superconductivity^[23], and thermoelectric effect^{[24, 25]}.

To date, there have been enormous experimental efforts on the precise measurement and control of QIS, such as molecular magnets, nanoclusters, quantum dots, etc. For instance, single molecular magnets can be used to construct spin qubits, and hence, they are expected to become the building blocks of quantum storage devices^[26]. Another example is that the manipulation of local quantum states in magnetic molecules has been realized by using scanning tunneling microscopy (STM)^{[16, 27, 28]} and mechanically controlled break junction (MCBJ) setups^{[29, 30]}. These intriguing experimental findings have aroused increasing interest and enthusiasm from researchers. To better investigate and understand the physical origin of experimental observations, the development of theoretical methods has become essential. However, accurate simulations of strongly correlated QIS still face multifold challenges. First, as an open quantum system, the exchange of energy and particles with the environment is inevitable, which leads to a huge computational complexity. Second, conventional perturbation methods have difficulty treating strong system-bath coupling problems and multielectron cotunnelling processes. Third, QIS itself may involve complicated electron-electron and vibronic coupling interactions. Moreover, the evolution of a system often reflects the non-Markovian memory of the environment. Therefore, an accurate, efficient and universal theoretical method capable of tackling all these challenges and exactly characterizing strongly correlated QIS is highly desirable.

A vast amount of theoretical effort has been devoted to developing formally explicit nonperturbative methods. These include the numerical renormalization group (NRG) method^{[31, 32]}, the density matrix renormalization group (DMRG) method^[33] and its time-dependent extension (TD-DMRG)^[34–36], the exact diagonalization (ED) method^[37], the quantum Monte Carlo (QMC) method^[38–40], the real-time path integral (PI) method^[41–43], the iterative summation of path integrals (ISPI) method^[44–46], Green’s function method^[47–50], the stochastic equation of motion (SEOM) method^[51–53], the multiconfiguration time-dependent Hartree (MCTDH) method^{[54, 55]} and its multilayer extension (ML-MCTDH) method^[55] and the second-quantized version (ML-MCTDH-SQR) method^{[56, 57]}. Although these methods have been successfully applied to some specific physical models, they are still restrained by certain intrinsic limitations. For example, the NRG method is only applicable to the normalizable Hamiltonian form and faces difficulties in characterizing systems far from thermal equilibrium^{[40, 49]}. The QMC method often suffers from the dynamical sign problem, and its performance in the low temperature regime is unsatisfactory^{[47, 49]}. The time-dependent versions of the DMRG, QMC, and PI methods have difficulties in obtaining converged results in the long-time regime, and the MCTDH method often encounters numerical problems when it is used to handle a complicated environment^[40]. In addition, because of the limited computer resources and memory storage space at present, the above methods are only suitable for a number of physical models^[50]. Therefore, it is vital to find an efficient, accurate, and universal approach to investigate QIS and provide important theoretical insights into experimental findings.

An alternative numerically exact approach called the hierarchical equations of motion (HEOM) method has emerged as a popular method, which is based on quantum dissipation theory (QDT) and formulated in terms of the reduced density matrix. It is a powerful numerical approach to describe the reduced system dynamics and steady-state properties of QIS coupled to a non-Markovian and nonperturbative environment. The mathematical form of the HEOM is a hierarchical set of linear differential equations, whose basic variables are the reduced system density operator and auxiliary density operators (ADOs); see Fig. 1. These variables contain the key information about the evolution of the quantum states of the systems. In principle, the HEOM formalism is formally rigorous, as long as the bath environment satisfies Gaussian-Wick statistics, such as free electron reservoirs. It is capable of revealing the underlying physics of the strong electron-electron interactions, the fluctuation-dissipation relation, and non-Markovian memory in a nonperturbative manner, and it can be directly applied to the numerical simulations of various QIS, thereby accurately evaluating system observables, including the static and dynamic response properties in both equilibrium and nonequilibrium situations.

Figure  1.  Schematic diagram for the hierarchal structure of the HEOM. Each green discs represents a density operator. $\rho^{(0)}$ is the reduced-system density operator and $\rho^{(n>0)}$ is the $n$th-tier ADOs. Purple double-arrows denote the couplings between two density operators. The size of the hierarchy space is set by the horizontal M and vertical L dimensional parameters. M is the number of basis functions to unravel the environmental memory and L denotes the truncation tier. Reprinted with permission from Ref. [114]. Copyright 2018, American Institute of Physics.

DownLoad: Full-Size Img PowerPoint

In 1989, Tanimura and Kudo first proposed the HEOM formalism to study an open quantum system coupled to a high temperature bath^[58]. Subsequently, they extended this method to lower temperatures by adding correction terms^[59–61]. The HEOM formalism has also been derived by Shao et al. based on stochastic decoupling of system-bath dissipative interactions^[62] and by Xu and Yan et al. based on the Feynman-Vernon path-integral formulation^{[63, 64]}. Later, Tanimura and coworkers developed real-time and imaginary-time HEOM to calculate thermodynamic variables^{[65, 66]}. Shi and Yan proposed an efficient filtering algorithm to reduce the number of ADOs, which significantly improves the computational efficiency of the HEOM method^[67]. Wu et al. developed an extended HEOM to study the dynamics of the spin-boson model in near-zero temperature regimes^[68–70]. Shi and coworkers formulated the HEOM method in the framework of matrix product states and improved its numerical stability^{[71, 72]}. The HEOM method has been widely applied to study open quantum systems coupled to bosonic environments, and considerable progress has been made in the study of nonlinear spectroscopy^[73–75], quantum information^[75–77], and quantum heat transport^[78]. In 2008, Jin, Zheng, and Yan extended the HEOM method to open quantum systems linearly coupled to fermionic reservoirs^[79], and later to strongly correlated QIS, which are the main subject of this paper^[80]. This has enabled comprehensive studies on the Kondo effect at low temperatures^{[19, 22, 81–83]}. In particular, Zheng and coworkers have developed a computer program—HEOM for QUantum Impurity with a Correlated Kernel (HEOM-QUICK)—which provides an accurate, efficient and universal tool for the characterization of strongly correlated QIS^[84]. Moreover, the HEOM-QUICK program has been utilized in conjunction with the density functional theory (DFT) method to conduct first-principles-based simulations on the local quantum states of single molecular junctions, which have offered useful insights into a variety of novel quantum phenomena^{[19, 81, 85, 86]}.

In addition to the above efforts, impressive progress has been achieved by other authors. For instance, Härtle et al. employed the HEOM method to simulate nonequilibrium electron transport with vibronic couplings through nanosystems^{[87, 88]}, and investigated the nonequilibrium steady-state properties^{[57, 89, 90]}. In recent years, B$\ddot {\rm{a}}$tge and Thoss have derived the HEOM formalism for open quantum systems coupled with both bosonic and fermionic environments and employed it to reveal the intricate interplay among electronic and vibrational degrees of freedom in quantum transport processes^[91]. Additionally, Klaus et al. implemented the HEOM method with shared memory parallel programming techniques to improve the computational efficiency and enable simulation for large complex systems^[92]. Shi and coworkers have applied the Tucker and hierarchical Tucker tensors to efficiently construct the HEOM^[93].

It is worth noting that different theoretical methods have inherent advantages and disadvantages. For instance, while the NRG method has been very successful in studying the low-energy equilibrium properties of QIS^[32], it often encounters problems in the exploration of high-energy nonequilibrium properties^[40]. The DMRG, TD-DMRG and MCTDH methods are capable of characterizing the steady-state properties of QIS and addressing quantum transport problems^[49]. However, the study of long-time dynamics behavior using these methods requires an intensive computational cost^[75]. The ISPI method allows for the investigation of long-time propagation and complex environments^{[40, 45, 46]}, but it often involves the cutoff of environmental memory^[94–96]. The HEOM method is capable of resolving the combined effects of time-dependent external sources, system-environment dissipation, non-Markovian memory, and many-body correlation in a nonperturbative manner^{[65, 66, 75]} and hence provides a highly accurate, efficient and universal tool for characterizing the stationary and dynamic properties of strongly correlated QIS^{[75, 84]}.

Despite the rapid progress and emerging applications of the HEOM method, it still faces several challenges. First, for extremely low temperature regimes where long-time memory should be considered, it is necessary to introduce high-order ADOs and use more memory basis functions to resolve multielectron cotunnelling process and non-Markovian memory effect. This inevitably leads to enormous computational costs^{[40, 49, 97–99]}. Second, the hierarchical space needs to be truncated at a finite size in practical calculations, which not only brings numerical errors, but also makes the Liouville matrix in the HEOM space slightly ill-conditioned so that the numerical instability problem emerges in the calculation of stationary-state properties and long-time dissipative dynamics^{[22, 100, 101]}. Over the past few years, researchers have made numerous efforts to develop novel schemes to overcome these challenges. These schemes generally fall into two categories: schemes for decomposing the environmental memory and schemes for truncating the hierarchy. The former is devoted to reproducing the exact reservoir correlation functions to minimize the horizontal dimension $M$ of the hierarchy in HEOM theory^{[22, 98, 102–105]}, while the latter is targeted at reducing the vertical dimension $L$ as much as possible without compromising the accuracy of the output^{[61, 100, 106]}; see Fig. 1. These developments in methodology significantly improve the efficiency and applicability of the HEOM method, especially in cryogenic temperature regimes. Moreover, their implementation greatly reduces the computational cost and alleviates the numerical instability problem in the time evolution calculation.

This paper aims to review the recent theoretical advances in the development of fermionic HEOM method and some of its representative applications. Emphasis will be placed on realizing highly efficient and accurate HEOM calculations at cryogenic temperatures as well as the characterization of strongly correlated QIS. The remainder of this paper is organized as follows. In Section 2, we review the construction of the HEOM formalism, along with recent advances in spectrum decomposition schemes and terminators. The usefulness and practicality of the fermionic HEOM method are exemplified by studies on nonequilibrium quantum transport, precise control of local quantum states, and quantum thermodynamics problems in Section 3. Concluding remarks and perspectives on the HEOM method are given in Section 4.

2.   Advances in the development of HEOM method

2.1   Modeling quantum impurities and electron reservoir environment

Let us consider a typical open quantum system embedded in a noninteracting fermionic reservoir environment consisting of free electrons. The total system-and-environment Hamiltonian reads (for brevity, we set $\hbar = k_{\rm{B}} = e = 1$ hereafter)

Here, $\hat H_{\rm{sys}}$ is the system Hamiltonian, and $\hat H_{\rm{res}} = \displaystyle \sum \limits_{\alpha}\hat h_{\alpha} = \displaystyle \sum \limits_{\alpha k s} \epsilon_{\alpha k s}\hat{d}^{\dagger}_{\alpha k s}\hat{d}_{\alpha k s}$ is the noninteracting reservoir Hamiltonian, with $\hat{d}^{\dagger}_{\alpha k s}(\hat{d}_{\alpha k s})$ being a spin-$s$ electron creation (annihilation) operator for the $k$th state with energy $\epsilon_{\alpha k s}$ in the $\alpha$-reservoir. $H_{\rm{coupling}} = \displaystyle \sum \limits_{\alpha v s}(\hat{a}^\dagger_{vs}\hat{F}_{\alpha v s}+{\rm{H.c.}}$) represents the system-reservoir coupling Hamiltonian, with $\hat{F}_{\alpha v s} = \displaystyle \sum \limits_{k}t_{\alpha k v s}\hat{d}_{\alpha k s}$, where $t_{\alpha k v s}$ is the coupling strength between the $v$th orbital of the impurity system and the $k$th state of the $\alpha$-reservoir. The impurity system usually involves a small number of degrees of freedom. A key quantity to characterize open quantum systems is the system-reduced density matrix $\rho$, which is calculated by taking the partial trace of the total density matrix $\rho_{\rm{T}}$ over the reservoir subspace:

Reservoirs are usually treated by quantum statistics methods; that is, the grand canonical density operator is utilized to describe the thermodynamic equilibrium state of the $\alpha$-reservoir^[79], which has the form of

Here, $\hat N_{\alpha} = \displaystyle \sum \limits_{k}\hat{d}^{\dagger}_{\alpha k s}\hat{d}_{\alpha k s}$ is the occupancy number operator of the $\alpha$-reservoir and $\mu_{\alpha}$ is the chemical potential of the $\alpha$-reservoir. $\beta_{\alpha} = 1/T_{\alpha}$, where $T_{\alpha}$ represents the temperature of the $\alpha$-reservoir.The reservoir spectral density functions assume a Lorentzian function form of

where $\Gamma_{\alpha}$ is the coupling strength between impurities and reservoir $\alpha$, and $\Omega_{\alpha}$ and $W_{\alpha}$ represent the band center and the band width of the $\alpha$-reservoir, respectively. It is evident that the linearly coupled noninteracting electron reservoirs satisfy the Gaussian statistics, and thus, the influence of the $\alpha$-reservoir on the dynamic properties of systems is completely captured by the reservoir correlation functions

Here, the notation $\sigma = +$ or − is introduced and $\bar \sigma \equiv -\sigma$ represents the opposite sign of $\sigma$. $\hat{F}_{\alpha v s}^{+}(t) \equiv {\rm{e}}^{i{\hat H}_{\rm{res}}t}\hat{F}_{\alpha v s}^{\dagger}{\rm{e}}^{-i{\hat H}_{\rm{res}}t}$ and $\hat{F}_{\alpha v s}^{-}(t) \equiv {\rm{e}}^{i{\hat H}_{\rm{res}}t}\hat{F}_{\alpha v s}{\rm{e}}^{-i{\hat H}_{\rm{res}}t}$ are the Heisenberg operators. The reservoir correlation functions obey the time-reversal symmetry and detailed-balance relations, i.e.,

More importantly, the expression of the reservoir correlation functions can further be associated with the spectral density functions $J_{\alpha s v v^{\prime}}(\omega)$ and the Fermi distribution function $f_{\alpha}^{\sigma}(\omega)$ via the fluctuation-dissipation theorem in the fermionic grand canonical ensembles,

where $f_{\alpha}^{\sigma}(\omega)\equiv 1/[1+{\rm{e}}^{\sigma \beta_{\alpha}(\omega-\mu_{\alpha})}]$ is the Fermi distribution function with the inverse temperature $\beta_{\alpha}$ and the chemical potential $\mu_{\alpha}$.

2.2   An overview of the fermionic HEOM formalism

In the following, the formally rigorous HEOM formalism is constructed through the Feynman-Vernon influence functional path-integral theory with Grassmann algebra. The construction starts with an initial factorization ansatz with the system decoupled from the reservoirs, which is expressed as $\rho_{\rm{T}}(t_0) = \rho(t_0)\rho_{\rm{res}}$. This ansatz is considered to be exact as long as the initial time $t_0$ takes the infinite past, i.e., $t_0\rightarrow -\infty$. Denote ${\cal{U}}(t,t_0)$ as the reduced propagator in the Liouville space. At the present time $t$, the system-reduced density operator $\rho(t)$ can be described via $\rho(t)\equiv{\cal{U}}(t,t_0)\rho(t_0)$. The path-integral expression of the reduced Liouville-space propagator has the form of^{[79, 107]}

where $S[\psi]$ is the reduced-system classical action function, with $\{\psi\}$ being an arbitrary basis set defined in the system space, and ${\boldsymbol{\psi}} \equiv(\psi,\psi^{\prime})$. ${\cal{F}}[{\boldsymbol{\psi}}]$ is the Feynman-Vernon influence functional, which captures the effect of the environment on the reduced system. By using Wick’s theorem on the thermodynamic Gaussian average, the influence functional can be exactly expressed as

Here, ${\cal{R}}$ is the dissipation functional. The Grassmann variable ${\cal{B}}_{\alpha v s}^{\sigma}$ is a generating functional related to a specific principal dissipative mode, which is expressed as

with ${B}_{\alpha v s}^{\sigma}[t;\{\psi\}] \equiv \displaystyle \sum \limits_{v^{\prime}} \int_{t_0}^{t}{\rm{d}}\tau C_{\alpha s v v^{\prime}}^{\sigma}(t,\tau)a^{\sigma}_{v^{\prime} s}[\psi(\tau)]$ and ${B}_{\alpha v s}^{\prime \sigma}[t;\{\psi^{\prime}\}] \equiv \displaystyle \sum \limits_{v^{\prime}} \int_{t_0}^{t}\nonumber {\rm{d}}\tau C_{\alpha s v v^{\prime}}^{\bar \sigma^{\ast}}(t,\tau)a^{\sigma}_{v^{\prime} s}[\psi^{\prime}(\tau)]$. Hereafter, a multicomponent index $j = \{\sigma \alpha s v v^{\prime} m\}$ is introduced to label a principal dissipation mode for brevity. To construct a formally closed HEOM, a sum-overpoles (SOP) representation is adopted to unravel the reservoir correlation functions, i.e.,

Here, the reservoir correlation functions are expanded by a series of exponential functions, where the complex parameters $\{\eta_{j}\}$ and $\{\gamma_{j}\}$ represent the expansion coefficients and exponents. $M$ is the total number of poles, which is given by summing the number of poles of the Fermi distribution function $P$ and the reservoir spectral density functions $Q$, i.e., $M = P+Q$. Each individual exponential term in Eq. (12) makes the generating functional ${\cal{B}}_{j}$ in Eq. (11) become

with ${ B}_{j}[t;\{ \psi\}] \equiv \sum_{v^{\prime}} \int_{t_0}^t {\rm{d}}\tau {\rm{e}}^{-\gamma_{j}(t-\tau)}a^{\sigma}_{v^{\prime} s}[\psi(\tau)]$. The time derivative of ${\cal{B}}_{j}$ satisfies

Similarly, a series of auxiliary density operators (ADOs) $\{\rho_{j_1 \cdots j_n}^{(n)}; n = 1,\cdots,L\}$ with the truncation tier $L$ can be defined via

where the path-integral expression of these auxiliary propagators ${{\cal{U}}}^{(n)}_{j_1 \cdots j_n}$ is

with ${\cal{F}}^{(n)}_{j_1 \cdots j_n}\equiv{\cal{B}}_{j_n}\cdots{\cal{B}}_{j_1}{\cal{F}}$.

Taking the time derivatives of Eq. (15), the final closed HEOM for $n$th-tier ADOs have a compact form of^{[79, 80, 84, 100]}

Eq. (17) clearly exhibits the relations between the $n$th-tier ADOs and its tier-up and tier-down ADOs. The structure of the hierarchy is shown in Fig. 1. System reduced density operator $\rho\equiv \rho^{(0)}$ and the ADOs $\{\rho_{j_1 \cdots j_n}^{(n)}; n = 1,\cdots,L\}$ in the hierarchy obey the Grassmann parity, and they are the key quantities to characterize the impurity system. In Eq. (17), ${\cal{L}}_{\rm{sys}}(t)$ represents the Liouville-space reduced system operator, defined as ${\cal{L}}_{\rm{sys}}(t)\hat O \equiv [\hat H_{\rm{sys}}(t), \hat O]$ for an arbitrary fermionic operator $\hat O$. The damping parameters $\gamma^{(n)}_{j_1 \cdots j_n}(t) = \displaystyle \sum \limits_{r = 1}^{n}\gamma_{j_r}(t)$, where $\gamma_{j_r}(t)$ is a characteristic dissipation rate corresponding to the relevant exponent in Eq. (12). In addition, ${\cal{C}}_{j}\equiv {\cal{C}}_{\alpha s v v^{\prime} m}^{\sigma}$ and ${\cal{A}}_{\bar j}\equiv {\cal{A}}_{\mu s}^{\bar \sigma}$ are superoperators that satisfy ${\cal{C}}_j \hat O = \eta_j a_{v s}^{\sigma}\hat O+\eta_j^{\ast}\hat O a_{v s}^{\sigma}$ and ${\cal{A}}_j \hat O = [a_{\mu s}^{\sigma},\hat O]$, respectively.

2.3   Decomposition of environmental memory

As shown in Fig. 1, the structure of the hierarchy is exhibited as a two-dimensional pyramid, which consists of the reduced-system operator and the ADOs. The horizontal dimension of the pyramid resolves all the available principal dissipation modes between the system and environment, and its width is determined by the number of basis functions $M$ in the SOP representation. The vertical dimension resolves all the possible dissipation events occurring through the available dissipation modes, and its height is determined by the value of truncation tier $L$. Generally, $M$ and $L$ together set the size of the hierarchy space, and thus determine the computational and memory cost of the HEOM calculations. It is highly desirable to minimize the computational cost while obtaining accurate results; thus, a vast amount of effective spectrum decomposition schemes and terminators are successively proposed to reduce the horizontal dimension $M$ and the vertical dimension $L$, respectively. It is well known that the decomposition of environmental memory is based on the expansion of the reservoir correlation functions $C_{\rm{res}}(t)$ (see Eq. (12)), i.e., the SOP unraveling of the Fermi function $f_{\alpha}^{\sigma}(\omega)$ and the reservoir spectral density function $J_{\alpha s v v^{\prime}}(\omega)$. In earlier years, various spetrum decomposition schemes have been proposed and developed, including the Matsubara spectrum decomposition (MSD) scheme^[79], the partial fractional decomposition (PFD) scheme^[108], a hybrid spectrum decomposition and frequency dispersion scheme^[109], and the Padé spectrum decomposition (PSD) scheme^{[102, 103, 110]}. Despite the success of these schemes, a common disadvantage of these schemes is that they have a rather short accuracy length at extremely low temperatures, thus hindering the investigation of strong correlation problems such as the Kondo effect. Therefore, an effective and accurate spectrum decomposition scheme to eliminate the low-temperature curse is highly desirable.

To this end, several novel schemes for the decomposition of environmental memory have been proposed to accurately reproduce the reservoir correlation functions. These include the low-frequency logarithmic discretization (LFLD) scheme^[104], the Fano spectrum decomposition (FSD) scheme^{[22, 98]}, the time-domain Prony fitting decomposition (t-PFD) scheme^[105], etc.

LFLD scheme. It was proposed in 2017 for the efficient decomposition of the fermionic reservoir spectrum^[104]. This scheme subtly combines the advantages of the minimum-dissipaton ansatz^[111] and the PSD scheme and requires a smaller number of exponential functions to expand the reservoir correlation functions to greatly improve the efficiency and accuracy of the HEOM method at low temperatures.

FSD scheme. Subsequently, the FSD scheme has been further put forward, which aims at overcoming the discontinuity of Fermi/Bose functions at $\omega = 0$ near zero temperature and thus allowing the HEOM method for cryogenic measurement^{[22, 98]}. This scheme is applicable for both fermionic and bosonic reservoir environments, although we only review the fermionic case later. The emphasis of the FSD scheme is placed on the efficient and accurate SOP expansion of Fermi distribution functions, that is, fragmenting the Fermi distribution function into two parts: a high-temperature reference and a low-temperature correction. The high-temperature reference part $f_{\rm{F}}^{\rm{PSD}}(\omega,T_0)$ is decomposed accurately via the standard PSD scheme as follows

where $T_0$ is an artificial reference temperature that is used to calibrate the boundary between the high-temperature reference and low-temperature correction. $P_0$ is the number of ${\rm{Pad\acute{e}}}$ first-order poles, which usually depends on the selection of the reference temperature $T_0$. The positive real parameters $\{\eta_p\}$ and $\{\varepsilon_p\}$ can be numerically given with high precision. The remaining part is called the low-temperature correction $\Delta f_{\rm{F}}(\omega,T,T_0)$ is mainly responsible for the discontinuity in the vicinity of $\omega = 0$ and influences the convergence speed of the results. It can be approximatively expressed as

where $\{ F_d^{j_d}(\omega,T_0;a_d,b_d); d = 1,\cdots,D\}$ are a set of modified Fano functions with the following expression

Here, the positive integers $\{j_d\}$ specify the orders of the Fano poles; the fitted real numbers $\{a_d\}$ and $\{b_d\}$ are dimensionless and temperature dependent. The resultant expansion of the Fermi function $f^{\rm{FSD}}_{\rm{F}}(\omega, T)$ is the sum of Eq. (18) and Eq. (19), has the form of

with the effective number of poles being $P_{D} = P_0 +\displaystyle \sum \limits_{d = 1}^{D}j_d$. It is worth noting that the high-order Fano poles with $j_d>1$ would cause extra polynomial-exponential terms to expand reservoir correlation functions, which improves the expression shown in Eq. (12). Consequently, with the FSD scheme, the reservoir correlation functions are unraveled by a series of polynomial exponential functions, which can be expressed as

Here, $\{m_p\}$ are nonnegative integers, and $P = P_J+P_{D}$ is the total number of memory basis functions, where $P_J$ represents the number of poles for unraveling the reservoir spectral density function $J(\omega)$. It has been confirmed by Fig. 2 that the FSD scheme is superior to the standard PSD scheme for low-temperature expansion of the Fermi distribution function. Fig. 2a depicts the expansion of the Fermi distribution function based on the PSD scheme and the FSD scheme, as well as its exact results at a rather low temperature. Apparently, relatively accurate results can be obtained by the FSD scheme with $P_0 = 20$ for the high-temperature reference part and $\displaystyle \sum \limits_{d = 1}^{D}j_d = 9$ for the low-temperature correction part, while this requires $P>600$ for the standard PSD scheme. Moreover, with the same number of total poles, the FSD scheme has a much longer accuracy length than the standard PSD scheme. Fig. 2b exhibits the comparison between the low-temperature correction part in the FSD scheme and the exact results. The negligible numerical error verifies the usefulness of the Fano function to expand the low temperature correction part. In general, these results show the superiority of the FSD scheme over other spectrum decomposition schemes for expanding the Fermi distribution function, especially in extremely low temperature regimes.

Figure  2.  (a) The expansion of the Fermi function based on the PSD scheme (green/blue dotted line), the FSD sheme (yellow/red solid line), and the exact results (empty triangles). (b) The low-temperature correction part in the FSD scheme (red solid line) expanded by the Fano functions and the numerical error (blue solid line). The original and reference temperatures are in unit of $\Gamma$: $T = 0.001$, $T_0 = 0.1$. Reprinted with permission from Ref. [98]. Copyright 2019, American Institute of Physics.

DownLoad: Full-Size Img PowerPoint

On the one hand, the FSD scheme requires fewer basis functions to unravel the non-Markovian memory of the environment, which accelerates computational speed and reduces the memory cost. On the other hand, the use of the FSD scheme indeed yields highly accurate results, even in extremely low temperature regimes. Hence, the FSD scheme is an advanced and reliable spectrum decomposition scheme to expand the reservoir correlation functions.

t-PFD scheme. Recently, a highly accurate, efficient and universal exponential series method called the time-domain Prony fitting decomposition (t-PFD) scheme has also been proposed^[105]. This scheme is capable of directly unraveling arbitrary reservoir correlation functions $C_{\rm{res}}(t)$ to the following form:

Here, $\{\zeta_m;\lambda_m\}$- and $\{\zeta_{m^\prime};\lambda_{m^\prime}\}$-parameters are obtained by fitting the real part $C^{(r)}_{\rm{res}}(t)$ and the imaginary part $C^{(i)}_{\rm{res}}(t)$, respectively, with only one of the parts being temperature dependent. The total number of poles is $M = M_{r}+M_{i}$, and all of them are first-order poles. Whatever the form of spectral density functions $J(\omega)$, the fitting results for the frequency-domain reservoir correlation functions in the vicinity of $\omega = 0$ via the t-PFD scheme always have smaller errors than those in conventional schemes. Moreover, the numerical efficiency and accuracy of the HEOM method based on the t-PFD scheme is significantly improved, especially in cryogenic temperature regimes. In Particular, the advances and superiority of the t-PFD scheme have been exemplified by demonstrating the HEOM calculations on the single-impurity Anderson model (SIAM). Therefore, the proposal of the t-PFD scheme is a huge breakthrough for the development of HEOM methodology.

In summary, the rapid progress of spectrum decomposition schemes not only realizes highly accurate and efficient measurements for QIS but also allows the HEOM method to access unprecedentedly low temperature regimes. However, the asymptotic instability problem remains inevitable in HEOM calculations, which will become one of the main issues of future research.

2.4   Terminators

The advances of the spectrum decomposition scheme referred to in the previous section target reproducing exact reservoir correlation functions and optimizing the horizatonal dimension of pyramid $M$. As mentioned before, the hierarchical structure of the HEOM method contains two dimensions. Naturally, the optimization of the vertical dimension is another choice to improve the HEOM method. To date, a variety of terminators have been developed to achieve the goal of exact and effective Lth tier truncation. These include:

Zero-value terminator. This conventional truncation scheme has been achieved in a simple and straightforward way, that is, to set all the ADOs beyond Lth tier $\{\rho_{j_1 \cdots j_n}^{(n>L)}\}$ to zero^{[107, 112, 113]}. Based on the zero-value terminator, Han et al. analyzed and studied the exact truncation tier of the fermionic HEOM method via a termination pattern originating from the Pauli exclusion principle^[114]. For a general system coupled to fermion reservoirs, the exact density operators up to the $n$th tier can be obtained by the $(n+\tilde L)$th-tier truncation with $\tilde L = 2N_\sigma N_v$, where $N_\sigma = 2$ for the fermionic environment and $N_v$ represent the system degrees of freedom. In contrast, for a noninteracting system, the exact truncation tier is only half of that for general systems, i.e, $\tilde L^{\prime} = N_\sigma N_v$. It is commendable that the termination pattern is always effective regardless of the number of reservoirs, the system-environment coupling strength, and the spectrum decomposition scheme. However, an intrinsic shortage of the zero-value terminator is that the influences from the higher-tier($n>L$) ADOs are entirely ignored.

Zero-derivative terminator. Instead of having $\{\rho_{j_1 \cdots j_{L+1}}^{(n>L)} = 0\}$, the zero-derivative terminator has been developed for the hierarchal structure, which sets the first-order time derivative of $(L+1)$th-tier ADOs $\{\dot{\rho}_{j_1 \cdots j_{L+1}}^{(L+1)}\}$ to zero^[106]. Meanwhile, the ${\rho}_{j_1 \cdots j_{L+1}}^{(L+1)}$ naturally satisfy the equations of motion (EOM) of $(L+1)$th-tier ADOs (see Eq. (17)). Combining these two equations, the value of ${\rho}_{j_1 \cdots j_{L+1}}^{(L+1)}$ can be evaluated by

This expression is then substituted into the EOM of the Lth tier ADOs, which gives rise to a truncation for the hierarchal structure. Clearly, the influences from higher-tier ADOs are partially considered in the terminal tier of the hierarchy. Moreover, the use of a zero-derivative terminator does not require more memory cost to store the ADOs since no extra basic variables are invoked during calculations. Usually, the ADOs are sparse matrices, which means that a large percentage of elements in the ADO matrix are exactly zero. By taking advantage of this feature, the memory cost can be greatly reduced. To this end, they developed an efficient algorithm consisting of two ansatz to quickly locate and screen out all zero elements. Consequently, the proposal of a zero-derivative terminator and the use of sparsity of ADO matrices significantly improve the efficiency of the HEOM method. Afterwards, Kasper and Thoss combined the zero-derivative terminator with an iterative approach to yield an efficient steady-state solver for the HEOM method and thus simulated nonequilibrium quantum transport processes in molecular junctions^[90].

Adiabatic terminator. Recently, an adiabatic terminator was proposed to further improve the numerical performance of the fermionic HEOM method^[100]. Zheng and coworkers found that the distribution of the poles of reservoir correlation functions unraveled by the FSD scheme is sufficiently discrete. This implies that dissipation rates of different dissipation modes are distinguishable. Inspired by the well-known Born-Oppenheimer approximation, which separates the motion of electrons and nuclei due to their difference in motion speed. The adiabatic terminator has been realized by decoupling the principal dissipation mode with the fastest dissipation rate from the other slow ones.

The HEOM hierarchy truncated at the first tier corresponds to a quantum master equation obtained by adopting the self-consistent Born approximation as follows:

with $\{\tilde \rho_j\}$ being the auxiliary dynamical variables and can be expressed as

where ${\cal{U}}(t-\tau)$ is the reduced system propagator. Recall the EOM for the zero-tier ADO:

Since Eq. (25) and Eq. (27) describe the same dissipative dynamics, Eq. (26) can be regarded as the adiabatic terminator for the first-tier HEOM. Defining a dissipative Liouville propogator ${\cal{W}}_{jv^{\prime}}\equiv \int_{t_0}^{t}{\rm{d}}\tau \eta_{j}(t-\tau)^{m_p}{\cal{U}}(t-\tau){\rm{e}}^{-\gamma_j(t-\tau)}$ and its conjugate ${\cal{W}}_{\bar jv^{\prime}}^\dagger \equiv \int_{t_0}^{t}{\rm{d}}\tau \eta_{\bar j}^\ast (t-\tau)^{m_p}{\cal{U}}(t-\tau){\rm{e}}^{-\gamma_{\bar j}^\ast (t-\tau)}$, the adiabatic terminator for the first-tier ADOs is then given by

Similarly, the adiabatic terminator for the Lth-tier ADOs can be expressed by a general form as

The hierarchy is finally truncated by substituting the values of ${\rho}_{j_1 \cdots j_n}^{(L+1)}$ into the EOM for the Lth-tier ADOs. Extensive numerical demonstrations have shown that the adiabatic terminator greatly enhances the convergence of calculation results and alleviates long-time numerical instability problems. For example, Fig. 3 depicts the impurity spectral function $A(\omega)$ calculated by the adiabatic terminator and zero-value terminator at different truncation tiers $L$. It is observed that the numerical results calculated by the HEOM method with the zero-value terminator converge at $L = 7$, while that only requires $L = 5$ with the adiabatic terminator. Moreover, in general, the calculation results are considered to be quantitatively accurate as long as convergence is achieved. Therefore, compared with the zero-value terminator, the adiabatic terminator substantially enhances the efficiency and accuracy for the calculations of $A(\omega)$.

Figure  3.  Impurity spectral function $A(\omega)$ of an SIAM at different truncation tiers $L$. The scattered data are obtained by the HEOM calculations with zero-value terminator and adiabatic terminator, respectively. Parameters here are all in unit of $\Gamma$: $T = 0.04$, $\epsilon_{\uparrow} = \epsilon_{\downarrow} = -U/2 = -4$, and $W = 20$. The inset shows the Kondo resonance peak in the vicinity of $\omega = 0$. Reprinted with permission from Ref. [100]. Copyright 2021, American Institute of Physics.

DownLoad: Full-Size Img PowerPoint

In summary, the proposals of these terminators greatly accelerate the convergence of calculation results and alleviate the computational burden in practical applications. However, the simulation of long-time dissipative dynamics in complex systems is still difficult. Hence, a goal for future research is to develop more sophisticated terminators and apply them to investigate more intricate dynamic processes.

3.   Applications of the fermionic HEOM method

Thanks to the continuous development of analytic formalism and numerical algorithms, the HEOM method has been applied to explore a broad range of systems and phenomena, including the local heating effect of quantum dots^{[24, 25, 115]}, nonequilibrium quantum transport through nanosystems^{[87, 88, 91, 116–118]}, strongly correlated Kondo states at molecule/substrate interfaces^{[19, 119, 120]}, dynamic and thermodynamic properties of QIS^{[21, 22, 121–123]}, quantum state evolution^{[70, 71, 106, 124–127]}, nonlinear and multidimensional spectroscopies^{[73, 128–133]}, precise control of spin states of magnetic molecules^{[3, 19, 81, 119, 120]}, and so on.

3.1   Characterization of nonequilibrium quantum transport

In recent years, nonequilibrium quantum transport in molecular junctions or quantum dots (QDs) has received wide attention and has exhibited various intriguing quantum phenomena such as Coulomb blockade^{[116, 134]}, Franck-Condon blockade^{[87, 135, 136]}, strong current fluctuations^{[116, 137, 138]}, and the quantum memristive effect^[20–22]. The HEOM method is one of the most powerful tools to accurately characterize nonequilibrium quantum transport.

For example, the nonequilibrium charge transport in nanojunctions with strong vibronic couplings was numerically exact investigated by Härtle and Thoss in 2016^[87]. The HEOM method has been applied to study a vibrationally coupled transport model with a broad spectrum of parameters, including resonant and off-resonant transport as well as ranging from the adiabatic to nonadiabatic regimes. The calculation results have shown that nonequilibrium effects play an important role in all these regimes. Particularly, for the off-resonant transport regime, the complex interplay of different transport processes and deviations from the commonly known $G_0/2$ rule of thumb are revealed via the analysis of the inelastic cotunneling signal for a vibrational mode in full nonequilibrium. Subsequently, to extend the HEOM method to efficiently treat QIS with a large nonequilibrium vibrational excitation, Thoss and coworkers changed the treatment for vibronic couplings within the HEOM formalism^[88]. Instead of considering the vibrational degrees of freedom of nanosystems as part of the reduced system, they are treated as part of the bath subspace. The HEOM method has been employed to solve a model nanosystem consisting of a single electronic level coupled to two fermionic leads and a single vibrational mode. On the one hand, the convergence behavior with respect to the truncation of electronic and vibrational hierarchy has been investigated in both the adiabatic and nonadiabatic transport regimes. On the other hand, the calculation results have elaborated the phenomenon of vibrational instability, that is, electronic-vibrational coupling decreases with increasing current-induced vibrational excitation.

In addition to the average current, the fluctuations and high-order cumulants of the distribution of transferred charge are also important to characterize and understand electron transport in nanosystems. In the framework of full counting statistics, the numerically exact evaluation of current fluctuations and high-order current cumulants for nonequilibrium charge transport in nanosystems has been realized by the HEOM method^[116]. Based on the calculations for a vibrationally coupled electron transport model, the influence of cotunneling in the off-resonant transport regime has been investigated. Meanwhile, it has been confirmed that the correction to the elastic noise upon opening of the inelastic transport channel is nonuniversal and greatly influenced by the nonequilibrium vibrational excitation. Generally, actual QIS are coupled with both bosonic and fermionic reservoirs and involve many interacting electronic and vibrational degrees of freedom. Based on this fact, Bärtge and Thoss recently reconstructed the HEOM formalism and demonstrated its usefulness for a nanosystem model coupled with multiple fermionic and bosonic reservoirs^[91]. The voltage and thermally driven quantum transport have been investigated by evaluating of the charge current, heat current, and their differential behaviors. Additionally, the calculation results have manifested the intricate interplay of electronic and vibrational degrees of freedom in the nonequilibrium transport scenario.

In addition, numerous efforts have been made to utilize the HEOM method to simulate real-time dynamic evolution processes. For instance, the real-time electric current dynamics of strongly correlated QDs driven by an external periodic voltage has been investigated^[21]. It was found that the dynamic current-voltage (I-V) characteristic pattern is sensitively dependent on the frequency $\Omega$ of the AC voltage. At a low $\Omega$, the dynamic I-V curve behaves close to the steady-state features, while at a high $\Omega$, the shape of the I-V curve is similar to an ellipse, which originates from the phase difference between the current and AC voltage. More importantly, at an intermediate $\Omega$, the dynamic I-V characteristics form a hysteresis loop with self-crossing in the first and third quadrants. This highlights the remarkable memory effects due to Kondo resonance. Moreover, it has been observed that the memory effects are suppressed as the temperature increases, which results in the QDs recovering classical behavior.

Afterwards, with the proposal of the FSD scheme and the adiabatic terminator, the description for the time evolution of the electric current at a relatively lower temperature with a high precision has also been achieved^[22]. Fig. 4a displays the voltage-driven time-dependent electric current calculated by adopting three sets of parameters for the FSD scheme. These arrows mark the time from which the calculation results start to diverge. Clearly, the calculated dynamic current would diverge after a certain time due to the notorious asymptotic instability problem, but the time starting to diverge can be delayed by adopting a more accurate set of parameters. Fig. 4b depicts the time evolution of the voltage-driven electric current calculated by the zero-value terminator and the adiabatic terminator at different truncation tier $L$. The numerical results obtained by HEOM-adia at $L = 5$ diverge at a later time than those obtained by HEOM-zero at $L = 5$. Additionally, the response of the electric current predicted by HEOM-adia at $L = 4$ remains convergent and stable for a long time, as shown in the inset of Fig. 4b. These calculation results again demonstrate the advantages of the novel adiabatic terminator.

Figure  4.  (a) The time evolution of ac voltage-driven electric current. The arrows represent the time of divergence start. The parameters are (in units of $\Delta$): $U = -2\epsilon_d = 12$, $\Delta_L = \Delta_R = 0.5$, $T = 0.05$, $W = 20$, $eV_0 = 1.5$, and $\omega = 0.3$. $\chi$ represents different parameter sets adopted for the FSD scheme, among which $\chi = 1000$ (marked by the star in the legend) is the most accurate. Reprinted with permission from Ref. [22]. Copyright 2020, American Institute of Physics. (b) The variation of real-time electric current $I_L(t)$ driven by ac voltage, flowing from impurity to left reservoir. The results are calculated by HEOM-adia and HEOM-zero at different truncation tier $L$, respectively. The arrows mark the time from which the results start to diverge. The inset shows the long-time evolution of $I_L(t)$ calculated by HEOM-adia truncated at $L$ = 4. The parameters adopted are all in unit of $\Gamma$: $V_0 = 3.0$, $\omega_0 = 1.2$, $T = 0.1$, $\epsilon_{\uparrow} = \epsilon_{\downarrow} = -U/2 = -12$, and $W = 40$. Reprinted with permission from Ref. [100]. Copyright 2021, American Institute of Physics.

DownLoad: Full-Size Img PowerPoint

Therefore, the HEOM method is capable of describing nonequilibrium quantum transport and the real-time evolution dynamics in nanosystems, which provides important theoretical insights into the relevant experiments.

3.2   Characterization of strongly correlated Kondo states and their interplay with local spin excitations

Over the past two decades, strongly correlated Kondo states in molecular junctions have attracted wide attention. The HEOM method has been successively applied to characterize the Kondo states and clarify their competition with local spin excitations^{[19, 82, 139–142]}. In 2014, the investigation of the Kondo phenomena for a magnetic absorption system consisting of a dehydrogenated CoPc molecule adsorbed on an Au(111) surface was realized by combining DFT and the HEOM method^[119]. In this work, the DFT approach is responsible for determining the geometric and electronic structures of the adsorption system, while the HEOM method is utilized to calculate the physical observables by adopting an appropriate SIAM with the relevant parameters extracted from the DFT calculations. A major finding is that the zero-bias peak in the differential conductance $dI/dV$ spectra has obvious temperature dependence in low-temperature regimes. This reveals the presence of strongly correlated Kondo states.

Another magnetic adsorption system, FePc/Au(111), has been investigated by several STM experiments, which confirmed that the measured Kondo features have strong adsorption site dependence^{[29, 139]}. To determine the underlying physics of the experimental observations, the combined DFT+HEOM method has been employed to calculate the $dI/dV$ spectra for the on-top and bridge adsorption configurations^[120]. The simulation results have shown that their Kondo features at zero bias have remarkable differences. For the on-top adsorption configuration, the $dI/dV$ curve resolves a Fano-Kondo antiresonance line shape, which is completely absent for the bridge adsorption configuration. A reasonable explanation for such a difference is that the Fano-Kondo signature in magnetic molecule/metal composites is strongly influenced by the environment-induced coupling.

Later, the combined DFT+HEOM approach was used to investigate the characteristics of local spin excitation in two bilayer composite junctions, i.e., FePc/FePc/Pb(111) and FePc/CoPc/Pb(111)^[140]. Fig. 5 depicts the simulated $dI/dV$ spectra for both composite junctions. In Fig. 5a, two step-like features near $V \simeq \pm 110\,{\rm{mV}}$ are clearly presented, which originate from the inelastic spin-flip excitations between the doubly degenerate ground states $|m_s = \pm1\rangle$ and the excited state $|m_s = 0\rangle$ of the molecule. In contrast, Fig. 5b exhibits four steps located at $V \simeq \pm 30\,{\rm{mV}}$ and $V \simeq \pm 140\,{\rm{mV}}$. The emergence of more steps is ascribed to the presence of nonzero in-plane magnetic anisotropy, which causes splitting of the doubly degenerate ground states $|m_s = \pm1\rangle$ and thus gives rise to two different spin excitation processes. These results provide important theoretical insights into experimental observations and highlight the capability and usefulness of the HEOM method to characterize spin excitations in various QIS.

Figure  5.  The $dI/dV$ spectra calculated by the DFT+HEOM method for the (a) FePc/FePc/Pb(111) and (b) FePc/ CoPc/Pb(111) composites at temperature T = 60 K. Reprinted with permission from Ref. [140]. Copyright 2018, Royal Society of Chemistry.

DownLoad: Full-Size Img PowerPoint

The intricate interplay between strongly correlated Kondo states and local spin excitations has also been studied by the DFT+HEOM method. Fig. 6 clearly demonstrates the differential conductance spectra of a composite junction consisting of a single iron(II) porphyrin molecule adsorbed on the Pb(111) substrate^[81]. The $dI/dV$ spectrum (red curve) depicted in Fig. 6a displays a three-peak lineshape. While the zero-bias peak characterizes the Kondo resonance, the two side peaks at $V = \pm D$ indicate the presence of local spin excitations. The attribution of these peaks can be confirmed by taking the magnetic anisotropy $D = 0$ (gray curve). Fig. 6b displays the $dI/dV$ spectrum (red curve) for a smaller distance between the tip and the magnetic molecule. In this case, the $dI/dV$ spectrum only resolves a single Kondo peak even though $D$ takes the same value as that in Fig. 6a. The absence of the two side peaks is due to the strong Kondo resonance, which overwhelms the spin excitations. Nevertheless, the weak spin excitations still reduce and broaden the Kondo peak, as indicated by a comparison with the gray curve corresponding to the case of $D = 0$. Therefore, the DFT+HEOM method is capable of disclosing the subtle competition between the Kondo resonance and the spin inelastic excitations.

Figure  6.  The differential conductance $dI/dV$ spectra for the SIAM corresponding to the tip/FeOEP/Pb(111) composite with different tip displacement (a) $\Delta z = -350\;{\rm{pm}}$, (b) $\Delta z = -450\;{\rm{pm}}$. The inset of (b) shows the schematic diagram for the energy levels of local spin states. Reprinted with permission from Ref. [81]. Copyright 2018, American Chemical Society.

DownLoad: Full-Size Img PowerPoint

Another example is the physical origin of the unconventional Kondo resonance signatures for the Ni-tip/Co/Cu(100) junction in spin-polarized scanning tunneling spectroscopy has been studied by the DFT+HEOM method^[82]. Fig. 7a and b depict the calculated and experimentally measured $dI/dV$ spectra at different tip displacements $z$, respectively. With $z$ taking a more negative value, the Kondo resonance peak at zero bias gradually splits into two asymmetric peaks, which correspond to two different electron cotunneling processes. The height of asymmetric peaks is sensitively dependent on the hybridization strengths $\Gamma_{\rm{s}}$ and $\Gamma_{{\rm{t}}\sigma}$ ($\sigma = \uparrow \text{or} \downarrow)$, where $\Gamma_{\rm{s}}$ represents the hybridization strength between Co and the substrate Cu(100), while $\Gamma_{{\rm{t}}\uparrow}$ ($\Gamma_{{\rm{t}}\downarrow}$) is the hybridization strength between the Ni tip and Co under a negative (positive) bias. Fig. 7c displays the simulated splitting width $\Delta V$ between two asymmetric peaks against the difference in hybridization strengths $\Delta E = \Gamma_{{\rm{t}}\uparrow}-\Gamma_{{\rm{t}}\downarrow}$. Their linear relationship has been expressed analytically by fitting these data points. The simulation results agree closely with the experimental observations, which highlights the usefulness of the DFT+HEOM method for understanding unconventional Kondo features.

Figure  7.  (a) Calculated and (b) experimentally measured $dI/dV$ spectra at several tip displacements $z$. (c) The splitting width $\Delta V$ between two asymmetric peaks verus the difference of hybridization strengths $\Delta E$. The gray line gives a linearly fitting result. Reprinted with permission from Ref. [82]. Copyright 2022, American Chemical Society.

DownLoad: Full-Size Img PowerPoint

In addition, several experiments have found that the appearance of the Kondo signature usually requires magnetic atoms to adsorb on the fcc or hcp domain of the substrate^[143]. However, in reality, such a requirement is difficult to satisfy because the substrate usually has some inhomogeneous regions. To regularize the Kondo resonance features, a novel strategy has been proposed, and its feasibility has been demonstrated in the Co/Au(111) composite system^[19]. With this strategy, a CoPc molecule is utilized to capture the dispersed Co adatoms absorbed on the Au(111) surface. The CoPc molecule can be regarded as a mold that regularizes the local spin distribution and the Kondo states of the captured Co adatoms. The resultant atom-mold complexes are highly symmetric and have strong $d_\pi$-$\pi$ bonding interactions between the Co adatom and the isoindole units. These are the key reasons to gain uniform and regular Kondo resonance features. More importantly, when the CoPc molecular mold simultaneously captures multiple Co adatoms, fine tuning of the Kondo states is realizable through the long-range superexchange interactions between Co adatoms separated by more than 12 Å. The manipulation processes are simulated by using the DFT+HEOM method, and the numerical results are highly consistent with the experimental observations. This result thus further affirms the practicality and importance of the DFT+HEOM method. Hence, we can conclude that the DFT+HEOM method is a reliable and powerful tool to characterize the Kondo resonance features and reveal their intricate interplay with local spin inelastic excitations.

3.3   Characterization of nonequilibrium thermodynamic properties

In recent years, the quantum thermodynamic properties of nanosystems have aroused great research interest^{[21, 24, 25, 115]}. A novel approach, achieved by measuring the local temperature of a QIS, has been proposed to investigate the local thermodynamic properties^{[144, 145]}. For instance, the zero current condition (ZCC)^{[146, 147]} and minimal perturbation condition (MPC)^{[24, 115, 148]} are two representative conditions that lay the theoretical foundation for determining the local temperatures of nonequilibrium QIS.

The ZCC usually requires the measurement of tiny heat currents, which is rather difficult to realize in experiments, and thus largely limits the applicability of the ZCC^{[144, 147]}. In contrast, MPC is a straightforward protocol that avoids the direct measurement of heat currents^{[24, 115, 144, 148]}. However, MPC is only applicable to single impurity systems. To overcome this limitation, a local MPC (LMPC) has been developed to extend the conventional MPC to open quantum systems consisting of multiple quantum impurities^[25]. In Ref. [25], the HEOM method is employed to calculate various physical quantities, such as the electric current and magnetic susceptibility, and thus determine the local temperature of the QIS. The calculation results suggested that the local temperature $T^{\ast}$ is unique for a single impurity system, and the value of $T^{\ast}$ depends on the specific choice of local observable $\hat O$, particularly in the near-resonance regions.

Fig. 8 depicts the local temperature profile based on the ZCC and the LMPC for a noninteracting four-impurity chain under a thermal bias. For the case of a strong impurity-lead coupling, as shown in Fig. 8a, both the LMPC and the ZCC give the almost linear variation of local temperature $T^{\ast}_i$ with respect to the impurity-$i$, which obeys the classical Fourier’s law. In contrast, Fig. 8b displays the local temperature profile with a weak impurity-lead coupling strength. The ZCC and the LMPC give rise to different predictions for the distribution of local temperature. The local temperature of all four impurities obtained by the ZCC nearly maintains the same value, while that determined by the LMPC exhibits conspicuous oscillations along the chain. The emergence of the oscillations originates from the quantum resonance effects. The inset illustrates that the peaks of the impurity spectral function $A(\omega)$ become higher and sharper as the impurity-lead coupling strength decreases, and the distance between neighboring peaks is clearly influenced by the coupling strength $t$ between two adjacent impurities. This example affirms that the HEOM method is a useful and reliable approach to simulate and understand the measurements of local temperatures, paving the way for the investigation of various thermodynamic properties.

Figure  8.  The distribution of local temperature $T^{\ast}_i$ along a linear chain consisting of four impurities based on the ZCC and the LMPC under a thermal bias with (a) a strong impurity-lead coupling strength $\Delta_\alpha = 0.5\Delta \,(\alpha = L, R)$ and (b) a weak impurity-lead coupling strength $\Delta_\alpha = 0.0125\Delta$. $t$ represents the coupling strength between two adjacent impurities. The inset shows that the impurity spectral function $A(\omega)$ calculated by the HEOM method at different $\Delta_\alpha$ and $t$. Reprinted with permission from Ref. [25]. Copyright 2021, American Physical Society.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions and perspectives

In conclusion, this paper reviews the theoretical construction of a formally rigorous and numerically exact fermionic HEOM method and elaborates its expression based on the Feynamn-Vernon path-integral formulation. With the continuous improvement of methodology over the past decade, the HEOM method has become one of the most popular theoretical approaches to characterize the strongly correlated QIS and reveal the physical origin of experimental phenomena. Specifically, on the one hand, the proposal of a variety of spectrum decomposition schemes including the LFLD, FSD, and t-PFD schemes, has enabled the HEOM method to accurately and efficiently investigate unprecedentedly low temperature regimes. On the other hand, the rapid progress of hierarchical terminators accelerates the convergence speed of calculation results with high precision and thus significantly improves the computational efficiency and saves memory costs. Numerous applications have also been exhibited to confirm the usefulness of the HEOM method and the superiority of the newly proposed schemes and terminators. For example, the HEOM method allows a highly accurate and efficient characterization of nonequilibrium quantum transport in nanosystems. The highly concerned Kondo resonance features and their competition with local spin excitations are clearly disclosed via the DFT+HEOM method. Moreover, the HEOM method plays an important role in the investigation of quantum thermodynamic properties. Therefore, the HEOM method is a universal, powerful and reliable theoretical tool to investigate diversified strong correlation effects in the QIS.

However, the HEOM method still faces several serious challenges, including the rapidly growing computational cost for the simulation of large complex systems and numerical instability problems. These problems severely limit the scope of application of the HEOM method and hinder the investigation of long-term evolutionary behavior. Hence, in future work, on the one hand, better spectrum decomposition schemes and hierarchal terminators are pursued; on the other hand, some novel theoretical frameworks, such as matrix product states proposed by Shi et al.^{[71, 72, 93, 149]}, will be explored to improve the numerical stability and computational efficiency of the HEOM method.

In terms of practical applications, the HEOM method is expected to provide a reliable and useful tool for more emerging fields, such as the precise measurement and control of local spin states in single molecule junctions. Recent technological advances, including spin-polarized STM and the combination of STM with electron spin resonance (ESR) spectroscopy, have enabled the measurement of spin-related properties and spin-spin interactions in single atomic or molecular junctions with unprecedentedly high spatial and energy resolution^[150–153]. Moreover, the combination of STM with pulsed ESR or pump-probe techniques has allowed the real-time control of local spin states in nanostructures^[154–156]. These impressive experimental measurements and observations present new challenges and opportunities for theoretical methods, including the fermionic HEOM method. Work along this direction is underway.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (21973086, 22203083) and the Fundamental Research Funds for the Central Universities (WK2060000018).

Conflict of interest

The authors declare that they have no conflict of interest.