J. Phys. Soc. Jpn. 90, 033001 (2021) [4 Pages]

Open Quantum Dynamics Theory for Non-Equilibrium Work: Hierarchical Equations of Motion Approach

+ Affiliations
Department of Chemistry, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan

A system–bath (SB) model is considered to examine the Jarzynski equality in the fully quantum regime. In our previous paper [J. Chem. Phys. 153, 234107 (2020)], we carried out “exact” numerical experiments using hierarchical equations of motion (HEOM) in which we demonstrated that the SB model describes behavior that is consistent with the first and second laws of thermodynamics and that the dynamics of the total system are time irreversible. The distinctive quantity in the Jarzynski equality is the “work characteristic function (WCF)”, 〈exp[−βW]〉, where W is the work performed on the system and β is the inverse temperature. In the present investigation, we consider the definitions based on the partition function (PF) and on the path, and numerically evaluate the WCF using the HEOM to determine a method for extending the Jarzynski equality to the fully quantum regime. We show that using the PF-based definition of the WCF, we obtain a result that is entirely inconsistent with the Jarzynski equality, while if we use the path-based definition, we obtain a result that approximates the Jarzynski equality, but may not be consistent with it.

©2021 The Author(s)
This article is published by the Physical Society of Japan under the terms of the Creative Commons Attribution 4.0 License. Any further distribution of this work must maintain attribution to the author(s) and the title of the article, journal citation, and DOI.

In thermodynamics, work, W, and heat, \(\Delta Q\), are thermodynamic process quantities, while the internal energy, \(\Delta U\), is an extensive quantity that cannot be measured directly. For a classical system beginning in an equilibrium state, it has been found that the work done under an arbitrary mechanical operation is related to the equilibrium free energy in accordance with the Jarzynski equality:114) \(-{\ln}(\langle\exp[-\beta W(\tau)]\rangle)/\beta =\Delta F_{A}(\tau)\). Here, \(\beta\equiv 1/k_{\text{B}}T\) is the inverse temperature with the Boltzmann constant \(k_{\text{B}}\), \(\Delta F_{A}(\tau)\) is the change in the free energy of the system, \(W(\tau)\) is the non-equilibrium work, and \(\langle\cdots \rangle\) is the ensemble average over all phase space trajectories under the time-dependent external perturbation from time \(t=0\) to τ. Although investigating this equality in the classical regime is straightforward in both theoretical and experimental contexts, doing so in the quantum regime remains challenging,1523) because the dynamics of a small quantum system itself are reversible in time and therefore the system cannot reach thermal equilibrium on its own without a system–bath (SB) interaction: We cannot assume a canonical distribution as the equilibrium state for the system itslef and the heat-bath, due to the presence of the SB interaction. In the present paper, we numerically evaluate the work characteristic function (WCF), \(\langle\exp[-\beta W(\tau)]\rangle\), and \(\Delta F_{A}(\tau)\) with the goal of extending the Jarzynski equality to the fully quantum regime.

A commonly employed model for this kind of investigation is described by a SB Hamiltonian, in which a small quantum system A is coupled to a bath B modeled by an infinite number of harmonic oscillators. We found that the behavior described by this model is consistent with the first and second laws of thermodynamics and provides an ideal platform to examine various fundamental propositions of thermodynamics in the fully quantum regime.24,25) In this model, the Hamiltonian of the total system is given by \begin{equation} \hat{H}(t) = \hat{H}_{A}(t) + \hat{H}_{I} + \hat{H}_{B}, \end{equation} (1) where \(\hat{H}_{A}(t)\), \(\hat{H}_{I}\), and \(\hat{H}_{B}\) are the Hamiltonian of the system, the interaction and the bath, respectively. The system Hamiltonian is given by \(\hat{H}_{A}(t) =\hat{H}_{A}^{0} +\hat{H}_{E}(t)\), with \(\hat{H}_{A}^{0} =\frac{1}{2}\hbar \omega_{0} (|e\rangle \langle e| - |g\rangle \langle g|)\) and \(\hat{H}_{E}(t)=0\) for \(t\leq 0\), where \(|e\rangle\) and \(|g\rangle\) are the excited and ground states of the system, and \(\hat{H}_{E}(t)\) is the interaction Hamiltonian with an external field. The bath Hamiltonian \(\hat{H}_{B}\) is expressed as \begin{equation} \hat{H}_{B} = \sum_{j} \left[\frac{\hat{p}_{j}^{2}}{2 m_{j}} + \frac{1}{2} m_{j} \omega_{j}^{2} \hat{x}_{j}^{2} \right], \end{equation} (2) where \(\hat{p}_{j}\), \(\hat{x}_{j}\), \(m_{j}\), and \(\omega_{j}\) are the momentum, position, mass and frequency of the jth bath oscillator, respectively. The SB interaction \(\hat{H}_{I}\) is given by \(\hat{H}_{I} =\hat{V}\sum_{j} g_{j}\hat{x}_{j}\), where \(\hat{V}\) is the system part of the interaction, and \(g_{j}\) is the coupling constant between the system and the jth bath oscillator. The effect of the bath is characterized by the noise correlation function, \(C(t)\equiv \langle\hat{X}(t)\hat{X}(0)\rangle_{B}\), where \(\hat{X}\equiv \sum_{j} g_{j}\hat{x}_{j}\), and \(\langle\cdots \rangle_{B}\) represents the average taken with respect to the canonical density operator of the bath. The noise correlation function is expressed as \begin{equation} C(t) = \hbar \int_{0}^{\infty} \frac{d \omega}{\pi} J(\omega) \left[\coth \left(\frac{1}{2}\beta \hbar \omega \right) \cos(\omega t) - i \sin(\omega t) \right], \end{equation} (3) where \(J(\omega)\equiv \sum_{j} (\pi g_{j}^{2}/2m_{j}\omega_{j})\delta(\omega-\omega_{j})\) is the spectral density and β is the inverse temperature of the bath.

When we apply the SB model to problems of thermodynamics, because the main system is microscopic and because the quantum coherence between the system and bath characterizes the quantum nature of the system dynamics, the role of the SB interaction has to be examined carefully. For example, although the factorized thermal equilibrium state, \(\hat{\rho}_{\text{tot}}^{\text{eq}}=\hat{\rho}_{A}^{\text{eq}}\otimes \hat{\rho}_{B}^{\text{eq}}\), where \(\hat{\rho}_{A}^{\text{eq}}\) is the equilibrium state of the system without the SB interaction, is often employed as an initial state when investigating open quantum dynamics, in actual situations, the system and bath are quantum mechanically entangled (a phenomenon referred to as “bath entanglement”).26,27)

In Refs. 24 and 25, we presented a scheme for calculating thermodynamic variables in the SB model on the basis of simulations including an external perturbation using the hierarchical equations of motion (HEOM).2632) The key quantity in this investigation is the change of the “quasi-static Helmholtz energy” at time τ, which is defined as25) \begin{equation} \Delta F_{A}(\tau) \equiv \int_{0}^{\tau} \mathop{\text{tr}}\nolimits_{A} \left\{\hat{\rho}_{A}^{\text{qeq}} (t) \frac{\partial}{\partial t} \hat{H}_{A}(t) \right\}\,dt, \end{equation} (4) where \(\hat{\rho}_{A}^{\text{qeq}} (t)\) is the “quasi-static” reduced density operator. Here, note that, as we demonstrated numerically, when \(\hat{H}_{A} (t)\) changes much more slowly than the relaxation time of the system, \(\hat{\rho}_{A} (t)\) can be evaluated within the HEOM approach as the quasi-thermal equilibrium state of the system, \(\hat{\rho}_{A}^{\text{qeq}} (\tau)\approx \mathop{\text{tr}}\nolimits_{B} \{e^{-\beta(\hat{H}_{A}(\tau) +\hat{H}_{I} +\hat{H}_{B})}\}/Z_{\text{tot}} (\tau)\), where \(Z_{\text{tot}} (\tau)\equiv \mathop{\text{tr}}\nolimits_{A+B} \{e^{-\beta(\hat{H}_{A}(\tau) +\hat{H}_{I} +\hat{H}_{B})}\}\) at time \(\tau=t\). Then the change of the “quasi-static Boltzmann entropy”, \(\Delta S_{A}(\tau)\), is given by \begin{equation} \Delta S_{A}(\tau) = k_{\text{B}} \beta^{2} \frac{\partial}{\partial \beta} \Delta F_{A}(\tau). \end{equation} (5) Although the increase of the internal energy and the Boltzmann entropy of the system arise not only from the change in the system Hamiltonian itself but also from the change in the system part of the SB interaction, the HEOM allows us to evaluate both variables accurately. Using the HEOM, we can also evaluate the change of the bath part of the SB interaction energy and bath energy, while the bath energy itself is treated as infinitely large, because the bath is regarded as possessing an infinitely large heat capacity. With this treatment, we found that the SB model describes behavior that is consistent with the first law of thermodynamics. Explicitly, we found that the relation, \(\langle W(\tau)\rangle =\Delta U_{A}(\tau) -\Delta Q(\tau)\), is satisfied, where \(\Delta U_{A}(\tau)=\Delta F_{A}(\tau)+ T\Delta S_{A}(\tau)\) is the internal energy of the system and \(\Delta Q(\tau)\) is the heat released from the bath.25) Moreover, we have numerically confirmed that the total entropy production is alway positive. With these results strongly supporting the validity of our approach, in this paper, we evaluate \(\langle\exp[-\beta W(\tau)]\rangle\) using the HEOM formalism with the goal of extending the Jarzynski equality to the fully quantum regime characterized by a non-Markovian and non-perturbative SB interaction.

In what follows, we examine two definitions of the WCF: (i) a definition based on the partition function (PF) (the PF-WCF)33) and (ii) a definition based on trajectory (path) (the path-WCF). For an isolated quantum system, the PF-WCF has been defined as34) \begin{align} \langle \exp[-\beta W(\tau)] \rangle_{\text{PF}} &\equiv \mathit{tr}\{\mathrm{e}^{-\beta \hat{H}^{(H)} (\tau)} \mathrm{e}^{\beta \hat{H} (0)}\hat{\rho}_{\text{tot}}(0)\}\notag\\ &= \frac{Z_{\text{tot}}(\tau)}{Z_{\text{tot}}(0)}, \end{align} (6) where \(\hat{\rho}_{\text{tot}}(0)\equiv \hat{\rho}_{\text{tot}}^{\text{eq}}\), and \(\hat{H}^{(H)}(\tau)\) is the Heisenberg operator of \(\hat{H}(\tau)\). We rewrite this expression in terms of the time-reversal Liouville operator as \(\langle\exp[-\beta W(\tau)]\rangle_{\text{PF}}=\mathop{\text{tr}}\nolimits \{\hat{P}_{\text{PF}}(\tau)\}\), where \(\hat{P}_{\text{PF}} (\tau) =\exp_{+} [-i\int_{\tau}^{0} dt\,\hat{H}^{\times}(t)/\hbar]\hat{Z}_{\text{tot}}(\tau)\), with \(\hat{Z}_{\text{tot}}(\tau)\equiv\exp[-\beta\hat{H}(\tau)]\), and we have introduced the time-ordered exponential \(\exp_{\pm}\). Here and hereafter we use the hyperoperator notation \(\hat{\mathcal{O}}^{\times}\hat{f}\equiv \hat{\mathcal{O}}\hat{f}-\hat{f}\hat{\mathcal{O}}\) and \(\hat{\mathcal{O}}^{\circ}\hat{f}\equiv \hat{\mathcal{O}}\hat{f}+\hat{f}\hat{\mathcal{O}}\) for any operator \(\hat{\mathcal{O}}\) and operand operator \(\hat{f}\). Unfortunately, because the heat bath possesses infinitely many degrees of freedom, we cannot evaluate \(\hat{P}(\tau)\). Instead, because \(\langle\exp[-\beta W(\tau)]\rangle_{\text{PF}}=(Z_{A}(\tau)Z_{B}(\tau))/(Z_{A}(0) Z_{B}(0))\), where \(Z_{A}(\tau)\) and \(Z_{B}(\tau)\) are the system and bath parts of the partition functions, we can evaluate it indirectly using Eq. (4) as \(\langle\exp[-\beta W(\tau)]\rangle_{\text{PF}}\approx \exp [-\beta(\Delta U_{A}(\tau) -\Delta Q(\tau))]\) with \(\Delta U_{A}(\tau)=\Delta F_{A}(\tau)+ T\Delta S_{A}(\tau)\), which, off course, is the first law of thermodynamics.25)

Alternatively, we can use the path-WCF on the basis of the non-equilibrium trajectories (paths) as \(\langle\exp[-\beta W(\tau)]\rangle_{\text{path}}=\mathop{\text{tr}}\nolimits \{\hat{K}(\tau)\}\), where \(\hat{K}(\tau)\equiv \hat{U}(\tau, 0)\hat{P}_{\text{path}}(\tau)\hat{U}^{\dagger}(\tau, 0)\), with \begin{equation} \hat{P}_{\text{path}}(\tau) = \mathrm{e}_{+}^{-\frac{\beta}{2} \int_{0}^{\tau} dt\,\dot{\hat{H}}_{A}^{(H)}(t)}\hat{\rho}_{\text{tot}}(0) \mathrm{e}_{-}^{-\frac{\beta}{2} \int_{0}^{\tau} dt\,\dot{\hat{H}}_{A}^{(H)} (t)}, \end{equation} (7) and \(\hat{U}(\tau, 0)\equiv \exp_{+} [-(i/\hbar)\int_{0}^{\tau} dt\,\hat{H} (t)]\). Here, instead of \(\hat{P}(\tau) =\exp_{+} [-\beta\int_{0}^{\tau} dt\,\dot{\hat{H}}_{A}^{(H)} (t)]\hat{\rho}_{\text{tot}}(0)\),34) we use the symmetric form in Eq. (7,), because otherwise \(\hat{P}(\tau)\) is not Hermitian and, as a consequence, the WCF may not be real valued. The equation of motion for \(\hat{K}(\tau)\) is given by \begin{equation} \frac{\partial}{\partial \tau} \hat{K}(\tau) = - \left[\frac{i}{\hbar}\hat{H}_{\text{tot}}^{\times}(\tau) + \frac{\beta}{2} \dot{\hat{H}}_{A}^{\circ}(\tau) \right] \hat{K}(\tau). \end{equation} (8) The path (or functional) integral form of \(\hat{K}(\tau)\) is expressed as \begin{align} K(\xi, \xi',\tau) &= \int \int d\xi_{0}\,d\xi_{0}' \int_{\xi(0)=\xi_{0}}^{\xi(\tau)=\xi} \mathcal{D}[\xi(t)] \int_{\xi'(0)=\xi'_{0}}^{\xi'(\tau)=\xi'} \mathcal{D}[\xi'(t)]\notag\\ &\quad \times \exp\left[\frac{i}{\hbar}S_{\text{tot}}[\xi; \tau]-\frac{\beta}{2}W[\xi;\tau] \right] \rho_{\text{tot}} (\xi_{0}, \xi_{0}', t_{0})\notag\\ &\quad \times \exp\left[-\frac{i}{\hbar}S_{\text{tot}}[\xi'; \tau]-\frac{\beta}{2}W[\xi';\tau]\right], \end{align} (9) where \(S_{\text{tot}}[\xi;\tau]\) and \(W[\xi;\tau]\equiv \int_{0}^{\tau} dt\,\dot{H}_{A}(\xi;t)\) are the total action and the work as a functional of \(\xi(t)=(\sigma(t), x_{1} (t),x_{2} (t),\ldots)\), i.e., the system coordinate appended to the bath coordinate. The above expression implies that \(\langle\exp[-\beta W(\tau)]\rangle_{\text{path}}\) is obtained as an ensemble average of possible Liouville space pathways \(\xi(t)\) and \(\xi'(t)\). In the framework of the classical Jarzynski equality, the paths \(\xi(\tau)\) and \(\xi'(\tau)\) in the work operators \(-\beta W[\xi;\tau]/2\) and \(-\beta W[\xi';\tau]/2\) are determined as the paths of minimal action for \(-iS_{\text{tot}}[\xi';\tau]/\hbar\) and \(iS_{\text{tot}}[\xi;\tau]/\hbar\), respectively. In the present case, however, the paths \(\xi(\tau)\) and \(\xi'(\tau)\) are altered by the presence of \(-\beta W[\xi;\tau]/2\) and \(-\beta W[\xi';\tau]/2\) in Eq. (9,): This violates the condition to satisfy the Jarzynski equality. Nevertheless, we use Eq. (7), because it is natural to assume that the measurement of the work cannot be carried out without disturbing the dynamics of the main system, because it is regarded as small.

For an open quantum system, we can derive the HEOM for Eq. (9,) using the same procedure as that used to obtain the HEOM for Eq. (1,),2732) because the only difference between the quantum Liouville equation and Eq. (8,) is the presence of the work operator \(-\beta \dot{\hat{H}}_{A}^{\circ}(\tau)/2\) in the latter. We assume that the spectral density is given by the Drude distribution, \(J (\omega) =\eta \gamma^{2}\omega/(\omega^{2} +\gamma^{2})\), where η is the SB coupling strength, and γ is the inverse correlation time of the bath-induced noise. Then, the noise correlation function takes the form of a linear combination of exponential functions and a delta function: \(C(t) =\sum_{k=0}^{L} (c'_{k} + ic''_{k})\gamma_{k} e^{-\gamma_{k} t} + 2\Delta_{L}\delta (t)\), where \(c'_{k}\), \(c''_{k}\), \(\gamma_{k}\), and \(\Delta_{L}\) are constants. Then the HEOM for Eq. (8,) is expressed as \begin{align} &\frac{\partial}{\partial t} \hat{K}_{(n_{0},\ldots,n_{L})}(t)\notag\\ &\quad = - \left[\frac{i}{\hbar} \hat{H}_{A}^{\times}(t)+\frac{\beta}{2} \dot{\hat{H}}_{A}^{\circ}(t) + \Delta_{L} \hat{\Phi}^{2} + \sum_{k=0}^{L} n_{k} \gamma_{k}\right] \hat{K}_{(n_{0},\ldots,n_{L})}(t)\notag\\ &\qquad + \sum_{k=0}^{L} n_{k} \hat{\Theta}_{k} \hat{K}_{(\ldots,n_{k} - e_{k},\ldots)}(t) + \sum_{k=0}^{L} \hat{\Phi} \hat{K}_{(\ldots,n_{k} + e_{k},\ldots)}(t), \end{align} (10) where \(\mathbf{e}_{k}\) is the unit vector along the kth direction. Each hierarchical density matrix is specified by the index \(\mathbf{n}= (n_{0},\ldots,n_{L})\). The density matrix for \(\mathbf{n}= 0\) corresponds to the actual work distribution operator, \(\hat{K} (\tau)\). The initial state is prepared by numerically solving Eq. (10,) with the fixed Hamiltonian \(\hat{H}_{A} (t=0)\) until all of the hierarchy elements reach a steady state and then these elements are used as the initial state. Because Eq. (10) is identical to the HEOM in the case that \(\partial\hat{H}_{A}(t)/\partial t\) has no explicit time dependence, the steady–state solution of the first hierarchy element is identical to the correlated thermal equilibrium state defined by \(\hat{K}_{(\mathbf{n}= 0)}^{\text{eq}} =\mathop{\text{tr}}\nolimits_{B} \{\exp(-\beta\hat{H} (0))\}/{\mathop{\text{tr}}\nolimits_{A+B}} \{\exp(-\beta\hat{H} (0))\}\).

We now report the results of our numerical computations of the PF-WCF, \(\langle\exp[-\beta W(\tau)]\rangle_{\text{PF}}\), defined in Eq. (6,), and the path-WCF, \(\langle\exp[-\beta W(\tau)]\rangle_{\text{path}}\), defined in Eq. (7,) [or Eq. (9,)], under the periodic external force described by \(\hat{H}_{E}(t) =\hbar \omega_{0}\theta(t)\sin(\Omega t) (|g\rangle\langle e|+|e\rangle\langle g|)/4\), where \(\theta(t)\) is the step function and Ω is the frequency of the external field. The SB interaction is defined as \(\hat{V} =|g\rangle\langle e|+|e\rangle\langle g|\). While \(\langle\exp[-\beta W(\tau)]\rangle_{\text{path}}\) is evaluated using Eq. (10,), the procedures for computing \(\Delta F_{A}(\tau)\) and \(\langle W(\tau)\rangle\) [\(=-{\ln}(\langle\exp[-\beta W(\tau)]\rangle_{\text{PF}})/\beta\)] on the basis of the HEOM are explained in Ref. 25. The effect of the bath on thermodynamic properties in this SB model is characterized by the SB coupling strength, the bath temperature, and the noise correlation time.24,25,35) Throughout this investigation, we fix \(\beta\hbar \omega_{0} = 1\) and \(\gamma=\omega_{0}\). These values correspond to intermediate temperature and moderately non-Markovian noise.

In Fig. 1(a), we display the time dependences of the PF-WCF, path-WCF, and \(\Delta F_{A} (\tau)\) for several values of the excitation frequency, Ω, in the weak (\(\eta=0.1\)) SB coupling case. Due to the production of heat, \(\Delta Q(\tau)\), the PF-WCF, i.e., \(-{\ln}(\langle\exp[-\beta W(\tau)]\rangle_{\text{PF}})/\beta=\langle W(\tau)\rangle\), increases as a function of time in accordance with the first law of thermodynamics, \(\langle W(\tau)\rangle =\Delta F_{A}(\tau)+ T\Delta S_{A}(\tau) -\Delta Q(\tau)\). This increase is largest in the resonant case \(\Omega =\omega_{0}\), because the excitation of the system is most efficient there. We calculated \(\Delta Q(\tau)\) separately and found that the cycle of the production of heat is similar to that of the WCF. Because the entropy production \(\Sigma_{\text{tot}}(\tau) =\Delta S_{A}(\tau) -\Delta Q(\tau)/T\) exhibits a time lag after the external excitation, we observe a phase delay in the PF-WCF results due to this contribution. The delay is largest at the resonant excitation, \(\Omega=\omega_{0}\), because the entropy production is largest there. The amplitudes of oscillations are suppressed because \(-\Delta Q(\tau)\) partially cancels out the contribution from the free energy. For the slowest modulation, \(\Omega=0.1\omega_{0}\), in this weak coupling case, the free energy is almost canceled by the heat production. In this case, the time evolution of the PF-WCF is dominated by \(T\Delta S_{A}(\tau)\).

Figure 1. (Color online) The quantity \(-{\ln}\langle\exp[-\beta W(\tau)]\rangle/\beta\) evaluated as the path-WCF (solid curves) and the PF-WCF (dashed curves), and the change in the free energy, \(\Delta F_{A}(\tau)\) (black dots), under the external perturbation \(\hat{H}_{E}(t) =(\theta(t)\sin(\Omega t)/4)\hbar \omega_{0} (|g\rangle\langle e|+|e\rangle\langle g|)\) are plotted as functions of time for fixed \(\beta\hbar \omega_{0} = 1\) in (a) the weak (\(\eta=0.1\)), (b) intermediate (\(\eta=1\)), and (c) strong (\(\eta=3\)) SB coupling cases. The colored dashed and solid curves represent the results for different frequencies: \(\Omega=0.1\omega_{0}\) (red curves), \(\Omega=\omega_{0}\) (green curves), and \(\Omega=5\omega_{0}\) (blue curves).

While the time evolution of the PF-WCF differs significantly from \(\Delta F_{A} (\tau)\), that of the path-WCF is quite similar. This similarity can be understood as follows. First, note that the ensemble average of \(W[\xi;\tau]\equiv \int_{0}^{\tau} dt\,\dot{H}_{A}(\xi;t)\) in Eq. (9,) is taken after the time integration of \({\dot{H}}_{A}(\xi;t)\) for a given Liouville path \(\xi(t)\). Then, because the contribution of \(\dot{H}_{A}(\xi;t)\) oscillates between positive and negative values rapidly in time, the heat production involved in the definition in Eq. (9) is suppressed. By contrast, the ensemble average of the PF-WCF is taken step by step in the time integration. Thus, \(-{\ln}(\langle\exp[-\beta W(\tau)]\rangle_{\text{PF}})/\beta=\langle W(\tau)\rangle\) becomes thermodynamic process function, while \(\langle\exp[-\beta W(\tau)]\rangle_{\text{path}}\) is not.

In Figs. 1(b) and 1(c), we present the results for the intermediate and strong SB coupling cases. As pointed out previously,35,36) the efficiency of the heat current is suppressed when the SB coupling becomes strong. Because the effective SB coupling depends on the characteristic time scale of the system, the increase of the PF-WCF is suppressed in the case \(\Omega\geq \omega_{0}\), while it is enhanced in the case \(\Omega <\omega_{0}\).37) In Fig. 1(b), due to the effect of the moderately strong SB coupling, the system closely follows its instantaneous equilibrium state, as the entropy production, \(\Sigma_{\text{tot}}(\tau)\), is suppressed. As a result, the amplitudes of oscillations and phase delay are small. In Figs. 1(b) and 1(c), because the eigenenergies of the system are significantly altered by the strong system–bath coupling, the deviation of the time profile of the PF-WCF from \(\Delta F_{A} (\tau)\) increases as Ω decreases.

For the path-WCF results represented by the solid curves, the deviation becomes larger as the strength of the SB coupling increases. This is because the calculated \(\Delta F_{A} (\tau)\) involves a contribution from the energy of the system part of the SB interaction,25) whereas \({\dot{H}}_{A}(\xi;t)\) does not include the contribution from \(H_{I}\), which is also time dependent in the reduced description, due to the non-Markovian nature of the noise that arises from the bath. The path-WCF approaches the free energy when the characteristic time scale of the system is shorter than \(1/\gamma \), because in this case, the contribution from the SB interaction is insignificant.

In the weak coupling case considered in Fig. 1(a), the difference between the path-WCF results and the free energy results is large in the resonant excitation case, \(\Omega=\omega_{0}\). This difference has a non-kinetic origin arising from the alteration of the Liouville path in Eq. (9) due to the presence of the work functional \(-\beta W[\xi;\tau]\equiv -\beta\int_{0}^{\tau} dt\,\dot{H}_{A}(\xi;t)\). In the case \(\Omega =\omega_{0}\), the contribution of the external perturbation, \({\dot{H}}_{A}(\xi;t)\), in the work operator is large, because the cyclic excitation with this excitation strongly perturbs the system dynamics. Thus the paths that should be determined by the total action are altered, and as a result, the calculated path-WCF exhibits time profiles that differ significantly from those in other cases. For this reason, the path-WCF results also differ significantly from \(\Delta F_{A} (\tau)\) in the low temperature case, because the contribution of \(-\beta W[\xi;\tau]\) is larger there (results not shown).

In this paper, we demonstrated a method for extending the Jarzynski equality to the fully quantum regime. We evaluated the WCF defined in two ways, the PF-WCF and path-WCF, using the numerally rigorous HEOM formalism. Although the path-WCF agrees with the free energy reasonably well, in particular in a weak SB coupling case or the fast excitation cases, while the PF-WCF exhibits very different time-dependence due to the heat production, the result is not equality but approximation. This discrepancy arises from the contribution of the SB interaction, which should also play a role in the classical case if the SB coupling strength is comparable to the system energy. Indeed, if we employ quantum hierarchical Fokker–Planck equations (QHFPEs) for a system described by Wigner distribution functions, we can investigate not only the quantum case but also the classical case by taking the classical limit: We can easily identify purely quantum mechanical effects by comparing the classical and quantum results for the Wigner distribution.37,38)

It should be mention that, although here we introduced the path-WCF, this is not physical observable,21,34) as seen from Eq. (8,). Moreover, we cannot determined the paths in the functional formalism of Eq. (9,), due to the limitation introduced by the uncertainty principle. Thus, in order to evaluate the free energy in the fully quantum regime, the path-WCF is not practical. Instead, Eq. (4) should be used to evaluate the free energy.

Although the present investigation is limited to spin-Boson systems for the specific definitions of the WCF, the applicability of our approach based on the HEOM formalism is in fact more general. Indeed, the same approach can be applied to all of the systems to which the HEOM formalism has been previously applied.26,27) Different definitions of the WCF should also be examined. We leave such extensions to future studies to be carried out in the context of the fluctuation theorem.