Quantum thermodynamics from the nonequilibrium dynamics of open systems: energy, heat capacity and the third law

TL;DR

This paper investigates the conditions under which quantum thermodynamics is well-defined for open quantum systems, using exact solutions of quantum Brownian motion to analyze energy, heat capacity, and the Third Law, challenging previous claims about negative specific heat and entanglement.

Contribution

It provides a systematic open system approach to quantum thermodynamics, deriving conditions for thermodynamic laws and analyzing a solvable model to clarify misconceptions about heat capacity and the Third Law.

Findings

Internal energy is extensive in the studied models.

Heat capacity remains positive down to zero temperature.

The Third Law holds in the examined quantum systems.

Abstract

We take the perspective of open quantum systems and examine from their nonequilibrium dynamics the conditions when the physical quantities, their relations and the laws of thermodynamics become well defined and viable for quantum many body systems. We first describe how an open system nonequilibrium dynamics (ONEq) approach is different from the closed combined system + environment in a global thermal state (CGTs) setup. Only after the open system equilibrates will it be amenable to conventional thermodynamics descriptions, thus quantum thermodynamics (QTD) comes at the end rather than assumed in the beginning. We see the open system approach having the advantage of dealing with nonequilibrium processes as many experiments in the near future will call for. Because it spells out the conditions of QTD's existence it can also aid us in addressing the basic issues in quantum thermodynamics from first principles in a systematic way. We then study one broad class of open quantum systems where the full nonequilibrium dynamics can be solved exactly, that of the quantum Brownian motion of $N$ strongly coupled harmonic oscillators, interacting strongly with a scalar field environment. Here we focus on the internal energy, heat capacity and the Third Law. We show for this class of physical models, amongst other findings, the extensive property of the internal energy, the positivity of the heat capacity and the validity of the Third Law from the perspective of the behavior of the heat capacity toward zero temperature. These conclusions obtained from exact solutions and quantitative analysis clearly disprove claims of negative specific heat in such systems and dispel allegations that in such systems the validity of the Third Law of thermodynamics relies on quantum entanglement.

Figures

Equations

$$S$$

$$+\int d^{4}x\frac{1}{2}{\partial }_{\mu }\varphi (x){\partial }^{\mu }\varphi (x),$$

$${\rho }_i$$

$$U(t_{f},t_{i})=\operatorname{T}\,\exp\biggl{[}-i\int_{t_{i}}^{t_{f}}\!ds\;H(s)\biggr{]}\,,$$

$$\rho^{(\chi)}(t_{f})=\operatorname{Tr}_{\phi}\Bigl{\{}U(t_{f},t_{i})\,\rho(t_{i})\,U^{\dagger}(t_{f},t_{i})\Bigr{\}}\,,$$

$$\langle O^{(\chi)}\rangle=\operatorname{Tr}_{\chi}\Bigl{\{}\rho^{(\chi)}(t_{f})O^{(\chi)}\Bigr{\}}\,,$$

$${\rho }^{(\chi )}({\chi }_f,{\chi }_f^{\prime },t_f)$$

$$\displaystyle\qquad\qquad\qquad\int_{\chi_{i}}^{\chi_{f}}\!\mathcal{D}\chi_{+}\int_{\chi^{\prime}_{i}}^{\chi^{\prime}_{f}}\!\mathcal{D}\chi_{-}\;\exp\Bigl{\{}i\,S^{(\chi)}[\chi_{+}]-i\,S^{(\chi)}[\chi_{i}]\Bigr{\}}\mathcal{F}[\chi_{+},\chi_{-}]\,,$$

$$m{\ddot{\chi }}_i(t)+m{\omega }_b^2{\chi }_i(t)+\sum _{j\ne i}m\sigma {\chi }_j(t)-e^2{\int }_0^{t}ds\sum _j{G}_R^{(\varphi )}(t-s,{\mathbf{z}}_i-{\mathbf{z}}_j){\chi }_j(s)=e{\xi }_i(t).$$

$$⟨{\xi }_i(t)⟩$$

$$G_R^{(\varphi )}(x-x^{\prime })$$

$$G_H^{(\varphi )}(x-x^{\prime })$$

$${\tilde{G}}_H^{(\varphi )}(\omega ;r)=\coth \frac{\beta \omega }{2}\mathrm{Im}{\tilde{G}}_R^{(\varphi )}(\omega ;r),$$

$$\tilde{f}(\omega )$$

$$\mathbf{\Xi }$$

$$\ddot{\mathbf{\Xi }}(t)+{\mathbf{\Omega }}_b^2\cdot \mathbf{\Xi }(t)-\frac{e^2}{m}{\int }_0^{t}ds{\mathbf{G}}_R^{(\varphi )}(t-s)\cdot \mathbf{\Xi }(s)=\frac{e}{m}\mathbf{\xi }(t).$$

$$\mathbf{\Xi }(t)={\mathbf{d}}_1(t)\cdot \mathbf{\Xi }(0)+{\mathbf{d}}_2(t)\cdot \dot{\mathbf{\Xi }}(0)+\frac{e}{m}{\int }_0^{t}ds{\mathbf{d}}_2(t-s)\cdot \mathbf{\xi }(s),$$

$$\tilde{\mathbf{d}}_{2}(\omega)=\Bigl{[}\boldsymbol{\Omega}_{b}^{2}-\omega^{2}\mathbf{I}-\frac{e^{2}}{m}\,\tilde{\mathbf{G}}_{R}^{(\phi)}(\omega)\Bigr{]}^{-1}\,.$$

$${\int }_{-\infty }^{\infty }dt|{\mathbf{d}}_i(t)|\nless \infty ,$$

$${\mathbf{\sigma }}_{\chi \chi }(t)=\frac{1}{2}⟨\{\mathbf{\Xi }(t),{\mathbf{\Xi }}^T(t)\}⟩$$

$$\underset{t\to \infty }{\lim }{\mathbf{\sigma }}_{\chi \chi }(t)$$

$${\tilde{\mathbf{d}}}_2(\omega )-{\tilde{\mathbf{d}}}_2^{†}(\omega )=2i\mathrm{Im}{\tilde{\mathbf{d}}}_2(\omega )$$

$$\mathbf{A}^{-1}-\mathbf{B}^{-1}=\mathbf{A}^{-1}\cdot\bigl{(}\mathbf{B}-\mathbf{A}\bigr{)}\cdot\mathbf{B}^{-1}\,,$$

$${\mathbf{\sigma }}_{\chi \chi }(\infty )=\frac{1}{m}{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\coth \frac{\beta \omega }{2}\mathrm{Im}{\tilde{\mathbf{d}}}_2(\omega )=\mathrm{Im}{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\coth \frac{\beta \omega }{2}{\tilde{\mathbf{G}}}_R^{(\chi )}(\omega ),$$

$${\tilde{\mathbf{G}}}_R^{(\chi )}(\omega )=\frac{1}{m}{\tilde{\mathbf{d}}}_2(\omega ).$$

$${\mathbf{\sigma }}_{\upsilon \upsilon }(t)=\frac{1}{2}⟨\{\dot{\mathbf{\Xi }}(t),{\dot{\mathbf{\Xi }}}^T(t)\}⟩,$$

$${\mathbf{\sigma }}_{\upsilon \upsilon }(\infty )=\mathrm{Im}{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }{\omega }^2\coth \frac{\beta \omega }{2}{\tilde{\mathbf{G}}}_R^{(\chi )}(\omega ).$$

$$E(t)$$

$$\displaystyle=\frac{m}{2}\operatorname{Tr}\,\Bigl{\{}\boldsymbol{\sigma}_{\upsilon\upsilon}(t)+\boldsymbol{\Omega}_{p}^{2}\cdot\boldsymbol{\sigma}_{\chi\chi}(t)\Bigr{\}}\,.$$

$${\omega }_p^2-{\omega }_b^2=-\frac{4\gamma \Lambda }{\pi },$$

$$\displaystyle E(\infty)=\frac{m}{2}\operatorname{Tr}\,\Bigl{\{}\boldsymbol{\sigma}_{\upsilon\upsilon}(\infty)+\boldsymbol{\Omega}_{p}^{2}\cdot\boldsymbol{\sigma}_{\chi\chi}(\infty)\Bigr{\}}$$

Full text

I Introduction

In this series of papers we take the perspective of open quantum systems (OQS) and examine from their nonequilibrium (NEq) dynamics the conditions when the physical quantities, concepts, constructs and the time-honored laws of thermodynamics (TD) become well defined and viable for quantum many body systems. We utilize one broad class of models where the nonequilibrium dynamics can be solved exactly – the Brownian motion of strongly coupled (SC) harmonic oscillators, interacting strongly with a scalar field environment – to explore a range of basic issues in quantum thermodynamics (QTD). The exact solutions possible in these OQSs enable us to examine and define these conditions more precisely in a quantitative, systematic and transparent way. This approach hopefully compensates for the rather loose, qualitative and at times contrived way thermodynamic descriptions for quantum systems are proposed because of the need to adhere to the dictum of classical thermodynamics, which is valid only under very special conditions.

A clarification in the meaning and contents of quantum thermodynamics (QTD) Mahler might be useful before we proceed: to us, it is the study of the thermodynamic properties of quantum many-body systems (MBS). Quantum now refers not just to the particle spin-statistics (boson vs fermion) aspects: the rather limited meaning of ‘quantum’ in traditional quantum statistical mechanics (QSM), but also includes in the present era the quantum phase aspects, such as quantum coherence, quantum correlations, and quantum entanglement. This is where quantum information has a hand in QTD ktln2 . The thermodynamics connotation can be extended to include systems not necessarily in equilibrium at all times, thus encompassing dissipative and relaxation processes for systems deviating from equilibrium, including linear or nonlinear response theories applied to quantum MBS 111Linear response theory considers small variations in the system while staying in thermal contact with the bath. This is the underlying assumption in the use of thermal Green’s functions, which is within the test-field approximation in quantum field theory terms. In a fully NEq treatment of the open system’s quantum dynamics, both the system and environment variables are dynamically determined. Thus it can cope with situations where the quantum system is small and the environment is finite., familiar in condensed matter or chemical physics and for classical MBS, topics incorporated in the traditional field of NEq TD NEqTDbooks . These considerations can be extended to weakly nonequilibrium conditions but not for far-from-equilibrium, fully arbitrary time evolutions. That is when an open quantum system treatment becomes necessary.

The issues addressed in this first paper encompass the nature of internal energy, heat capacity and the third law for a fully nonequilibrium (NEq) system. We demonstrate what it takes for it to evolve to an equilibrium (Eq) condition, and from that point establish the connection with traditional TD theory. The conditions for traditional TD theory to be well-defined and operative for a classical or quantum system are very specific despite its wide ranging applicability: A system of relatively fewer degrees of freedom in the presence of a thermal bath of a huge number or infinite degrees of freedom (we shall consider only heat but no particle transfer here and thus the TD refers only to canonical ensembles), the coupling between the system and the bath is vanishingly small, and the system is eternally in a thermal equilibrium state by proxy with the bath which is impervious to any change in the system 222This means that action of the system on the bath is excluded from TD considerations. In fact in TD the bath variables are not dynamical variables determined consistently by the interplay between the system and the bath through their coupled equations of motion, they only provide TD parameters such as temperature or chemical potential.. Already for classical systems, there is a difference between equilibration and thermalization. Equilibration refers to the system evolving to a steady state after relaxation. It is broader than thermalization, which refers to the system approaching a state described by the Boltzmann distribution. When the system-bath coupling is nonvanishing, such a difference is clearly discernible. For example, the potential of mean force Kirkwood is introduced to deal with such a situation. Details can be found in Appendix D.

For quantum systems this difference between equilibration and thermalization certainly remains, (see, e.g., GE16 ). New challenges at zero or very low temperatures posed by non-Markovian environments and in the treatment of non-Markovian dynamics can become prominent. By virtue of its ability to provide a first principles derivation of noise from quantum fluctuations (e.g., for Gaussian noise via the Feynman-Vernon identity, instead of being put in by hand), and linking fluctuations and noise with dissipation and relaxation by dynamical relations (such as the fluctuation-dissipation relation which can be traced to the unitarity in the original closed system before one coarse grains the environment to a description of mean field dynamics and its fluctuations), the open quantum system approach is also a natural setting for incorporating stochastic thermodynamics stochTD , which has seen a wide range of chemical and biological science applications 333Many physical systems show two intermediate stages between quantum and classical, namely, stochastic and semiclassical. Conventional stochastic thermodynamics starts from classical or macroscopic physics. Noise is added in phenomenologically for the consideration of fluctuations phenomena under different circumstances for specific purposes. Being rooted in classical physics conventional stochastic thermodynamics cannot capture the quantum features so easily. Open quantum systems approach on the other hand starts from microphysics at the quantum level. One can derive the stochastic equations including quantum and thermal noises: Langevin, Fokker-Planck or master equations for the description of fluctuations phenomena. Thus in the quantum open system approach the pathway from the quantum regime to the stochastic regime is well laid out. Taking the distributional average of noise yields the mean field theories at the semiclassical level. To go from quantum to classical physics one needs to add decoherence considerations, but the pathway is completely accessible. The challenge is, can we come up with an appropriate quantum microphysics model for the macroscopic phenomena of interest? .

The set-up: In our opinion, the NEq dynamics of open quantum systems, e.g., in the tradition of Feynman-Vernon, Caldeira-Leggett et al FeyVer ; CalLeg ; HPZ even though requiring more work, is the preferred setting for addressing new issues in quantum thermodynamics for future challenges 444Similar viewpoint has been expressed by a few others, notably, Kosloff Kosloff .. This is in comparison with a popular set-up which has been studied more in the literature namely, that of a global thermal state (CGTs) assumed for the combined or closed system (C) = system (S) + bath (B) 555The CGTs set up is used by many authors, notably GelTho ; KimMah10 ; HTL11 ; SIW15 ; Seifert16 ; PhiAnd ; GSI88 ; HIT08 ; Srednicki94 ; GLTZ06 ; PSW06 ; LPSW09 ; SF12 ; Reimann08 ; Rigol09 ; PSSV11 ; Deutsch91 ; CR10 ; GE16 .. In the CGTs set-up the initial and final states of C are the same, namely the combined system remains in an equilibrium global thermal state, because the dynamics of the combined closed system is unitary. This is visibly closest to the setting of thermodynamics and thus naturally convenient for exploring small extensions of thermodynamics. By contrast the open system NEq (ONEq) approach deals with time evolution of the open system. It requires the specification of the initial conditions and the derivation of the late time behavior of the open system. For those systems that upon interaction with a bath equilibrate at late times, one may then connect its behavior with the descriptions of thermodynamics. For sure this is a many-to-one relation - many different initial conditions can produce the same final steady or equilibrium state, or that there is no common final steady or equilibrium state. A lot depends on the structures of the system, the properties of the bath and the way they interact. All the above mentioned factors need to be considered for interacting quantum many-body systems before we construct thermodynamical quantities, address thermodynamical issues and invoke (or hasten to claim success in revoking 2LawSaga ) the well-established thermodynamical laws. We will elaborate on their differences in the following.

I.1 Main Contents

There are three main components in this paper:

1. Set-up and conditions: The physical differences of the set-ups, comparing our open system nonequilibrium (ONEq) approach (level 2) with traditional TD (level 0) on the one hand, and with the global thermal (CGTs) state set-up (level 1) on the other. In TD, as mentioned above, the system-bath coupling has to be vanishingly small whereas in both the Level 1 and 2 treatments the system-bath coupling can be strong. We will mention the CGTs approach as many existing works are based on this setup, but focus more on how to use an open system approach to define and quantify quantum thermodynamics. In a companion paper GTOS , we will attempt to build some bridges between these two approaches, via generating functional and reduced density matrix formulations. The hope is that from the open system perspective, one may be able to identify which entities and concepts are more suitable for treating new problems in QTD and which are residues of the old which may hinder new developments. Other authors using an open system approach to quantum thermodynamics include Duarte and Caldeira DuaCal who treated a coupled-oscillator system by the influence functional method, Carrega et al. CarrWeiss15 who treated a two level system via moment-generating functionals, and Esposito et al. EspositoPRL15 using nonequilibrium Green’s functions.

2. Model with exact solutions: We use a quantum Brownian motion (QBM) model of harmonic oscillators with strong-coupling both within the system ($\sigma$) and interacting with a scalar field bath ($\gamma$). The merit of this model, which represents a rather broad class of physical problems, is that being a Gaussian system it can yield exact solutions which enable us to cross examine the relevant issues, leaving little room for speculation. Even when familiar quantities like energy and entropy can be defined in different ways under different conditions, since we are treating NEq dynamics, if we make precise specific conditions, these quantities are defined. There is no worry about ambiguity. The results from this model study are used for addressing the following issues:

3. Issues and consequences:

Energy extensivity.

Thermodynamic functions are well defined under the conditions when thermodynamics theory is viable, namely, that the system is very weakly coupled to the bath, the bath being a passive source which provides a temperature parameter, not a dynamical variable which can back-react on the system. It is a meaningful question to ask if the nice properties we are accustomed to in conventional thermodynamics, e.g., the extensive property of internal energy, will still hold for strongly interacting quantum systems. In the model we studied here we answer this question in the affirmative, that the internal energy remains extensive under strong coupling.

Heat capacity.

From the internal energy we calculate the heat capacity and examine its behavior toward $T=0$. We find a power law, not an exponential decay. This has significant implications. This aids us to address a version of the third law and to resolve some puzzles raised in the literature such as the claimed negative specific heat near absolute zero even in well behaved systems HasegawaJMP .

Third law.

There are several formulations and statements of the third law. We approach it from the behavior of the heat capacity near absolute zero, which aids us to resolve some puzzles raised in the literature such as the claimed negative specific heat near absolute zero even in well behaved systems HasegawaJMP , and address some concerns expressed by Hanggi, Ingold, Talkner, Weiss, et al HanIng06 ; HIT08 ; IHT09 ; AIW14 ; SIW15 ; FordOC .

Vedral et al Vedral ; Vedral1 invoked heat capacity as an indicator of entanglement, and raised the issue of how the entanglement at a system’s ground state bears on the third law. For the (spin) system they studied they made the claim that “the validity of the third law of thermodynamics relies on quantum entanglement”. Using the behavior of the heat capacity at $T=0$ we derived here, combined with our earlier results on the entanglement between two coupled-oscillators interacting with a zero temperature bath HHPRD , we show that this is not the case at least for the coupled oscillator system. There is no connection between entanglement in the system and the third law.

The ONEq approach we adopt for the dynamics of the system provides means to calculate entropy production, but not before the meaning and definition of entropy for interacting quantum systems can be understood and clarified. We say this because even the most commonly invoked von Neumann entropy has problems if not used and understood properly. We shall mention this issue at the end of this paper but leave a proper treatment of heat, entropy, entanglement, and from it the first and second Law, to the second paper QTD2 in this series.

I.2 Closed-system Global Thermal State versus Open-system Evolved Equilibrium State

We begin by stating a few basic facts connecting the three levels of treatments: level 0 thermodynamics (TD), level 1 closed system (system and environment combined) in a global thermal state (CGTs) and level 2 open system evolving to an equilibrium state (ONEq).

1. Traditional statistical mechanics treats many-body systems in thermal (canonical distribution) and chemical (grand canonical) equilibrium. The starting point of quantum statistical mechanics (QSM) is probability density, no quantum phase information is invoked. This is encoded in the two fundamental postulates of quantum statistical mechanics: equal a priori probability to all accessible states and random phase approximation. Thus from a quantum information viewpoint, the system of interest to QSM is already fully decohered in the energy basis and behaves classically in an effective way – what is quantum in QSM only pertains to quantized energy levels and particle spin-statistics.

2. Partition function is well defined only for systems in thermal equilibrium. It *is ill-defined for systems under nonequilibrium conditions when the notion of temperature is lacking. * Pathologies may ensue if it is forced upon even perfectly normal systems (in contradistinction to systems for which the canonical ensemble does not exist and the heat capacity is negative in the microcanonical ensemble, such as gravitating systems). As noted in EspositoJSM , if one proceeds from assuming that the combined system + environment is in a thermal state the behavior of the heat capacity of the system is different when it is derived from the energy of the central system at equilibrium or from a partition function approach HanIng06 ; HIT08 . By examining the open system nonequilibrium dynamics with no reference to the partition function, one could avoid these pathologies. Likewise, old notions such as the Hamiltonian of mean force Kirkwood are only meaningful in the conceptual framework of equilibrium systems HTL11 as in CGTs.

3. The oft-heard statement, that the generating functional (in quantum field theory) is equal to the partition function (in equilibrium statistical mechanics) is true only for thermal fields, i.e., there exists a canonical distribution where a thermal state is well defined for all times. This statement arises from treating thermal (finite temperature) fields with imaginary (Matsubara) time quantum field theory. If one uses real time representation to describe the NEq dynamics of open systems the generating functional remains well defined but it is not the (canonical) partition function defined in imaginary time.

4. If an open system upon interaction with its environment can equilibrate at late times, and if it is further thermalized, it enters a thermal state. But this equilibrium state is different from that of a system in contact with a heat bath which behaves in a totally passive and non-dynamical way, in particular, with no back-action on the system. The latter is where a theory of quantum thermodynamics is often constructed, namely, from a simple extension of conventional classical thermodynamics. The difference lies in the dynamical correlations between the system and the bath, which conventional thermodynamics ignores completely by assuming a vanishingly small coupling.

5. There are important differences between the ONEq and the CGTs setups in their goals, approaches and consequences. There are also key differences between CGTs and thermalization in a closed quantum system in the vein of eigenstate thermalization hypothesis. Since the latter is an active topic in the last decade with many important contributions, we can only focus on the differences from the open system approach on the specific issues of interest to us here and cite some representative references and reviews for interested readers to appreciate the scope Srednicki94 ; GLTZ06 ; PSW06 ; LPSW09 ; SF12 ; Mahler ; Reimann08 ; Rigol09 ; PSSV11 ; Deutsch91 ; CR10 . We highlight some key features below. In Appendix C, we will illustrate some aspects of the ONEq and the CGTs setups with a simple model calculation from the ubiquitous QBM model.

Set-up and goals. In this work we assumed the field to be in a thermal state prior to its coupling to the system oscillators, which initially can be in an arbitrary state. Thus the system oscillators and the field are generically out of equilibrium before and after the interaction is turned on. Our focus is on the subsequent dynamics and relaxation of the system oscillators, without assuming the coupling to be weak. The ‘pure state quantum statistical mechanics’ assumes the whole system is in a pure state throughout. The main goal in GE16 is to derive statistical mechanics and thermodynamics from quantum mechanics without resorting to the notion of ensembles. It aims to show that even pure quantum states of interacting many-body systems can display relaxation to equilibrium and, in special cases, thermalize.

The methods developed in the ‘pure state quantum statistical mechanics’ literature are usually applied to closed systems without an intrinsic system-bath distinction. For instance, in Cramer et al. CE09 , a one-dimensional harmonic lattice is shown to locally relax to Gaussian states for arbitrary choice of subsystem and a wide class of initial states. The authors note the following “Every part of the system forms the environment of the other…”. Here we are only concerned with the relaxation of the system oscillators and do not require that any part the environment relaxes (In fact, in Appendix C we show that the bath modes never reach a steady state).

There is an important distinction in the meaning of equilibration. In the pure state quantum statistical mechanics paradigm equilibration is used more broadly to indicate relaxation to a steady state. For instance, depending on the context, the relaxation of the expectation values of certain operators to fixed values or of the reduced density matrix is considered equilibration. In our open system approach equilibration has a very specific meaning, see Eq. (162). In other words, there exists an environment and an interaction Hamiltonian such that the equilibrium state is obtained by tracing out the environment in the global thermal state.

Integrability. It has been discussed in GE16 that integrable quantum models indeed equilibrate to a suitable generalized Gibbs ensemble. Furthermore CEF11 examine the behavior of the one- and two-point correlation functions after a quench in various models, and it is found that the relaxation dynamics and equilibrium values can be well understood by means of a generalized Gibbs ensemble.

I.3 Key Results

1. Energy extensivity. In conventional thermodynamics, when the intra-system coupling is negligible, the internal energy is extensive in term of the number of the oscillators, like the case of the dilute gas. When this coupling is finite, we may instead understand the extensive property of the internal energy in terms of the normal modes of the coupled oscillators. We have shown that with this definition of extensivity the internal energy becomes extensive after the system reaches equilibrium, as implied by (38). It is interesting to note that the degrees of freedom of the oscillators used to describe the extensive property of the internal energy are neither the original degrees of freedom associated with each oscillator, nor the modes that decouple their equations of motion. Rather they are the degrees of freedom that diagonalize the oscillation frequency matrix $\boldsymbol{\Omega}_{p}^{2}$. In this regard, the extensive property of the internal energy in the final equilibrium state is the same as that of coupled oscillators in conventional thermodynamics, that is, in the vanishing system-bath coupling limit. This offers an explicit theoretical justification, from the open-system viewpoint, of conventional thermodynamics when applied to such a many-body system.

2. Heat capacity. When the system of $N$ coupled oscillators in a shared scalar field bath reaches equilibrium, its heat capacity is shown to be always non-negative for all nonzero bath temperatures, and it moves towards zero only if the bath temperature approaches zero. These properties are independent of the spatial arrangement of the oscillators, the inter-oscillator coupling and the system-bath interaction strength, as long as the collective non-Markovian motion of the system is stable.

3. The third law. Therefore from the viewpoint of behavior of the heat capacity at $T=0$ for this class of systems in an equilibrium state the third law is not violated. In this connection we also addressed the issue of entanglement and the third law pertaining to heat capacity. It was stated in Vedral the following: “ One may therefore say that in these systems the validity of the third law of thermodynamics relies on quantum entanglement …”. Our view is that The third law depends on the nondegeneracy of the ground state manifold and has nothing to do with entanglement directly. Indeed it has been shown HHPRD in the case of two spatially separated but coupled oscillators in a zero-temperature shared bath that the equilibrated state of this two oscillator system is not always entangled. For example, with sufficiently strong oscillator bath interaction, the reduced state of the two oscillators is separable. (See Fig. 3 in HHPRD .) Here, we have shown that the heat capacity of the coupled harmonic oscillator system goes to zero independent of the system-bath interaction strength. Thus it offers a counterexample to the above claim, that the validity of the third law of thermodynamics relies on quantum entanglement.

II Brownian motion of systems of oscillators strongly coupled to an environment

We now begin our detailed model study for considering the viability in the establishment of a thermodynamics theory of open quantum systems.

Consider a collection of coupled quantum harmonic oscillators in a shared finite-temperature $\beta^{-1}$ bath modeled by a massless scalar field in 1+3 Minkowski spacetime. The action of such a system is given by

[TABLE]

with $x=(t,\mathbf{x})$. Each oscillator is located at a fixed spatial coordinate $\mathbf{z}_{i}$, and has the same mass $m$ and bare natural frequency $\omega_{b}$. The “current” $j(x)$ in the oscillator-bath interaction term takes the form $\displaystyle j(x)=e\sum_{i}\chi_{i}(t)\,\delta^{(3)}(\mathbf{x}-\mathbf{z}_{i})$, with $e$ the coupling strength between the oscillator and the bath. The parameter $\sigma$ is the strength of direct coupling between two oscillators and assumed to be positive for concreteness 666It can take either sign which only affects the interpretation of the normal modes. In addition, the numerical values of $e$ and $\sigma$ are confined to ranges where instability in the dynamics is avoided. We will comment on this point later..

Here we suppose that the initial state of the combined system is a factorized state, given by

[TABLE]

where $H^{(\phi)}$ is the free Hamiltonian of the scalar field. While the field is initially prepared in a thermal state, the initial state $\rho_{i}^{(\chi)}$ of the system can be quite arbitrary. Thus in the beginning the system and the bath are not in equilibrium, nor correlated. We will let them interact and evolve in time. We will explore and make explicit the conditions when the system can and will reach equilibration 777The equilibration issue for classical coupled oscillator systems was studied before by, e.g., Agarwal Agarwal . This equilibrium state in general will have no resemblance to the thermal state of the combined system, nor of the reduced system. Thus the setup here is in strong contrast to the closed system globally-thermal state (CGTs) often adopted in the discussions of quantum thermodynamics. There, for the total Hamiltonian of the combined system $H=H^{(\chi)}+H_{\mathrm{int}}+H^{(\phi)}$ it is assumed that

[TABLE]

In the global thermal state, the system has already established correlation with the bath, and the interaction between the system and the bath is such that it maintains this correlation throughout. Since they are in thermal equilibrium, the combined system will remain in the global thermal state unless an external disturbance is introduced to bring the system out of equilibrium.

The evolution of the combined system is governed by the unitary evolution operator $U$

[TABLE]

where $\operatorname{T}$ denotes chronological time-ordering and $H$ is the Hamiltonian operator of the combined system that corresponds to the action (II). Given the initial state (2) of the total system, the density matrix of the reduced system of interest is then given by

[TABLE]

after we trace out the degrees of freedom of the bath. The reduced density matrix of the system enables us to calculate the quantum expectation values of the operators, say $O^{(\chi)}$, associated with the system by

[TABLE]

from which we may construct the quantum thermodynamics of the system in a nonequilibrium setting.

When the initial state (2) is Gaussian, Eq. (5) can be evaluated analytically and exactly for the combined system described by (II). Using a path integral representation of $U$ and $U^{\dagger}$, the reduced density matrix elements in (5) become

[TABLE]

where $S^{(\chi)}$ is the action of the system alone and $\chi_{\pm}$ denotes the system variable in the respective forward and backward time branches. This is where the ‘closed-time-path’ integral (CTP) / ‘in-in’ formalism CTP or its close kin, the Feynman-Vernon FeyVer influence functional (IF) $\mathcal{F}[\chi_{+},\chi_{-}]$ becomes particularly useful.

For a Gaussian bath, its influence on the system can be understood as caused by a classical noise by way of the Feynman-Vernon Gaussian identity: the imaginary part of the IF can be represented by a stochastic source term which inherits the quantum statistics of the bath. Using techniques from the CTP formalism, a revised imaginary part combined with the original real part of the influence functional, together with the action of the system, form a new effective action which is real, known as the stochastic effective action. Variation of this stochastic effective action produces a Langevin equation which describes the evolution of the reduced system. For details and working examples in this functional approach to open quantum system dynamics, please refer to Appendix A. In what follows, we will adopt the Langevin equation approach in the discussion of quantum thermodynamics at strong coupling.

II.1 Langevin Equation for the Reduced System

Following this well-established procedure the Langevin equation for the stochastic dynamics of the $i^{\mathrm{th}}$ oscillator (strongly) interacting with other oscillators and their shared bath with action (II) is given by

[TABLE]

In addition to the drag force and the quantum fluctuations of the bath found in a single oscillator system a new factor entering in the present coupled-oscillator shared-bath system is the induced interaction between the oscillators through their respective interaction with the scalar field bath. The field-environment mediated effect is non-Markovian in nature (see e.g., ASH06 ; LinHu09 ; HH15 ), often absent 888Unless the spatial information is retained in the system-bath interaction. in a shared bath modeled by a collections of oscillators 999The differences between a oscillator bath and a field bath will be discussed in Appendix E. (see, e.g., CHY ; PazRec ). This feature introduces additional complications and brings forth new physics in analyzing the stochastic dynamics of the quantum many-body system, as well as its quantum thermodynamics.

The statistics of the Gaussian noise field $\xi_{i}(t)=\xi(t,\mathbf{z}_{i})$ is determined completely by the first two moments

[TABLE]

All the higher even moments can be expressed by the second moment with the Wick expansion while all the odd moments vanish. The $\langle\cdots\rangle$ notation denotes either an ensemble average or expectation value, depending on whether the variable under consideration is stochastic or quantum. The two kernel functions $G_{R}^{(\phi)}(x-x^{\prime})$ and $G_{H}^{(\phi)}(x-x^{\prime})$ are most relevant for our present study: They are the retarded and the Hadamard functions of the scalar field $\phi$ in its thermal state, defined by

[TABLE]

with $\tau=t-t^{\prime}$ and $r=\lvert\mathbf{x}-\mathbf{x}^{\prime}\rvert$. Since they are time-translation invariant, their Fourier transforms with respect to the $\tau$ variable satisfy the well-known relation,

[TABLE]

where the Fourier transformation of the function $f(\tau)$ is defined by

[TABLE]

Introducing the matrix representation of the equation of motion (8),

[TABLE]

where $[\mathbf{G}(\tau)]_{ij}\equiv G(\tau,\mathbf{z}_{i}-\mathbf{z}_{j})$, we obtain a matrix equation

[TABLE]

The solution generically takes the form

[TABLE]

where $\{\boldsymbol{\Xi}(0),\dot{\boldsymbol{\Xi}}(0)\}$ are the initial conditions and $\mathbf{d}_{i}(t)$ are a special set of homogeneous solutions to (15). The actual form of $\mathbf{d}_{1}$ is not important but the Fourier transform of $\mathbf{d}_{2}(t)$ is

[TABLE]

Later it will be shown that for certain choices of parameters, the solution to (15) can exhibit instability and grows indefinitely when $t$ approaches infinity. In these case, the homogeneous solutions $\mathbf{d}_{i}(t)$ are not integrable,

[TABLE]

so their Fourier transforms do not exist in the usual sense. Thus when results are expressed in terms of $\tilde{\mathbf{d}}_{i}(\omega)$, it pays to be careful about their interpretations.

Finally we note from (16) that the oscillator is driven not only by the local noise at its very location, but is also affected by the quantum fluctuations of the bath at the locations of the other oscillators. This intriguing feature is essential in keeping the energy balance of the reduced system after it equilibrates. This will become clearer when we calculate the energy balance in Sec. III.3.

II.2 Covariance Matrix

From (16), if there exists an equilibrium state 101010The existence of the equilibrium state is related to the fact that the complex poles of $\tilde{\mathbf{d}}_{2}(\omega)$ lie on the upper half of the complex $\omega$ plane. See Sec. III.5. This is also the very basis on which we can discuss the fluctuation-dissipation relation of the reduced system and the energy balance among the dissipative, retarded and noise force terms. for the reduced system, then the moment

[TABLE]

at late time, after the reduced system completely relaxes, is well defined and given by

[TABLE]

where the superscripts $T$ and $\dagger$ denote the transposition and Hermitian conjugate of the matrix, respectively.

Since $\mathbf{d}_{2}$ is a symmetric matrix, we observe that

[TABLE]

with the help of the matrix identity

[TABLE]

for two nonsingular matrices $\mathbf{A}$, $\mathbf{B}$. Thus we can write $\boldsymbol{\sigma}_{\chi\chi}(\infty)$ in (20) as

[TABLE]

where the retarded Green’s function $\tilde{\mathbf{G}}_{R}^{(\chi)}(\omega)$ of the reduced system is in fact

[TABLE]

Similarly we introduce

[TABLE]

and at late times it becomes

[TABLE]

This integral in general is not well-defined due to the presence of ultraviolet (UV) divergence so regularization is needed.

II.3 Internal Energy

We define the internal energy of the system as its total mechanical energy. The total mechanical energy $E$ of the coupled oscillators is

[TABLE]

Here $\operatorname{Tr}$ is the matrix trace, and the matrix $\boldsymbol{\Omega}_{p}^{2}$ is defined in a way similar to $\boldsymbol{\Omega}_{b}^{2}$ except that the elements $\omega_{b}^{2}$ in $\boldsymbol{\Omega}_{b}^{2}$ are replaced $\omega_{p}^{2}$, where $\omega_{p}$ is the renormalized or physical frequency, which will be determined by the system preparation at the experimental energy scale. The difference between them is not necessarily large and depends on the choice of the cutoff frequency $\Lambda$, such that

[TABLE]

where the damping constant $\gamma$ is equal to $\gamma=e^{2}/(8\pi m)$. In the equilibrium state (note it is not necessarily the Gibbs state SFTH ), the total mechanical energy becomes

[TABLE]

The heat capacity $C$ is then given by

[TABLE]

The evaluation of $E(\infty)$ can be trickier than expected if regularization is not properly introduced.

At this point it may be desirable to get some physical feel of the dynamics and thermodynamics of the system. In Appendix B we treat a simpler system of one and two oscillators so that we can see the subtleties involved in the non-Markovian dynamics and thermodynamics of a strongly interacting open quantum system. Otherwise, we may proceed to the formal development for the $N$ oscillator system.

III Thermodynamics of Open Quantum Systems

As mentioned in the beginning, in this paper we use the model of an $N$ coupled-oscillator system interacting with a scalar field bath to address the energy and heat capacity issues and discuss the third law of thermodynamics.

Before we launch our studies of the $N$ coupled-oscillator model in full rigor, it would be useful to gain some feeling of the anticipated physical results for simpler cases. Hence, we summarize what we have learned from the one and two coupled harmonic oscillators examples below. Details for these two cases are placed in Appendix B.

For a system containing just one harmonic oscillator coupled to a thermal bath with finite coupling strength the heat capacity behaves qualitatively different at low temperatures from traditional thermodynamics, which assumes that the system-bath coupling is vanishingly small. A new scale associated with the coupling strength $\gamma$ appears. As shown in (136), the heat capacity approaches zero following a power law when the bath temperature is lowered to zero. In contrast, quantum statistical mechanics calculations assuming vanishing system-bath coupling predict in (141) that the heat capacity approaches zero exponentially as the bath temperature is lowered to zero. This is mostly transparently seen in Fig. 1 and 2.

With a mere increase of the number of system oscillators from one to two, the physics of the reduced system becomes more intricate because the two oscillators will have, on top of their direct coupling, also an indirect coupling mediated by the ambient scalar field bath, which introduces non-Markovian effects in the reduced system dynamics. As for quantum entanglement, in addition to the system-bath entanglement in the one-oscillator case, one needs to consider also entanglement between the constituent oscillators. Noteworthy on this issue is, as shown in Eq. (158) and Fig. 3-(c) of Appendix B, the behavior of the heat capacity for the two-oscillator system near absolute zero temperature does not depend on the presence or the absence of quantum entanglement between the two system oscillators. The heat capacity still approaches zero no matter what, and has a qualitatively similar behavior as the one-system-oscillator case.

III.1 System of $N$ Coupled Oscillators in a Common Bath

Now we consider a system that contains $N$ coupled harmonic oscillators in a shared thermal bath. Their spatial locations, specified by $\mathbf{z}_{i}$ with $i=1$, $\cdots$, $N$, is arbitrary, and their initial states can be far from equilibrium. From the previous discussions, we have learned that their motion is highly non-Markovian and intertwined, so it is not obvious whether systems that contain a large number of constituents always equilibrate. This would be the most important issue to address, namely, identify the conditions, or lack thereof, for an $N$ coupled harmonic oscillators in a shared thermal bath to reach in time an equilibrium (note different from thermal) state. We will show that indeed it exists. Then in this equilibrium state, we can discuss for this non-Markovian system the fluctuation-dissipation relation and the energy balance. We then advance towards the thermodynamics issues, beginning with a proof of the extensivity of the internal energy, the positivity of the heat capacity and finally, the behavior of the heat capacity as the temperature approaches to zero, pertaining to the issues of the third law.

The number $N$ of the system constituents can be arbitrary but cannot be infinite, because when it is comparable with the number of degrees of freedom of the bath, 1) it may lose the character of a system in contradistinction to its environment, as the basic definition of open systems calls for. 2) The system and environment should in this situation be considered as two equal subsystems interacting with each other which have a very different dynamics from open systems, e.g., recurrence; More seriously 3) the system may never equilibrate because any oscillator will be continually perturbed by the non-Markovian influences from its faraway counterparts all the time. This will make the motion of the system difficult to settle down.

For a finite $N$, following our earlier analysis outlined in the two-oscillator case, we note there are exceptional cases that equilibration may not be always possible. For example, we exclude those arrangements where some of the oscillators are placed remotely from all others because such a setup can render the relaxation time unusually long. From these considerations we assume the number $N$ is much smaller than the number of degrees of freedom of the bath, and that the oscillators are all localized within a finite region. As a reminder, this still does not exclude the possibility that when the non-Markovian effects are sufficiently strong, albeit not enough to induce instability, the system might still have an extraordinarily long relaxation time so as to behave almost like an undamped one.

Since analytical results for the $N$-oscillator system are unavailable, we will provide a qualitative but general analysis based on the mathematical properties of positive matrices. We start with two simpler topics by first examining the fluctuation-dissipation relation of the reduced system in the final equilibrium state, and then the energy balance between the reduced system and the bath. They provide the basis for extensivity of internal energy and positivity of heat capacity. We will save the discussion on the existence of this equilibrium state for the end.

III.2 Fluctuation-Dissipation Relation, Stationarity

We direct our attention now to the correlation function of $\chi_{i}(t)$ and derive the corresponding fluctuation-dissipation relation when the reduced system reaches equilibrium. From (16), we find the correlation function, namely, the Hadamard function of $\boldsymbol{\Xi}(t)$ given by

[TABLE]

Again $[\mathbf{G}_{H}^{(\phi)}(s-s^{\prime})]_{ij}=G_{H}^{(\phi)}(s-s^{\prime},\mathbf{z}_{i}-\mathbf{z}_{j})$. It is not invariant in time translation so the intermediate state is not an equilibrium state. If we choose the parameters of the configuration in such a way that no runaway solution is allowed, then $\mathbf{d}_{i}(t)$ exponentially decays with time. Thus in (31), those terms that are not inside integrals will be exponentially small at late times. The double integrals in (31) can be written as

[TABLE]

among which we have ignored terms that are exponentially small at late times and have used the approximation that when $t$ is sufficiently large,

[TABLE]

with $\alpha$ being some positive number to describe the generic decaying behavior of $\mathbf{d}_{2}$ with time. Thus we see the nonstationary components in $\mathbf{G}_{H}^{(\chi)}(t,t^{\prime})$ becomes negligibly small as $t$, $t^{\prime}\to\infty$. We can then focus on the stationary component,

[TABLE]

where we have invoked the fluctuations-dissipation relation of the free (stand-alone) scalar field

[TABLE]

In (34), we notice that the integrand in fact is $\tilde{\mathbf{G}}_{H}^{(\chi)}(\kappa)$ by the definition of the Fourier integral, and thus we arrive at

[TABLE]

when the reduced system reaches equilibrium. Thus, from the derivation we see that the correlation function of the reduced system is not stationary in time during the nonequilibrium evolution, but dissipation causes the nonstationary component of the correlation to decay with time such that when the dynamics of the reduced system is relaxed, the correlation becomes stationary. This reflects the presence of a final equilibrium state.

Stationarity enables us to express the fluctuation-dissipation relation of the reduced system in the frequency domain, similar to that of the bath. However, even though they appear deceptively similar in structure, they are utterly different in physical contents. Essentially (35) is based on the initial thermal state of the bath, while (36) is established only because there exists a final equilibrium state, which by no means is necessarily a Gibbs thermal state; however it still inherits the information of the initial thermal state of the bath. This is related to the fact the late-time statistics of the reduced system is governed by the bath. A similar behavior is also observed for the case when a charged oscillator interacts with a quantized electromagnetic field, initially prepared in a squeezed vacuum WHL12 . The only difference is that the proportionality constant in the fluctuation-dissipation relation like (36) takes a different form and depends on the squeeze parameters of the bath’s initial squeezed vacuum state.

III.3 Energy Balance in the Equilibrium State

Now we turn to the energy balance of the reduced system described by (8) in the equilibrium state,

[TABLE]

where complexity arises from the frequency renormalization and the nonlocal causal influence among oscillators. In the single oscillator case, when equilibrium is reached, the net energy flow between the oscillator and the bath stops. The energy flowing in from the noise force of the bath is counterbalanced by the energy flowing out of the oscillator due to the frictional force, as captured by the fluctuation-dissipation relation. In the multi-oscillator case, it is then interesting to ask whether only the same two factors are needed to balance the energy flow in the course of equilibration, or other mechanisms are also involved? If so, what are their roles in the fluctuation-dissipation relation?

We will show that

[TABLE]

that is, the energy transfer mediated by the shared bath ceases after the motion of the reduced system reaches equilibrium.

We first rewrite the third term in (III.3) as,

[TABLE]

where $\Gamma^{(\phi)}(t)$ vanishes for a scalar-field bath at late times, and we have introduced a new kernel function $\Gamma^{(\phi)}(\tau)$,

[TABLE]

In what follows we will calculate the power delivered to the $i^{th}$ oscillator in the equilibrium state.

Each of the terms in (39) gives a contribution with a distinct physical interpretation. The first term on the righthand side of (39) will be absorbed into the bare frequency $\omega_{b}$ to form the physical frequency $\omega_{p}$

[TABLE]

The second term on the righthand side of (39) thus represents the dissipative force whose mean power delivered to the $i^{\mathrm{th}}$ oscillator is

[TABLE]

The mean power exerted by the noise force on the $i^{th}$ oscillator is

[TABLE]

Finally the net power delivered by the other oscillators to the $i^{th}$ oscillator via the nonlocal causal influence transmitted by the field is given by

[TABLE]

These three contributions look very distinct in nature, but we will show that at late times after the system of oscillators relaxes, their sum vanishes. Let us rewrite (42)–(44) in the limit $t\to\infty$,

• •

$P_{\gamma}^{(i)}(\infty)$ : It is given by

[TABLE]

where we have used several facts

1. (a)

In general $\mathbf{G}_{H}^{(\chi)}(t,s)=\langle\{\chi_{i}(t),\,\chi_{i}(s)\}\rangle/2$ is not invariant with time translation unless $t$, $s$ are sufficiently large. That is, the non-stationary components will decay with time, so when $t$, $s\to\infty$, we can write $\mathbf{G}_{H}^{(\chi)}(t,s)$ into $\mathbf{G}_{H}^{(\chi)}(t-s)$. 2. (b)

The real part of the Fourier transform of a retarded Green’s function $\tilde{\mathbf{G}}_{R}^{(\phi)}(\kappa)$ is an even function in $\kappa\in\mathbb{R}$ , but the imaginary part is an odd function. 3. (c)

\operatorname{Im}\bigl{[}\tilde{G}_{R}^{(\phi)*}(\kappa,\mathbf{0})\bigr{]}_{ii}=-\operatorname{Im}\bigl{[}\tilde{G}_{R}^{(\phi)}(\kappa,\mathbf{0})\bigr{]}_{ii}.

• •

$P_{\xi}^{(i)}(\infty)$ : It is given by

[TABLE]

where we have made use of the fluctuation-dissipation relation (36) for the reduced system.

• •

$P_{c}^{(i)}(\infty)$ : It is given by

[TABLE]

We observe that unlike the one-oscillator case,

[TABLE]

so in the multi-oscillator case, the energy balance is more delicate. On the other hand, the contribution $P_{\gamma}^{(i)}(\infty)$ can be combined with $P_{c}^{(i)}(\infty)$ to form

[TABLE]

which turns out to be the negative of $P_{\xi}^{(i)}(\infty)$. We thus see in fact we should have

[TABLE]

if both of the fluctuation-dissipation relations

[TABLE]

hold.

Eq. (50) immediately implies (38). Here we see additional mechanisms are at play in the energy transfer between coupled oscillators. The motion of any oscillator is, apart from direct coupling, causally affected by all the other oscillators via the shared bath. These coherent and correlated contributions from the other oscillators, depending on their individual evolution history, do not necessarily induce a drag nor a push force on that very oscillator. The net effects of the retarded influence are thus highly complicated, hinging on the distance between any two oscillators and their states of motion. It is not obvious how they participate in balancing the energy flow between each system oscillator and the bath. However, we have mentioned earlier that each oscillator, in addition to experiencing the disturbance from the noise of the bath locally, is also affected by the bath fluctuations at the locations of the other oscillators. We can see that the correlations of the bath fluctuations will be passed on to the oscillators such that their motions are also correlated. These correlated noises can be the counterparts of the causal influences, both of which are the off-diagonal elements in the fluctuation-dissipation relation of the bath (51), in the roles of either the fluctuation-dissipation relation or the energy balance of the system. Moreover once we observe that the retarded influence is in fact related to the Liénard-Wiechert-type radiation of the scalar field as a consequence of the oscillators’ motion, it is clear that the damping force is of the same physical origin as the non-Markovian causal influences. Thus grouping $P_{\gamma}^{(i)}(\infty)$ and $P_{c}^{(i)}(\infty)$ together in (49) becomes natural, and from that (50) follows.

Eq. (38) says that even with the presence of non-Markovian influences in the motion of the reduced system, the interaction of the system with the bath is such that when the system settles down in its equilibrium state, its total mechanical energy becomes constant in time. Hereafter the reduced system acts as a collection of coupled undamped oscillators, oscillating at the physical frequency $\omega_{p}$, and can be completely described by the final equilibrium density matrix. That is, the reduced system becomes self-contained and free from any further intervention from the bath. This motivates us to assign the total mechanical energy as the internal energy of the system.

Next we will discuss the extensive property of the internal energy of a system of $N$ coupled oscillators in a shared bath.

III.4 Extensivity of Internal Energy

Before proceeding to the coupled system in a nonequilibrium configuration, we first delineate the extensivity of the internal energy in the simpler equilibrium case.

Formally, equilibrium thermodynamics is realized in the limit $\gamma\to 0$, whereby the matrix $\tilde{\mathbf{G}}_{R}^{(\chi)}(\omega)$ reduces to

[TABLE]

where in order to preserve the retarded property of $\mathbf{G}_{R}^{(\chi)}(\tau)$, we have let $\omega\to\omega+i\,\epsilon$ with $\epsilon>0$. Since the matrix $\boldsymbol{\Omega}^{2}_{p}$ is real and symmetrical, we may find a real orthogonal constant matrix $\mathbf{U}$, independent of $\omega_{p}$ and $\sigma$, to diagonalize it, that is,

[TABLE]

The matrix $\mathbf{W}_{p}$ is real and diagonal and we have assumed that its diagonal elements remain positive definite, with the appropriate choice of $\omega_{p}$ and $\sigma$ to avoid instability in motion. Thus we write (52) as

[TABLE]

with the diagonalized $\boldsymbol{\mathfrak{G}}_{R,\,0}^{(\chi)}(\omega)$ matrix given by

[TABLE]

Now since the symmetric matrices $\boldsymbol{\sigma}_{\chi\chi}(\infty)$ and $\boldsymbol{\sigma}_{\upsilon\upsilon}(\infty)$ are related to $\tilde{\mathbf{G}}_{R}^{(\chi)}(\omega)$, according to (23) and (26), we can write them into the diagonal forms as well with the help of $\mathbf{U}$,

[TABLE]

in which the diagonal matrices $\mathbf{n}$, $\mathbf{m}$ are

[TABLE]

So far what we have done is equivalent to expressing the results in terms of the normal modes of the coupled oscillators when their interaction with the shared bath is almost nonexistent. The matrices $\mathbf{n}$, $\mathbf{m}$ are nothing but the position and velocity uncertainties of the normal-mode coordinates. This decomposition implies that the mean mechanical energy in the equilibrium thermal state, due to the presence of the trace, is invariant under the orthogonal transformation acted by $\mathbf{U}$,

[TABLE]

The advantage of the form (59) is that since every matrix in it is diagonal, (59) can be literally and formally written as

[TABLE]

where

[TABLE]

is essentially the mechanical energy associated with each normal mode. That is, the total mechanical energy is the sum of the mechanical energy of each normal mode. Thus when the interaction between the coupled oscillators and the shared bath is negligible, the mechanical or internal energy is extensive, at least with respect to the normal modes. This is the limiting condition underlying conventional thermodynamics.

When the oscillator-bath interaction is not negligible, the Green’s function matrix $\tilde{\mathbf{G}}_{R}^{(\chi)}(\omega)$ of the oscillators contains the contribution from the retarded Green’s function matrix $\tilde{\mathbf{G}}_{R}^{(\phi)}(\omega)$ of the free scalar field,

[TABLE]

Since the values of the elements of the matrix $\tilde{\mathbf{G}}_{R}^{(\phi)}(\omega)$ depend on the locations of the coupled oscillators,

[TABLE]

the orthogonal matrix $\mathbf{U}$ that can diagonalize $\boldsymbol{\Omega}_{b}^{2}$ in general cannot diagonalize $\tilde{\mathbf{G}}_{R}^{(\phi)}(\omega)$, because the latter two matrices do not commute in general unless the locations of the oscillators are especially arranged. That is, in general the matrices

[TABLE]

are not diagonal, neither are $\mathbf{n}$ and $\mathbf{m}$ in this case. Even though (59) always holds, and it will give an impression that the total mechanical energy can still be expressed as a sum like (60), here $E_{i}$ does not enjoy the special significance of the mechanical energy of each normal mode for the full equation of motion (15)

[TABLE]

This may be most easily understood if we apply the transformation $\mathbf{U}$ to the coupled equations of motion (65), and it becomes

[TABLE]

where $\mathbf{X}=\mathbf{U}^{T}\cdot\boldsymbol{\Xi}$ represents the coordinates of the normal modes of the coupled oscillators in the absence of the bath as is discussed earlier, but is not the normal modes of the coupled oscillators in the presence of the shared bath. The off-diagonal elements $\boldsymbol{\mathfrak{G}}_{R}^{(\phi)}(t-s)$ will link up any given element in $\mathbf{X}$ with all other elements.

Therefore, from (38) we arrive at some interesting conclusions. When the coupled oscillators interact with a shared bath, after the coupled system reaches equilibrium, the internal energy of the system oscillators becomes extensive; however this extensivity is expressed by neither its original degrees of freedom nor the decoupled degrees of freedom. Instead, the extensive property of the system’s internal energy is only manifested by a specified set of modes obtained from the orthogonal transformation that diagonalizes $\boldsymbol{\Omega}^{2}_{p}$, as can be seen from Eq. (38). Moreover, before the motion of the reduced system equilibrates, this special extensiveness property does not hold, as is implied by (38). Thus we are not able to discuss the extensive property of the system’s internal energy during the nonequilibrium evolution of the reduced system, until the final equilibrium state of the reduced system is attained.

In particular, this seemingly mundane conclusion, together with (38), justifies or explicitly demonstrates, in the weak oscillator-bath coupling, why conventional thermodynamics (at least for the system that constitutes coupled oscillators) works, why we need only the density matrix of the system to describe the behaviors of the system, and why we need not be concerned with renormalization, relaxation, damping, bath noise.

III.5 Positivity of Heat Capacity and Existence of the Equilibrium State

Now we would like to discuss the positivity of the heat capacity for a system of $N$ coupled oscillators in a shared bath in the context of nonequilibrium thermodynamics. The positivity of heat capacity, the decaying behavior and the retarded nature of $\mathbf{d}_{2}(t)$ all hinge on the existence of the equilibrium state. Thus in this section, we will also address the conditions that a nonequilibrium system settles into an equilibrium state at late times.

Given the internal energy (29) when the system reaches equilibrium, we proceed to examine the positivity property of the heat capacity $C(\infty)$, given by

[TABLE]

from (29) and (30). The damping factor $e^{-\frac{\lvert\kappa\rvert}{\Lambda}}$, with $\Lambda\to+\infty$, is necessary to regularize the integral.

If we consider only the case that there is no runaway solution in the motion of the system, such as with the inverted oscillator, then this requires that the matrix $\boldsymbol{\Omega}_{p}^{2}$ should be at least positive definite 111111Positive semi-definiteness can be too weak because the non-Markovian contributions can easily induce instability in the strong system-bath coupling regime or in the limit of extremely short separations among the oscillators.. This allows us later to define a matrix that would be the square root of $\kappa^{2}\mathbf{I}+\boldsymbol{\Omega}_{p}^{2}$. The imaginary part of $\tilde{\mathbf{d}}_{2}(\kappa)$ can be written as

[TABLE]

Using (17) we know that $\tilde{\mathbf{d}}_{2}(\kappa)$ takes the form

[TABLE]

with the help of (22). The matrix $\tilde{\not{\mathbf{G}}}_{R}^{(\phi)}(\kappa)$ is $\tilde{\mathbf{G}}_{R}^{(\phi)}(\kappa)$ with its diagonal elements removed,

[TABLE]

The diagonal elements of $\tilde{\mathbf{G}}_{R}^{(\phi)}(\kappa)$ account for the usual damping term and the frequency renormalization. Now we introduce the matrix $\tilde{\mathbf{D}}_{2}(\kappa)$ by

[TABLE]

such that

[TABLE]

where we have used the cyclic property of the matrix trace. The product of the last pair of matrices is positive, so we would like to examine whether the matrix sandwiched by $\tilde{\mathbf{d}}_{2}^{\vphantom{\dagger}}(\kappa)$ and $\tilde{\mathbf{d}}_{2}^{\dagger}(\kappa)$ in (68) is positive as well. In general, it takes the form

[TABLE]

where $\ell_{ij}=\ell_{ji}=\lvert\mathbf{z}_{i}-\mathbf{z}_{j}\rvert$. Here we assume that the choices of the parameters $\omega_{p}$, $\sigma$ and $\ell_{ij}$ are such that the matrix $\boldsymbol{\Gamma}$ is strictly diagonally dominant, that is, its elements satisfying

[TABLE]

At the first sight, this assumption looks pretentious; however we observe that the strictly diagonally dominant matrix has a nice property of being positive definite HJ85 . That is, its eigenvalues are all positive. Thus this assumption, together with positive definiteness of $\boldsymbol{\Omega}_{p}^{2}$, implies that the integrand in (67) is always positive. Since the integral is well defined, we conclude the value of the heat capacity $C(\infty)$ remains positive for all temperatures $\beta^{-1}$ with one exception that $\beta\to\infty$. In that limit the factor

[TABLE]

goes to zero 121212We can generalize the more sophisticated arguments in the context of the two-oscillator system to the current case., so it indicates that the heat capacity for the system of coupled oscillators will be zero at zero temperature; otherwise it is always positive.

Physically the assumption (73) amounts to the existence of the effective damping constants for all modes of motion, and thus the motion of the system, described by $\tilde{\mathbf{d}}_{2}$ reduces to that of a collection of coupled damped oscillators. This can be read off from the denominator of $\tilde{\mathbf{d}}_{2}(\kappa)$. Suppose the real matrix $\boldsymbol{\Gamma}$ can be diagonalized by the orthogonal matrix $\mathbf{V}$

[TABLE]

The matrix $\mathbf{V}$ in general cannot diagonalize $\boldsymbol{\Omega}_{p}^{2}$ unless $\boldsymbol{\Omega}_{p}^{2}$ commutes with $\boldsymbol{\Gamma}$, but it will transform $\boldsymbol{\Omega}_{p}^{2}$ to another symmetric, positive matrix, which we denote by $\boldsymbol{\mathcal{W}}^{2}_{p}$. The diagonal elements of $\boldsymbol{\mathcal{W}}^{2}_{p}$ describe the same physical frequencies $\omega_{p}$ of the transformed modes and the off-diagonal ones account for the coupling among them. The mode-mode couplings are usually different among pairs of modes. Explicitly the denominator of $\tilde{\mathbf{d}}_{2}(\kappa)$ is transformed to

[TABLE]

and the zeros of its determinant

[TABLE]

identify the eigen-modes of the motion of the system. The signs of the imaginary part of the solutions to (76) provide information about the stability of the motion. If there exists a solution whose imaginary part is positive, then instability of the collective motion 131313Unless the unstable mode is not excited, and that is highly unlikely for a generic initial state. will occur. Eq. (III.5) in fact corresponds to a simultaneous set of equations of motion that describes a system of coupled damped oscillators,

[TABLE]

Thus the stability condition associated with (77) is equivalent to whether the characteristic polynomial (76), when $\kappa=-i\,s$, is a (strict) Hurwitz polynomial AH95 , whose zeros are all located on the left half of the complex $s$ plane. In other words, the motion described by (77) is stable if the characteristic polynomial associated with the Laplace transformation of the lefthand side of (77)

[TABLE]

is Hurwitz. In general, a sufficient and necessary condition is provided by the Routh-Hurwitz stability criterion PD85 , which states that all principal minors of the Hurwitz matrix associated with $p(s)$ are positive.

This criterion becomes computationally cumbersome as $n$ grows, and it is very hard to establish an apparent connection between this criterion and the physically meaningful matrices $\boldsymbol{\Gamma}^{\prime}$ and $\boldsymbol{\mathcal{W}}^{2}_{p}$. For this reason we turn to finding arguments to directly relate the properties of the matrices $\boldsymbol{\Gamma}^{\prime}$ and $\boldsymbol{\mathcal{W}}^{2}_{p}$ with the stability condition of motion described by (77). These arguments, although mathematically less rigorous, are physically more transparent. The idea is that solving the polynomial $p(s)=0$ is equivalent to finding the eigenvalue $s$ of the system VE11

[TABLE]

with the normalized column eigenvector $\mathbf{x}$, with $\mathbf{x}^{T}\cdot\mathbf{x}=1$. We multiply (79) from the left with $\mathbf{x}^{T}$, transforming the matrix expression (79) to an ordinary quadratic equation of $s$,

[TABLE]

so that

[TABLE]

Since we have required that $\boldsymbol{\Gamma}^{\prime}$ and $\boldsymbol{\mathcal{W}}^{2}_{p}$ are (strictly) positive definite, the variables $b$ and $c$ are also positive by construction. This implies that

[TABLE]

Thus the positive definiteness of $\boldsymbol{\Gamma}^{\prime}$ and $\boldsymbol{\mathcal{W}}^{2}_{p}$ is sufficient to ensure the stability of the motion (77) which in turn signals the existence of an equilibrium state. In addition, the expressions in (80) resemble those we have seen for the case of one oscillator interacting with a bath, where $\sqrt{c-b^{2}}\propto\sqrt{\smash[b]{\omega_{p}^{2}-\gamma^{2}_{\vphantom{p}}}}$ is related to the resonance frequency.

In summary, the requirement that $\boldsymbol{\Gamma}^{\prime}$ and $\boldsymbol{\mathcal{W}}^{2}_{p}$ are positive matrices implies that $\mathbf{d}_{2}(\tau)$ is indeed a retarded Green’s function and $\tilde{\mathbf{d}}_{2}^{\vphantom{\dagger}}(\kappa)$ does not have any pole along the real axis of $\kappa$ and on the upper half of the complex $\kappa$ plane. Therefore the integrand in (67) is positive and bounded, so the heat capacity (67) is positive and approaches zero as $\beta\to\infty$. In addition it ensures the existence of the equilibrium state, which is needed a) for the reduced system to have a meaningful fluctuation-dissipation relation, b) to show the energy balance between the system and the bath, c) to ensure the extensive nature of the internal energy of the system and finally d) for the associated heat capacity to be positive definite in our framework of open system nonequilibrium dynamics approach to quantum thermodynamics.

IV Summary and Discussions

IV.1 Summary of major results

As a preamble we bring up the rather special conditions whereupon the foundation of thermodynamics is laid, from an open system perspective: A small open system interacting with a vast environment (whose thermal properties can be captured by a few physical parameters, its temperature, chemical potential), it is in the limit of vanishing coupling between them, only when the system can equilibrate and thermalize at late times, that thermodynamics makes sense. These considerations can be extended to nonequilibrium conditions but not for far from equilibrium, fully arbitrary time evolutions. We mentioned the important differences in the setups for treating quantum thermodynamics (QTD), namely, between level 1 Assuming the closed system (comprising the system and its environment) remains in a global thermal state (which we call CGTs) and level 2 an open system approaching equilibrium at late times (we call it ONEq), which is the preferred approach we adopt for the discussion of level 0. The centroid of this paper is a detailed model study, that of a system of $N$ coupled, spatially separated quantum oscillators interacting with a common scalar quantum field bath at finite temperature, where the existence of exact solutions can provide unambiguous quantification of physical variables, thermodynamic relations and help to clarify many basic issues in QTD. The set of issues we addressed include the following.

IV.1.1 Gateway to thermodynamics: The existence of an equilibrium state

1. Equilibrium state at late times: Let the system initially be prepared in a state that is not in thermal equilibrium with the shared bath, it has been known that if the coupling between the system and the bath is vanishingly weak, the reduced system will equilibrate at late times. This is the pre-condition for talking about its thermodynamic behavior. The new challenge is whether the system will equilibrate for strong coupling. This point has been emphasized in e.g. SFTH who used the quantum Brownian motion model where the system consists of $N$ quantum harmonic oscillators and the environment is an infinite-oscillator bath.

2. Equilibration, not thermalization: The strong coupling regime poses new challenges: Allowing the coupling between the system oscillators and the interaction between the system and the bath to be strong, and assuming that the dynamics of the system remains stable, the first and foremost statement is that due to non-weak system-bath interaction, this final state (of the system) is not described by a density matrix of the Gibbs form with respect to the system Hamiltonian. Therefore one should refrain from using the word thermalization to describe the end result, and note that conventional thermodynamics need not apply. The tough question is, when will TD remain a viable theory for this equilibrated strongly coupled system.

3. Environment-induced non-Markovian inter-oscillator interaction: The newer challenge which we need to take on here is to show equilibration for a system of strongly coupled $N$ quantum oscillators at finite spatial separation and strongly interacting with an environment composed of a quantum scalar field. The case of $N$ oscillators in the same spatial location is easier to prove because one needs not worry about the field-induced non-Markovian effects. However, beware of the pathology of even two oscillators stacked up at the same spatial location, as described in Appedix. E. The added complication is due to the non-Markovian nature of the induced interaction discovered in ASH06 ; LinHu09 ; HHPRD amongst the system oscillators (or qubits) mediated by the field environment. This issue has not been dealt with in this context before, as far as we know.

4. The existence of an equilibrated state for the case of two coupled oscillators has been demonstrated. The conditions for $N$ oscillators are discussed in Sec. III.5. We argue that certain positive-definiteness requirements must be satisfied to the effect that the effective damping constants of the oscillators stay positive and the effective oscillator frequencies remain real.

5. With the assurance of an equilibrated state, many nice properties follow. Specifically, the extensivity of the internal energy and the positivity of the heat capacity. The absence of such a state for open quantum systems severs the linkage to thermodynamics. QTD in the form described here does not exist for these systems.

IV.1.2 Internal energy, heat capacity and the third law

1. The internal energy for certain strongly bounded systems may not be straighforward to define (e.g., the presence of self energy as when gravity is involved) but fortunately not so in the model we studied: it is the sum of the kinetic energy of each oscillator, the harmonic potential energy, and their coupling energy. Heat capacity is the derivative of the internal energy with respect to the bath temperature.

2. We examine the third law from the behavior of the heat capacity at low and zero temperatures. We are concerned with a) low temperature behavior, b) the positivity and c) the extensivity of heat capacity.

3. The internal energy and the heat capacity for a system consisting of only one harmonic oscillator have been derived before in the CGTs setup HanIng06 . They are derived here in an open quantum system ONEq setup, which in the epoch after equilibration, can be compared, in the weak oscillator-bath coupling limit, with the quantities derived in the CGTs and in the conventional thermodynamics. They all agree with each other.

4. Complexity arises when the system has more than one constituent. The bath-induced non-Markovian effects cannot be properly described in conventional thermodynamics. This also brings in question the validity of energy extensivity because the system constituents are not only directly coupled (which is easy to deal with by normal mode separation) but are indirectly coupled or intertwined in a non-Markovian way by the induced interaction through their common environment.

5. With the proven existence of an equilibrated state for $N$ spatially separated but mutually coupled system oscillators we have shown that even a strongly coupled system can still have asymptotic extensivity of the internal energy, and the heat capacity remains positive as long as the motion is stable.

6. To compare results calculated in the three different levels, CGTs, ONEq and conventional thermodynamics, we need to beware of their respective regimes of validity and identify their common denominators, e.g., the common physical quantities and the states they are in. We work here with an open system ONEq set up, namely, we allow the system to evolve from a nonequilibrium initial state to a final equilibrium state. If the system thermalizes then the results obtained in the ONEq setup can be compared to results in conventional thermodynamics, as we did. Since the reduced density matrix of the ONEq after equilibration and that of the system in the CGTs setups are the same SFTH , the results from these two setups can be compared using this quantity.

IV.2 Heat capacity and third law

For the same system we have studied here there are claims of negative heat capacity in variance to our findings. For example, the author of HasegawaJMP working with a global thermal state CGTs setup claims that the heat capacity of the system of multiple quantum harmonic oscillators can be negative at low temperatures for 1) stronger system-bath interaction, or 2) smaller number of system constituents. However, we have demonstrated under rather general conditions that after the system reaches equilibration, the heat capacity of the system is always positive and approaches zero, for the full range of system-bath interaction strength. Thus the third law, from the aspect of low-temperature behavior of the heat capacity of the system, is not violated. Our findings of the extensivity of the energy and the positivity of heat capacity for $N$ coupled oscillator system should hold in all three levels of inquiry.

IV.3 On entanglement witnesses and heat capacity

In Ref. Vedral Wieśniak et al have “shown that the low-temperature behavior of the specific heat can reveal the presence of entanglement in bulk bodies in the thermodynamical equilibrium”. They drew this conclusion by showing that the heat capacity is an entanglement witness for some spin models. This involves finding a lower bound on heat capacity that can be achieved by separable states of the system. The third law requires the heat capacity to approach zero asymptotically with the temperature. This means that at low enough temperatures the behavior of heat capacity is not compatible with separable states, and is thus an indicator of entanglement.

Here we find that specific heat is not a reliable indicator of entanglement for our model.141414It is worth mentioning that the models studied in Ref. Vedral differ from the ones studied in this paper in two important aspects. First, we are studying harmonic oscillators with infinite dimensional Hilbert spaces as opposed to the finite dimensional spin systems. Second, we are taking an open system approach, allowing the system-bath coupling to be finite, and we are dealing with equilibrium states as opposed to thermal states. We see no requirement HHPRD that the zero temperature state is always entangled. For example, when the coupling between the constituents in the system is sufficiently strong, the system tends to relax to an entangled state at zero temperature. However, for sufficiently strong system-bath coupling, the system can relax to a separable state even at zero temperature. In this equilibrium state, the system oscillators are disentangled among themselves, but can get entangled with the bath oscillators. This can be understood as a consequence of the monogamy of entanglement.

At the same time, we have shown that the third law is valid for our model in all parameter regimes. There is no connection between the third law and entanglement. Moreover, even in the models or parameter regimes in which this connection exist, we do not interpret this observation as the third law relying on quantum entanglement. The third law stands. The claim that the third law implies the existence of entanglement could well be affected by the use of entanglement witness as a criterion.

IV.4 Relation with global thermal state formulation, sequel on heat, energy and entropy

IV.4.1 Relation to global thermal state formulation: Seifert’s systematics of energy heat and entropy for quantum systems

Much work on QTD has been done under the closed system in a global thermal state (CGTs) set up. It would be useful to find a link between it and our open system nonequilibrium (ONEq) approach. We have carried out a first step towards this goal. We focused on Seifert’s rendition of energy and entropy for classical thermodynamics Seifert16 in a closed system global thermal state (CGTs) set up. In a companion paper GTOS we have generalized his results for the thermodynamics of classical systems to quantum systems. This may enable us to use his systematics for discussing the First and Second Laws for quantum systems.

The diversity of how thermodynamic functions are defined is both a resource of adaptivity and at times a source of confusion. To show how the thermodynamic functions are used and how they enter into the TD relations, in this same companion paper we studied the approach of Gelin and Thoss GelTho who also work in the CGTs setup but adopt a different set of thermodynamics functions from Seifert’s. Consulting Seifert’s systematics and allowing for varying thermodynamic functions we hope to construct more links between our open systems ONEq formulation and the prevailing CGTs formulations of quantum thermodynamics.

IV.4.2 Sequel: On heat, entropy, entanglement and the second law

The issues of heat, entropy and entanglement in strongly coupled open quantum systems will be the center of attention in our second paper. Notice the subtle yet important difference between energy and heat. In this paper we have focused on the internal energy and heat capacity of the system, but not heat, which is the energy transfer between the system and the environment. This is because heat transfer as energy change contains ambiguities in an open system context. For example Esposito et al EspositoPRB15 found that any heat definition expressed as an energy change in the reservoir energy plus any fraction of the system-reservoir interaction is not an exact differential when evaluated along reversible isothermal transformations, except when that fraction is zero. Even in that latter case the reversible heat divided by temperature, namely entropy, does not satisfy the third law of thermodynamics and diverges in the low temperature limit.

We also have reservations in some claims of violation of the Second Law for quantum systems 2LawSaga . For example if one uses the Clausius inequality representation for the Second Law, we know it is only valid for classical systems at high temperatures. One needs to scrutinize the different definitions of entropy for strongly coupled quantum systems to make sure they are physically sound in the quantum regimes (such as at low temperatures), including possible non-Markovian behaviors, before adopting them to address foundational issues.

The relevant issues ranging from quantum correlations, entanglement, information to entropy production and heat can be sampled in these references of the last 15 years EisPle ; HB05 ; HBJSP08 ; HiltLutz ; ELvdB ; DeffLutz ; Kim12 ; EspositoJSM ; AurEic15 ; AnkPek14 ; EspositoPRB15 ; CarrWeiss16 , which span the scope of our sequel studies using the paradigm established in this paper.

Acknowledgement

B.L.H enjoyed the hospitality of the Physics Theory Center of Fudan University when this work was conceived while J.T.H enjoyed the hospitality of the Physics Theory Center at the University of Maryland when it was set in its final form. B.L.H thanks E. Bianchi for discussions related to his recent work on eigenvalue thermalization. C.H.C acknowledges the support from Ministry of Science and Technology, Taiwan under Grant No. MOST 105-2112-M-006-011 and National Center for Theoretical Sciences (NCTS). Y.S. acknowledges that this work was performed under the auspices of the U.S. DOE Contract No. DE-AC52-06NA25396 through the LDRD program at LANL.

Appendix A Influence Functional Formalism for Quantum Thermodynamics

For the study of the thermodynamical properties of a system one mainly focuses on the dynamics and thermal properties of the system under the influence of a thermal bath it interacts with, not particularly about the bath itself. This is the arena where the open system conceptual framework is most suited. When the system of interest interacts with the bath strongly, one needs to take into account not only the influence of the bath on the system but also the backreaction of the system on the bath. The influence functional formalism is particularly adept for this description because it respects the self-consistency of the system-bath evolutionary dynamics (an example being the fluctuation-dissipation relation) which has increased significance for treating systems with non-vanishing coupling with its environment.

In this appendix we give a self-contained description of the influence functional formalism used in this series of papers for obtaining the thermodynamic properties of physical quantities of a system interacting strongly with a thermal bath. Over and beyond the standard classic sources of Feynman-Vernon FeyVer and Caldeira-Leggett CalLeg , we also provide additional materials in the use of the coarse-grained effective action CGEA and stochastic effective action JH1 developed in the 1990s and 2000s. And, from the stochastic effective action, as a further development, we follow a recent work HH15 to expound the different advantages of using the Langevin equation route which is more intuitive versus the more formal route via the reduced density operator, which can account for the full quantum dynamics of the reduced system and enforce the operator ordering. This system of tools was used to describe the thermodynamics of quantum many-body systems in a nonequilibrium steady state. The reader can find more details of this formalism and examine its application to a more complex problem in HH15 . In this paper we shall use it for the development of quantum thermodynamics at strong coupling, starting with the third law.

A.1 Influence functional for open quantum systems

Consider for our system a quantum harmonic oscillator (called an Unruh-DeWitt detector in relativistic quantum information) moving along a prescribed spatial trajectory $\mathbf{z}$ in $1+3$ Minkowski spacetime. We can call $\mathbf{z}$ its external degree of freedom, while its internal degree of freedom is the oscillator’s displacement $\chi$. The bath is represented by a massless scalar field $\phi$. The action of the combined system is

[TABLE]

where $S_{\chi}$, $S_{\phi}$, $S_{I}$ are the actions which describe the free quantum oscillators, the bath field and their interaction, given respectively by

[TABLE]

where $x=(t,\mathbf{x})$, $\omega_{b}$ is the bare natural frequency of the oscillator, and an overdot denotes taking the time derivative of a variable. The internal degree of freedom of the quantum oscillator is assumed to be linearly coupled to its bath with coupling strength $e$, which can take on a finite (non-vanishing) value.

Assume that the combined system at time $t=0$ is in a product state 151515For a discussion of the physical consequences of factorizable initial conditions and generalizations, see e.g., HPZ ; JH1 ; PazRom ; FRH ..

[TABLE]

where $\rho^{(\chi)}$ is the initial density operator for the internal degree of freedom, and has a Gaussian form

[TABLE]

The parameter $\varsigma$ is the width of the wavepacket, and $\chi_{i}$, $\chi_{f}$ are shorthands for $\chi$ at the initial time $t=0$ and the final time $t$ respectively, that is, $\chi_{i}=\chi(0)$ and $\chi_{f}=\chi(t)$. This subscript convention will be adopted for the other variables. The bath is initially in its own thermal state at temperature $\beta^{-1}$, with the density matrix

[TABLE]

Here $H^{\phi}[\phi]$ is the free scalar field Hamiltonian associated with the action $S_{\phi}[\phi]$.

The time evolution of the density matrix operator of the combined system is then described by the unitary evolution operator $U(t,0)$ associated with the action (83),

[TABLE]

In the path-integral representation the total density matrix at time $t$ is related to its values at an earlier moment $t=0$ by

[TABLE]

The variables evaluated along the forward and backward time paths, respectively corresponding to $U$ and $U^{-1}$ in (87), will be distinguished by the subscripts $+$, $-$.

A.2 Reduced Density Operator and Coarse-Grained Effective Action

When we are interested only in the dynamics of the system as influenced by the bath, we can work with the reduced density matrix of the system, obtained by tracing out the microscopic degrees of freedom of the bath in the total density matrix, namely,

[TABLE]

where the retarded Green’s function $G_{R}^{(\phi)}$ of the scalar field $\phi$ is defined by

[TABLE]

The Hadamard function is the expectation value of the anti-commutator of the quantum field $\phi$, and is hence temperature dependent. The retarded Green’s function, on the other hand, does not have any temperature dependence. The exponential containing $G_{R}^{(\phi)}$ and $G_{H}^{(\phi)}$ in (89) is called the Feynman-Vernon influence functional $\mathcal{F}$, and we may write it as

[TABLE]

where $S_{IF}$ is called the influence action. It consistently incorporates all the influences of the bath on the system of our interest.

From (89), we may view the time evolution of the reduced density matrix as a map from its initial value $\rho_{\chi}(0)$ to its final value $\rho_{\chi}(t)$ by a superoperator $J$,

[TABLE]

The matrix elements of the superoperator $J$ are expressed by

[TABLE]

where $q$, $r$ are respectively the relative coordinate and the center-of-mass coordinate,

[TABLE]

Here, $S_{CG}$, called the ‘coarse-grained’ effective action, governs the dynamics of the reduced system under the influence of the bath. It contains the actions of the system plus the influence of the bath on the system described by the influence action $S_{IF}$,

[TABLE]

The path integral in the evolutionary operator $J$ can be evaluated exactly because the coarse-grained effective action (96) is quadratic in $q$ and $r$.

A.3 Stochastic Effective Action and Langevin Equations

With the help of the Feynman-Vernon identity, we can express the imaginary part of the coarse-grained effective action $S_{CG}$ in (96) in terms of a Gaussian noise $\xi$,

[TABLE]

with

[TABLE]

Here the angular brackets represent the ensemble average over the probability distribution functional $\mathcal{P}[\xi]$. Thus the exponential of the coarse-grained effective action $S_{CG}$ can be expressed as a distributional integral

[TABLE]

where $S_{SE}$ is the stochastic effective action

[TABLE]

At this point, we may use the stochastic effective action to either derive the Langevin equation, or to construct the stochastic reduced density matrix. We first take the former route.

A.3.1 Langevin Equations

Taking the variation of $S_{SE}$ with respect to $q$ and then letting $q=0$, we arrive at the Langevin equation,

[TABLE]

It describes the time evolution of the reduced system under the non-Markovian influence of the bath. In this case, the influence is manifested in the form of the local stochastic driving noise $\xi$ and the nonlocal dissipative force,

[TABLE]

In general, this nonlocal expression implies the evolution of the reduced system is history-dependent. However, in this one-oscillator example, the retarded Green’s functions matrix has a simple form

[TABLE]

so the Langevin equation reduces to a purely local expression

[TABLE]

where the renormalized frequency $\omega$ is obtained by lumping the divergence term of $G_{R}^{(\phi)}(s,s^{\prime})$ with the original bare frequency $\omega_{b}$, and $\gamma=e^{2}/8\pi m>0$ is the damping constant, which serves as a convenient measure for the system-bath coupling strength. Eq. (103) is seen to be the Langevin equation for a driven, damped oscillator, as anticipated.

The reduced system is superficially non-conservative with the presence of friction and noise forces, which originate from the interaction between the system and its environment. These two processes are, however, connected by the fluctuation-dissipation relation. This relation plays a fundamental role in the energy flow balance between the system and the bath: fluctuations in the bath show up as noise and its backaction on the system gives rise to dissipative dynamics.

The general solution to (101) or (103)

[TABLE]

are expanded in terms of fundamental solution matrices $d_{1}$ and $d_{2}$. They are simply the homogeneous solutions of the corresponding equation of motion but satisfy a particular set of initial conditions,

[TABLE]

This can be the starting point for computing the physical observables of the reduced system. For example, the power delivered by the noise $\xi$ at any given time $t>0$ is given by

[TABLE]

since initially the oscillator’s displacement $\chi$ is not correlated with the noise force $\xi$. Likewise the power output at time $t$ due to the dissipative self-force is given by

[TABLE]

where the velocity uncertainty $\langle\dot{\chi}^{2}(t)\rangle$ takes the form

[TABLE]

if initially $\chi(0)$ and $\dot{\chi}(0)$ are not correlated. The first two terms in (110) exponentially decay with time, so at late times the third term dominates. This reflects the fact that at late times the dynamics of the reduced system is governed by the bath.

A.3.2 Stochastic Reduced Density Matrix

The Langevin equation approach illustrated above to obtain the desired physical quantities associated with the dynamics of the reduced system is less formal, but more flexible and physically intuitive. It is particularly convenient if the quantities at hand involve noise either from the bath or externally introduced. Alternatively, a more systematic and complete approach is by means of the reduced density operator, which accounts for the full quantum dynamics of the reduced system and enforces the operator ordering. The drawback is that since the influence functional does not have explicit dependence on the noise, it is not straightforward to implement it for the cases that explicitly depend on the bath noise. An example is the average power input by the noise shown in (108). Only after invoking the Feynman-Vernon identity could the bath noise be made explicit. Next we will show a way to combine the advantages of these two approaches, by incorporating the noise from the bath in the reduced density matrix.

Let us rewrite the reduced density matrix (93) in terms of the stochastic effective action $S_{SE}$ in (100),

[TABLE]

The $\rho^{(\chi)}(q_{f},r_{f},t;\xi]$ is called the stochastic reduced density matrix, which has explicit dependence on the noise $\xi$ of the bath

[TABLE]

with the stochastic effective action given by (100). In this representation, we see that the reduced system, now driven by a classical stochastic force of the bath, is described by the stochastic density matrix. For each realization of the bath noise, the reduced system evolves to a state described by the density matrix (112). Different realizations make the system end up at different final states with probability given by $\mathcal{P}[\xi]$.

To compute the quantum and stochastic average of a dynamical variable, say, $f(\chi;\xi]$ at time $t$, which depends on both the stochastic variable $\xi$ and the quantum operator $\chi$ of the reduced system, we simply evaluate the trace associated with the system variables and the ensemble average associated with the bath noise,

[TABLE]

The procedure in (113) is understood as follows: for each specific realization of the stochastic noise $\xi$, we first calculate the expectation value of the quantum operator $f(\chi;\xi]$ for the state described by the reduced density operator $\rho^{(\chi)}(t;\xi]$. The obtained result, still dependent on the stochastic variable, will then be averaged over according to the probability distribution $\mathcal{P}[\xi]$ of the noise.

As an example, we will compute the same average power $P_{\xi}$ delivered by the stochastic force $\xi$ from bath as in (108). Once we note that $p=m\dot{\chi}$, the power $P_{\xi}$ is then give by

[TABLE]

where the momentum $p$ canonical to the coordinate $\chi$ is given by

[TABLE]

and the trace over the dynamical variables of the reduced system is defined as

[TABLE]

Since the initial state of the reduced system is a Gaussian state and the stochastic effective action is quadratic in the system’s variables, the final state will remain Gaussian and the corresponding reduced density operator thus can be evaluated exactly. To derive the explicit form of the reduced density matrix, we first evaluate the path integrals in (112),

[TABLE]

where $\mathcal{N}$ is the normalization constant, and can be determined by the unitarity requirement. It is given by

[TABLE]

Note that the mean trajectories $\overline{q}$, $\overline{r}$ are solutions to the Langevin equation (103) with the boundary conditions $\overline{q}(t)=q_{i}$, $\overline{q}(0)=q_{i}$ and $\overline{r}(t)=r_{b}$, $\overline{r}(0)=r_{i}$. Thus they and their time derivatives are functionals of the stochastic noise $\xi$. Explicitly, in terms of the boundary values, we can write $\overline{r}(s)$ as

[TABLE]

for $0\leq s\leq t$. The functions $\mu(s)$, $\nu(s)$ are

[TABLE]

and the current $\mathcal{J}_{r}(s)$ is given by

[TABLE]

Moreover, we can write the partial derivative $\partial/\partial\chi$ as

[TABLE]

Now we are ready to evaluate the power delivered by the stochastic force $\xi$, Eq. (114) becomes

[TABLE]

The expressions in the square brackets can be reduced to

[TABLE]

Thus the power delivered to the system from the bath is equal to

[TABLE]

This is exactly the same as (108).

Appendix B Thermodynamics of simple systems in a common bath

B.1 System of One Harmonic Oscillator

We first examine the one-oscillator system. The absence of mutual interaction between the constituents of the system renders it exactly and fully solvable without approximation. It can readily be compared with the corresponding case in conventional thermodynamics, where the system is extremely weakly coupled to and assumed to be always in thermal equilibrium with a heat bath, whose dynamics is of no concern to the system beyond its being in a thermal state with a temperature parameter. In this case the mechanical energy is given by

[TABLE]

where $\tilde{d}_{2}(\kappa)$ is given by

[TABLE]

The integral in (126) is logarithmically divergent, so we will introduce a regularization scheme in due course. In addition since it is much more difficult to evaluate the integral of the hyper-trigonometric function, we write the factor $\coth\dfrac{\beta\kappa}{2}$ as the summation of the algebraic function of the Matsubara frequency $\nu_{n}$,

[TABLE]

in hope that the resulting integral contains only the algebraic function.

On account of regularization, we may assume it valid to exchange the order of integration and summation, so that (126) becomes

[TABLE]

The evaluation of the integral in (129) is straightforward except for the contribution of the zero mode $n=0$, which needs a separate treatment from the $n\neq 0$ case. The contribution from the $(n=0)$ zero mode has an infrared (IR) divergence. We introduce an IR cutoff $\delta$ and obtain

[TABLE]

For $n\neq 0$, the integral in the summation is well defined and it gives

[TABLE]

Now Eq. (129) becomes

[TABLE]

The summation including infinitely high Matsubara frequencies will give an inevitable UV logarithmic divergence, as expected from (126). We insert a damping factor $e^{-\nu_{n}/\Lambda}$, with $\Lambda>0$, to regularize the summation and arrive at

[TABLE]

where $w_{\pm}=(W\pm i\,\gamma)/2\pi$ and $W$ is the resonance frequency, given by $W=\sqrt{\omega_{p}^{2}-\gamma^{2}}$. The hypergeometric function ${}_{2}F_{1}(a,\,b,\,c;\,z)$ is defined by

[TABLE]

and has a branch-cut on the complex $z$ plane along the real axis from 1 to $\infty$.

In the limit $\Lambda\to+\infty$, the mechanical energy of the oscillator becomes

[TABLE]

where $\mathbb{H}(n)$ is the $n^{\mathrm{th}}$ harmonic number. The cutoff parameter $\Lambda$ defines the highest energy scale in the problem, Its presence can be understood as the consequence that the oscillator couples with a bath that contains a huge number of degrees of freedom. It results from the bath contribution on very short length scales. Since $\ln\Lambda$ is accompanied by the damping constant $\gamma$, the cutoff-dependent term in (135) is negligible for weak oscillator-bath interaction while it can have a significant contribution in the strong interaction limit. Thus the internal energy $E(\infty)$ in principle can depend on the cutoff scale. Note that since the cutoff-dependent term does not depend on temperature, it will not appear in the heat capacity.

This yields for the heat capacity

[TABLE]

when $\gamma<\omega_{p}$. We thus see the heat capacity grows algebraically from zero at low temperature and then saturates to unity at high temperature. The function $\psi^{(n)}(z)$ is the $n^{\mathrm{th}}$ derivative of the digamma function.

Next we examine the weak oscillator-bath coupling limit $\gamma\to 0$. In this limit the mean mechanical energy (135) becomes

[TABLE]

where $\mathcal{O}(\gamma)$ contains the finite contribution of the order $\gamma$. 161616It also contains the cutoff-dependent contributions but they are of the order $\Lambda^{-1}$ and higher. contribution of the order $\gamma$. In principle as $\gamma\to 0$, the physical frequency $\omega_{p}$ will approach to the bare value $\omega_{b}$. This seems innocuous at first sight for the weak oscillator-bath interaction regime, but their values can be drastically different in the strong interaction and the large cutoff scale limit. In particular since $\omega_{p}$ and $\omega_{b}$ are related by (28), the choice of $\gamma$ and $\Lambda$ must be restricted so that $\omega_{p}^{2}$ remains positive definite to prevent unstable dynamics. In fact the physical frequency $\omega_{p}$ is determined by experimental preparation; thus from the operational viewpoint we can let the physical frequency be fixed at the energy scale of measurement. Alternatively, in the action (II), we can assume the system parameters take the physical values, and we introduce counter terms to cancel contributions due to interactions with the bath DC06 . Both approaches eventually produce equivalent results. We can compare results for different values of the system-bath interaction strength $\gamma$.

Following this protocol, the corresponding heat capacity is then given by

[TABLE]

It is divergence-free due to the fact that the logarithmic divergence does not depend on temperature.

Let us now compare this expression with the corresponding heat capacity in conventional thermodynamics. Consider a harmonic oscillator in its thermal state

[TABLE]

Its mean mechanical energy is given by

[TABLE]

It does not depend on the coupling between the oscillator and the bath (because of this, one may not realize that conventional thermodynamics is an open-system theory) and in fact this expression is the same as the dominant term of (137) in the weak coupling limit. For this system the heat capacity in conventional thermodynamics is exactly given by

[TABLE]

and is equal to the leading term of (138) in the weak coupling limit. Note that it has a different low temperature asymptote from that in (136). This is central to the consideration of the third law. Similar results have been obtained for the same configuration, which is nonetheless initially prepared in the equilibrium global thermal state HanIng06 , but the physical contents are different. The similarity in outcomes based on the nonequilibrium initial state and equilibrium global thermal state is not a coincidence, as has been discussed in SFTH . Essentially it is the consequence of the damped, stable motion of the reduced system due to the interaction with the bath.

For nonvanishing oscillator-bath coupling, there is a stark difference between the conventional thermodynamical equilibrium and our open-system nonequilibrium approaches. In general, the results in the open-system framework may have cutoff-dependent contributions, as a result of the huge number of degrees of freedom in the bath. Since this cutoff-dependent term is proportional to the oscillator-bath interaction strength, it tends to be ignored in the weak coupling approximation. The second distinction is related to the observation that there is one more scale $\gamma$ in the open-systems nonequilibrium framework, in addition to $\omega$ and $\beta^{-1}$ already existent in the conventional equilibrium thermodynamics framework. This introduces an additional subtlety in defining the low temperature limit $\beta\omega\gg 1$. In the open-systems nonequilibrium framework, there may be a more stringent criterion such as $\beta\gamma\gg 1$ when $\gamma/\omega<1$, or $\beta f(\gamma,\beta)\gg 1$, where $f(\gamma,\beta)$ is a dimensional function of $\gamma$, $\omega$ and $[f]=L^{-1}$, that is, inversely proportional to the length scale $L$. The presence of this additional scale contributes to different predictions of the heat capacity between the two frameworks in the low temperature regime, as can be seen from (136) and (141).

The low-temperature behavior of the heat capacity in (136) has been argued HZ95 , for the global thermal state case, to be related to the density of the state of the harmonic oscillator in the thermal bath. There it has been shown that if the combined system is initially in the equilibrium global thermal state, then the original discrete energy spectrum of the undamped oscillator will become a continuous one with a unique ground level. That is, the oscillator-bath interaction renders the oscillator a gapless system. This interesting observation has not been proven for a nonequilibrium initial state. It may still apply because the spectrum depends on the effective Hamiltonian (Lagrangian) instead of the prepared initial state.

B.2 System of two coupled oscillators in a common bath

When the system has two or more oscillators, the dynamics of the reduced system becomes much more intricate. Other than their direct coupling, the oscillators also interact with each other indirectly through their shared bath. This indirect influence by one oscillator will propagate in the form of the bath or field disturbance exerted onto all the other oscillators. In turn, more and more subsequent repercussions will be proliferated among the perturbed oscillators. Thus the total effect on the system as a whole depends on the history of each oscillator, leading to very complex evolution. In addition to this, further complication in the interpretation of the results arises from the reduced (environment-influenced) system’s parameter space containing regions where its motion is unstable. We need to identify and exclude this case, then expound the results of the reduced system in an equilibrium state after its motion is fully relaxed.

B.2.1 Dynamics

From (8), the equations of motion for two coupled oscillators are

[TABLE]

where the oscillators 1, 2 are respectively located at $\mathbf{z}_{1}$, $\mathbf{z}_{2}$ so that $\ell=\lvert\mathbf{z}_{1}-\mathbf{z}_{2}\rvert$. It is convenient to reorganize the coupled motion of these two oscillators to an uncoupled motion of a fast mode $\Sigma=(\chi_{1}+\chi_{2})/2$ and a slow mode $\Delta=\chi_{1}-\chi_{2}$,

[TABLE]

with $\omega^{2}_{\pm}=\omega^{2}_{p}\pm\sigma$. Although in appearance the variables $\Sigma$ and $\Delta$ satisfy separate equations (144) and (145), it does not mean that their motions are decoupled. This is because the righthand sides of (144) and (145) indicate that the noises from two locations get mingled together. This clearly shows that the two oscillators are correlated due to arbitration of the ambient quantum field bath.

Their solutions are most easily found if we perform the Laplace transformation over this set of equations of motion and turn them into a simultaneous set of algebraic equations

[TABLE]

where $\xi_{+}=(\xi_{1}+\xi_{2})/2$ and $\xi_{-}=\xi_{1}-\xi_{2}$. Unstable motion occurs when the solutions to

[TABLE]

have the positive real parts. The solution that corresponds to unstable (runaway) motion can be shown to be always located on the positive real axis of the complex $z$ plane. Since (148) is not a simple algebraic equation, it will have an infinite number of solutions, symmetrically distributed on both sides of the real axis, except for the previously mentioned runaway solution. It reflects the mutual undulant disturbance mediated by the field from each oscillator. Thus we expect the spatial non-Markovianity renders the motion of two oscillators much more intricate than that of one oscillator. Finding the solutions to (148) is then nontrivial but since the details about locating the perturbative or the asymptotic solutions have been discussed in HHPRD ; HWL08 ; BC63 ; LinHu09 , they will not be repeated here.

However, a word of caution about the choice of $\ell$: When $\ell$ is extremely small such that $2\gamma/\ell\gg 1$ the contribution from the retardation term is comparable with the frequency renormalization due to the interaction of the oscillator with the bath. Therefore the expression for the equation of motion like (148) becomes dubious in the sense that 1) a point particle model is not always feasible in the context of the self-force, as was long pointed out by Rohrlich and others roh , 2) the equation of motion has a different damping term, proportional to the third-order time derivative, instead of the first-order one HHPRD . Thus the effect of finite size of the oscillator must be taken into consideration.

In summary, instability of motion occurs when the formal effective oscillating frequencies of at least one of the two modes

[TABLE]

become negative. It is likely to happen when 1) the oscillating frequencies of the normal mode $\omega_{\pm}$ become imaginary, and/or 2) the non-Markovian field-induced effect becomes too extreme. Finally we remark that the ratio $\varsigma=\sigma\ell/(2\gamma)$ measuring the relative strength between the direct inter-oscillator coupling and the indirect environment-induced non-Markovian effect is a useful quantity for this consideration, as was introduced in HHPRD for expounding the competing physical mechanisms determining the quantum entanglement between two coupled oscillators in a shared bath.

B.2.2 Internal energy, Heat capacity

In the matrix notations (15) for this two-oscillator system, we have

[TABLE]

where $\ell=\lvert\mathbf{z}_{1}-\mathbf{z}_{2}\rvert$, and from (10), we can find $\tilde{G}_{R}^{(\phi)}(\omega;\,r)$ given by

[TABLE]

In the diagonal elements of $\dfrac{e^{2}}{m}\,\mathbf{G}_{R}^{(\phi)}(\kappa)$, the divergent or cutoff-dependent term will be absorbed with the diagonal elements $\omega_{b}^{2}$ in $\boldsymbol{\Omega}_{b}^{2}$ to form the physical frequency $\omega^{2}_{p}$. The remaining term in the diagonal elements of $\dfrac{e^{2}}{m}\,\mathbf{G}_{R}^{(\phi)}(\kappa)$ will then give $i\,2\gamma\kappa$. Thus $\tilde{\mathbf{d}}_{2}(\kappa)$ in (17) becomes

[TABLE]

This is essentially what we need to compute the total mechanical energy (29) of two oscillators when their motion reaches equilibrium after relaxation.

We have assumed somewhat artificially that $\sigma$ is a constant independent of the separation between the two oscillators. We may relax this restriction to allow it to be a function of $\ell_{ij}=\lvert\mathbf{z}_{i}-\mathbf{z}_{j}\rvert$, namely,

[TABLE]

Here $f(z)$ is a monotonically decreasing function of $z$, except for the case $z=0$, where we require $f(0)=0$, i.e, no self-interaction.

The total mechanical energy for the two coupled oscillator system is then given by

[TABLE]

and the corresponding heat capacity is

[TABLE]

Let us examine the sign of the generic expression

[TABLE]

for all $\kappa$ because it will determine the sign of the heat capacity. Explicitly it takes the form

[TABLE]

which is obviously positive for all $\kappa$. Note that the integrand in (B.2.2) is nowhere negative and the integral is well defined. Thus the heat capacity $C(\infty)$ in (B.2.2) is always positive for all nonzero temperatures even with the presence of bath-induced non-Markovian effects between the oscillators. We observe that in the low temperature limit $\beta\to\infty$, the factor

[TABLE]

goes to zero. Thus it implies that the heat capacity must vanish at zero temperature since the integral in (B.2.2) is regular. This argument may be too simplistic. In fact, when $\beta\to\infty$, the major contributions to the integral in (B.2.2) comes from the interval $\lvert\kappa\rvert<\mathcal{O}(\beta^{-1})$, within which the rest of the integrand is slowly varying. Thus we may pull the slowly changing component out of the integral and write (B.2.2) as

[TABLE]

where $K$ is a finite positive constant independent of $\beta$. Thus indeed the heat capacity vanishes algebraically fast as the temperature approaches zero. Most important of all, the heat capacity vanishes for all permissible choices of oscillator separation, inter-oscillator coupling, and oscillator-bath interaction strength.

The analytical expression of the heat capacity for the two-oscillator system cannot be given without resort to approximation due to the retarded non-Markovian effect. However its low-temperature behavior is expected to be more complicated than its one oscillator counterpart because of additional scale $\ell$ dependence. A numerical example is given in Fig. 3. In particular in Fig. 3-(c), we compare the temperature dependence of heat capacity of this system between the conventional thermodynamics and the present nonequilibrium approaches.

An interesting feature appears when the interaction strength between the oscillators and the bath is sufficiently strong. The heat capacity in the two-oscillator system does not increase monotonically as a function of $\gamma$ as is the case in the one-oscillator system. In fact, as shown in Fig. 4, when the damping constant is greater than a critical value $\gamma_{c}$, the heat capacity will increase monotonically afterwards with increasing $\gamma$. This non-monotonic behavior is most easily seen at high temperatures. In this limit the critical value of $\gamma_{c}$ is approximately given by $\sigma\ell/4$ but slightly falls off with lower temperature. This explains why the non-monotonic behavior of the heat capacity is not seen at low temperatures.

The non-monotonicity in the heat capacity results from the presence of the term $\frac{2\gamma}{\ell}\,e^{i\,\kappa\ell}$ in (B.2.2), which in turn corresponds to the retarded terms in the equations of motion (142) and (143). That is, the non-monotonicity in the heat capacity is due to the strong nonlocal effects mediated by the bath. The critical value $\gamma_{c}=\sigma\ell/4=\gamma\varsigma/2$ has a special physical significance, namely, it indicates the relative significance between the direct coupling and the indirect bath-induced causal influence among the oscillators.

When $\gamma>\gamma_{c}$, i.e., $\varsigma<\mathcal{O}(1)$, the bath-induced effect dominates, whereas when $\gamma<\gamma_{c}$, i.e., $\varsigma>\mathcal{O}(1)$, the direct inter-oscillator coupling wins over. This can be further seen from Fig. 5. In Fig. 5–(a) the curves of heat capacity at high temperature $\beta\to 0$ are plotted for different choices of damping constants $\gamma$ with respect to the oscillator separation $\ell$, while in (b) they are drawn for various selections of direct coupling strengths $\sigma$. They all show a monotonic behavior once $\ell$ is smaller than the critical separation $\ell_{c}=4\gamma/\sigma$. Similar behaviors also appear in thermal entanglement HHPLB15 .

Mathematically it is not difficult to understand why this non-monotonicity in heat capacity is readily seen at high temperatures. From (B.2.2), we immediately see that

[TABLE]

and since the remaining terms in the integrand does not depend on temperature $\beta^{-1}$, it implies that the temperature will not have any effect on the contributions from the causal influence $\frac{2\gamma}{\ell}\,e^{i\,\kappa\ell}$. Physically the rise of heat capacity with increasing $\gamma$ at short inter-oscillator separation can be understood as a consequence of increasing thermal fluctuations of the scalar-field bath, since heat capacity is a measure of the variance of internal energy at least in conventional equilibrium thermodynamics 171717In fact it is used as a criterion for the validity of canonical ensemble and the thermodynamic stability of the system.

After the system equilibrates its dynamics is governed by the noise force, thus large quantum or thermal fluctuations from the bath can in principle induce large fluctuations in the energy of the system, which shows up as large heat capacity. To be specific, the Hadamard function of the scalar field (11) depicts the correlation of the noise forces in (142) and (143), and Fig. 6 shows the corresponding power spectrum, which is the temporal Fourier transform of the Hadamard function. In general the spectrum of the scalar-field bath shows larger values for shorter inter-oscillator separation $\ell$ and higher bath temperature $\beta^{-1}$. They reflect the mere facts that when we probe into a smaller spatial region, we see larger quantum field fluctuations (simply seen on account of the uncertainty principle), and at higher temperatures thermal fluctuations become more pronounced. These large fluctuations will then produce large variances in the variables of the system via, for example, (20).

This non-monotonic behavior of heat capacity at short distance may not be easy to observe because the transition occurs on a scale smaller than $\ell_{c}$ at moderate temperatures. This scale may be shorter than the minimal separation at which the motion of the system remains stable or fall below the physical size of realistic systems modeled by oscillators.

Appendix C Difference between ONEq and the CGTs Setups

Here we use a simple calculation based on the ubiquitous quantum Brownian motion to illustrate the difference between ONEq and the CGTs setups.

Consider a system of an harmonic oscillator bilinearly coupled to an environment consisting of harmonic oscillators, also referred to as the quantum Brownian motion (QBM) model. The Hamiltonian of the combined system is:

[TABLE]

where the variables in capital are associated with the system, and $M_{0}$, $\Omega$ are the mass and the oscillating frequency of the system oscillator. The parameter $g_{n}$ is the coupling strength between each bath mode $q_{n}$ and the system oscillator $Q$, while $p_{n}$ and $P$ are their respective conjugate momenta. Such a combined system is sufficiently simple to serve our current purpose. At some initial time $t_{i}$ the system and environment are assumed to be uncorrelated. Moreover, the environment is in the thermal state with respect to its isolated Hamiltonian $H_{E}$:

[TABLE]

In the above equation $Z_{E}=\text{Tr}_{E}\left\{e^{-\beta H_{E}}\right\}$. If the environment is very large and the system-environment coupling is such that the system relaxes to a unique steady (or equilibrium) state irrespective of the initial state, we call the environment a thermal reservoir or a heat bath. Under these assumptions of equilibration it has been shown that SFTH the unique steady (or equilibrium) state $\rho_{S}^{eq}$ has the form:

[TABLE]

where $Z_{C}$ is defined analogous to $Z_{E}$ above, and $\rho_{S}(t)$ is the reduced state of the system, evolved out of the initially uncorrelated global state (161). This state is referred to as the equilibrium state to distinguish it from the thermal state given by the Boltzmann-Gibbs distribution with respect to the system Hamiltonian alone.

It is important to note that although the system relaxes to the equilibrium state, the global state of the system+environment is not in equilibrium. In particular,

[TABLE]

The density matrix $\rho_{C}$ here is the state of the combined system and environment, evolved from the initial state (161); thus $\rho_{S}(t)=\operatorname{Tr}_{E}\rho_{C}(t)$. Eq. (163) is true because the thermal state of the combined system+environment is a stationary state of the full Hamiltonian and thus cannot be reached from a nonstationary state under Hamiltonian dynamics. What (162) says is that the reduced system state is consistent with a global thermal state and that no further information about the nonequilibrium state of the combined system+bath can be obtained if one has access to the system only. In other words the information on the non-equilibrium state of the combined system+environment is not stored in the system but rather in the environment and in the correlations between the system and the environment.

To demonstrate this point, we focus our attention to a single bath mode. The solutions to the equation of motion are:

[TABLE]

Here $G(t)$ is the Green’s function of the system oscillator, the exact form of which is not important for this discussion (see Ref. SFTH for details). Under the assumptions of equilibration $G(t)\rightarrow 0$ as $t\rightarrow\infty$. This is the reason why the system cannot keep any of the memory of its initial state. However, since the closed system dynamics is unitary, that information is not lost. Rather, it is distributed over the bath modes. The first two terms in the expression for $q_{k}(t)$ do keep the memory of the initial thermal state of the bath. The last term is the only place where the memory of the initial state of the system survives in the bath.

[TABLE]

The important observation is that although $G(t)$ and $\dot{G}(t)$ are decaying functions, the above integrals are oscillatory. (This can be seen explicitly by choosing $G(t)\propto e^{-\gamma t}$ (Ohmic case) and doing the integrals explicitly.) As a result, the information on the initial state of the system survives forever in the bath degrees of freedom. Similarly, information on initial state also survives in the correlations between the system and bath modes as can be seen by studying the late-time limit of $\langle Q(t)q_{k}(t)\rangle$. However, the coupling to each mode is extremely small under the assumption of equilibration, and this makes it practically impossible to extract this information in reality.

Under the assumption of equilibration, the reduced system relaxes to the equilibrium state at late times, irrespective of the initial state of the system. As a result any thermodynamic formulation that only relies on the reduced state of the system for its definitions, will be independent of the initial state of the system. Moreover, all the information about the reduced state of the system can be obtained by assuming a global thermal state without introducing any errors.

In the weak coupling limit, $|H_{I}|\rightarrow 0$, the equilibrium state approaches the thermal state. But in general, the equilibrium state differs from the thermal state. It is common to define the so called “Hamiltonian of mean force” to quantify this difference. In essence, Hamiltonian of mean force is the operator with respect to which the equilibrium state has the Boltzmann-Gibbs form, where the temperature is that of the bath.

In conclusion, a) for strong coupling between the system and bath, if the system can approach the equilibrium state, then the reduced density matrix of the open system is the same as the reduced density matrix in the CGTs framework upon integrating out the bath. b) However, the final global states are different, despite the fact in both cases the dynamics is generated by the same Hamiltonian ($H_{C}=H_{S}+H_{I}+H_{B}$). This is because two global systems start out with different initial states, and the unitary evolution does not change the distinguishability between the states. c) Note that in the CGTs set up, even though the closed system is assumed to be in a global thermal state, the system is not necessarily in a thermal state. One needs extra assumptions, such as the system is very weakly coupled to the bath. It is under the same condition that the open system can reach thermalization. One should be mindful of these presumptions when comparing the thermal states in both setups. We mention these facts as a cautionary note to remind ourselves and for those who want to compare quantities calculated in these three different setups.

Appendix D Difference between Equilibrium States and Thermal States in Classical Thermodynamics

To help understand the important differences between equilibrium states and thermal states, we give a short yet general derivation below where one can see how their differences can be quantified.

Consider a classical system ($S$) described with position and momentum variables $x=(r,p)$ in the phase space. The system might consist of many particles in arbitrarily many dimensions. In that case $x$, $r$, $p$ are vectors. The thermal state of this system is described by the following probability distribution on the phase space

[TABLE]

where $H_{S}(x)$ is the system Hamiltonian. Now assume the system is coupled to an environment whose position and momenta we denote by $y=(r^{\prime},p^{\prime})$. The Hamiltonian of the environment is $H_{E}(y)$ and the interaction Hamiltonian is $H_{I}(x,y)$. Consider the thermal state of the combined system plus environment. The probability distribution on the combined phase space is given by:

[TABLE]

If we are only interested in the state of the system ($S$), we integrate this distribution over the environment variables $y$. The result is the equilibrium distribution for the system alone:

[TABLE]

It is conventional to lump the second term on the righthand side into the exponent by defining the potential of mean force as

[TABLE]

where upon the probability density of the system becomes

[TABLE]

It is clear that $H(x)\neq H^{*}(x)$ in general. As a result, the thermal state differs from the equilibrium state classically as well as quantum mechanically.

A more rigorous way to quantify the difference between equilibrium and thermal distributions for classical systems is to use the Kullback–Leibler divergence which is a measure of distinguishability for probability distributions. Its quantum analog is the quantum relative entropy.

Appendix E Causality Issues with the $M$-oscillator Bath

Here we add some cautionary comments on the differences between an oscillator bath which has been used extensively in Brownian motion studies, e.g., FeyVer ; CalLeg ; GSI88 ; HPZ and the scalar-field bath, which also has been used by many, e.g, LebSpo77 ; UnrZur89 ; LinHu09 ; HHPRD . We show, with the assumption of two or more system oscillators placed in the same spatial location, e.g., in CHY ; PazRec which is often assumed in the $M$-oscillator bath situation, there is a causality issue. No such issues exist with system oscillators in spatially separate locations in a scalar-field bath.

Recall how one proceeds from a harmonic oscillator representation to field theory HM94 . Consider the plane-wave expansion of a massless scalar quantum field

[TABLE]

The creation and annihilation operators $a^{\dagger}_{\mathbf{k}}(t)$, $a^{\vphantom{\dagger}}_{\mathbf{k}}(t)$ satisfy an equation of motion similar to that of the ladder operators for the quantum harmonic oscillators

[TABLE]

with $\omega=\lvert\mathbf{k}\rvert$. In particular when $\mathbf{x}=0$, Eqs. (172) and (173) together represent a collection of quantum harmonic oscillators. The Hamiltonian of the massless scalar field in the plane-wave expansion (172)

[TABLE]

is to be compared to the Hamiltonian of a collection of quantum harmonic oscillators,

[TABLE]

where $m$ is the mass and $\chi_{i}$ is the displacement of $i^{\text{th}}$ oscillator, whose natural frequency is given by $\omega_{i}=i\,\Delta$. The parameter $\Delta$ is some fundamental mode frequency. The ladder operators $a^{\vphantom{\dagger}}_{i}$, $a^{\dagger}_{i}$ in the latter case are defined by

[TABLE]

where $p_{i}$ is the canonical momentum conjugate to $\chi_{i}$. For this and other reasons most people would not make a distinction between two bath models. For example, this is so if one does not ask where each of the system oscillators is located. If they are stacked up at one particular spatial location or when the dipole approximation $e^{i\,\mathbf{k}\cdot\mathbf{x}}\simeq 1$ is applicable, then there is no difference whether the bath is described by a bunch of oscillators or a massless quantum scalar field, since one can arbitrarily shift the location of the system such that $\mathbf{x}=\mathbf{z}_{i}=0$. However, there are still fundamental differences between these two models.

In quantum field theory, the creation and annihilation operators $a_{\mathbf{k}}^{\vphantom{\dagger}}$, $a_{\mathbf{k}}^{\dagger}$ depend on the mode functions we use to expand the quantum field. In the plane wave expansion case, the mode function takes a simple form $e^{i\,\mathbf{k}\cdot\mathbf{x}}$ so the vacuum annihilated by the corresponding $a_{\mathbf{k}}^{\vphantom{\dagger}}$ and the associated number states have a definite three-momentum $\mathbf{k}$. Thus, even merely using the uncertainty principle argument, one can see that these quanta are highly nonlocalized, distributed over the whole configuration space. By contrast for harmonic oscillators, the ground state and the excited states of each bath oscillator are essentially confined by the corresponding harmonic potential $m\omega^{i}\chi_{i}^{2}/2$. Thus the higher the value of the natural frequency $\omega_{i}$ is, the more localized the corresponding mode.

The difference between a field and a collection of harmonic oscillators shows up, for example, when there is a spatial boundary present. The boundary will modify the mode functions of the field and alter the two-point functions of the bath to recognize the effects of the boundary. It is not obvious how this boundary-induced spatial dependence can be naturally implemented in the bath-oscillator model.

Additional complexity emerges when the system contains more than one oscillator which are spatially separated from one another, since the $M$-oscillator bath model does not have the dynamical degrees of freedom to register the locations of the spatially-separated constituent oscillators in the system unless extra input of the spatial information of the system constituents and how it enters in their interaction with the bath is provided. More often than not, one simply ignores the spatial distances between the system constituents by assuming that this separation is so small that we can essentially view them as being situated at the same spatial location or by taking the dipole approximation. This creates a causality problem. For relativistic quantum fields, the influence of an object at one spacetime point on another object at another spacetime point is affected by the former’s imprint on the field, which propagates causally to the latter when it begins to exert its influence. The dynamics of the latter will in turn trigger a new disturbance, on top of the previous one it received, in the field, which propagates at finite time and exerts its influence on the motion of all the other components. This field-induced interaction among the spatially separated components of the system depends on the history of its constituents and is thus fundamentally non-Markovian in nature (see, the plots in e.g., LinHu09 ). Notably this causal propagation feature is mostly lacking in the oscillator-bath mode if the spatial information of the system constituents is not properly accounted for 181818There may exist a nonlocal dynamics due to the appropriate choices of the system-bath coupling constants, but the motion is at best temporally nonlocal, unlike the field-bath where the dynamics is in general spatio-temporally nonlocal.. The field-bath model being relativistic naturally incorporates the spatial correlation and respects the causality.