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

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 (2020) 234107], 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)", $\langle \exp[-\beta W] \rangle$, where $W$ is the work performed on the system and $\beta$ 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.

the system cannot reach thermal equilibrium on its own without a system-bath interaction: We cannot assume a canonical distribution as the equilibrium state for either the system of the heat-bath, due to the presence of the system-bath interaction. In the present paper, we numerically evaluate the work characteristic function (WCF), exp[−βW(τ)] , and ∆F A (τ) 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 system-bath (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 bŷ whereĤ A (t),Ĥ I andĤ B are the Hamiltonian of the system, the interaction and the bath, respectively. The system Hamiltonian is given byĤ A (t) =Ĥ 0 A +Ĥ E (t), withĤ 0 A = 1 2 ω 0 (|e e|− |g g|) andĤ E (t) = 0 for t ≤ 0, where |e and |g are the excited and ground states of the system, andĤ E (t) is the interaction Hamiltonian with an external field. The bath Hamiltonian H B is expressed asĤ wherep j ,x j , m j and ω j are the momentum, position, mass and frequency of the jth bath oscillator, respectively. The SB interactionĤ I is given byĤ I =V j g jx j , whereV 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, whereX ≡ j g jx j , and . . . B represents the average taken with respect to the canonical density operator of the bath. The noise correlation function is expressed as where J(ω) ≡ j πg 2 j /2m j ω j δ(ω − ω 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,ρ eq tot =ρ eq A ⊗ρ eq B , 2/11 whereρ eq A 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). [26][27][28][29][30][31][32] The key quantity in this investigation is the change of the "quasi-static Helmholtz energy" at time τ, which is defined as 25) whereρ qeq A (t) is the "quasi-static" reduced density operator. Here, note that, as we demonstrated numerically, whenĤ A (t) changes much more slowly than the relaxation time of the system,ρ A (t) can be evaluated within the HEOM approach as the quasi-thermal equilibrium state of the system,ρ qeq Then the change of the "quasi-static Boltzmann entropy", ∆S A (τ), is given by 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 these energies themselves are 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, is the internal energy of the system and ∆Q(τ) 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 exp[−βW(τ)] 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 3/11 path-WCF). For an isolated quantum system, the PF-WCF has been defined as 34) whereρ tot (0) ≡ρ eq tot , andĤ (H) (τ) is the Heisenberg operator ofĤ(τ). We rewrite this expression in terms of the time-reversal Liouville operator as exp[−βW(τ)] PF = tr P PF (τ) , , and we have introduced the time-ordered exponential exp ± . Here and hereafter we use the hyperoperator notationÔ ×f ≡Ôf −fÔ andÔ •f ≡Ôf +fÔ for any operatorÔ and operand operatorf . Unfortunately, because the heat bath possesses infinitely many degrees of freedom, where Z A (τ) and Z B (τ) are the system and bath parts of the partition functions, we can , 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 andÛ(τ, 0) 34) we use the symmetric form in Eq. (7), because other-wiseP(τ) is not Hermitian and, as a consequence, the WCF may not be real valued. The equation of motion forK(τ) is given by The path (or functional) integral form ofK(τ) is expressed as where S tot [ξ; τ] and W[ξ; τ] ≡ τ 0 dtḢ A (ξ; t) are the total action and the work as a functional of ξ(t) = (σ(t), x 1 (t), x 2 (t), . . .), i.e., the system coordinate appended to the bath coordinate. 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), [27][28][29][30][31][32] because the only difference between the quantum Liouville equation and Eq. (8) is the presence of the work operator −βḢ • A (τ)/2 in the latter. We assume that the spectral density is given by the Drude distribution, J(ω) = ηγ 2 ω/(ω 2 + γ 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: where c k , c k , γ k , and ∆ L are constants. Then the HEOM for Eq.(8) is expressed as where e k is the unit vector along the kth direction. Each hierarchical density matrix is specified by the index n = (n 0 , . . . , n L ). The density matrix for n = 0 corresponds to the actual work distribution operator,K(τ). The initial state is prepared by numerically solving Eq.  24,25,35) Throughout this investigation, we fix β ω 0 = 1 and γ = ω 0 .
These values correspond to intermediate temperature and moderately non-Markovian noise. 6/11 In Fig. 1 (a) 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 Ω ≥ ω 0 , while it is enhanced in the case Ω < ω 0 . 37) In addition, due to the effect of the moderately strong SB coupling, the system closely follows its instantaneous equilibrium state, as the entropy production, Σ tot (τ), is suppressed. As a result, the amplitudes of oscillations and phase delay are small. In Figs. 1 (b) and (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 ∆F A (τ) 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 ∆F A (τ) involves a 7/11 contribution from the energy of the system part of the SB interaction, 25) whereasḢ A (ξ; 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 t c , 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, Ω = ω 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 −βW[ξ; τ] ≡ −β τ 0 dtḢ A (ξ; t). In the case Ω = ω 0 , the contribution of the external perturbation,Ḣ A (ξ; 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 ∆F A (τ) in the low temperature case, because the contribution of −βW[ξ; τ] 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, 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 8/11 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. 9/11