Nonequilibrium Steady State in Open Quantum Systems: Influence Action, Stochastic Equation and Power Balance
NNonequilibrium Steady Statein Open Quantum Systems : Influence Action,Stochastic Equation and Power Balance
J.-T. Hsiang
Department of Physics, National Dong Hwa University, Taiwan andCenter for Theoretical Physics, Fudan University, Shanghai, China
B. L. Hu
Joint Quantum Institute and Maryland Center for Fundamental Physics,University of Maryland, College Park, Maryland 20742, USA
Abstract
The existence and uniqueness of a steady state for nonequilibrium systems(NESS) is a fundamental subject and a main theme of research in statisticalmechanics for decades. For Gaussian systems, such as a chain of harmonicoscillators connected at each end to a heat bath, and for anharmonic oscilla-tors under specified conditions, definitive answers exist in the form of proventheorems. Answering this question for quantum many-body systems poses achallenge for the present. In this work we address this issue by deriving thestochastic equations for the reduced system with self-consistent backaction fromthe two baths, calculating the energy flow from one bath to the chain to theother bath, and exhibiting a power balance relation in the total (chain + baths)system which testifies to the existence of a NESS in this system at late times. Itsinsensitivity to the initial conditions of the chain corroborates to its uniqueness.The functional method we adopt here entails the use of the influence functional,the coarse-grained and stochastic effective actions, from which one can derivethe stochastic equations and calculate the average values of physical variablesin open quantum systems. This involves both taking the expectation values of
Email addresses: [email protected] (J.-T. Hsiang), [email protected] (B. L. Hu)
Preprint submitted to Annals of Physics May 28, 2014 a r X i v : . [ c ond - m a t . s t a t - m ec h ] J un uantum operators of the system and the distributional averages of stochasticvariables stemming from the coarse-grained environment. This method thoughformal in appearance is compact and complete. It can also easily accommodateperturbative techniques and diagrammatic methods from field theory. Takenall together it provides a solid platform for carrying out systematic investiga-tions into the nonequilibrium dynamics of open quantum systems and quantumthermodynamics. Keywords:
Nonequilibrium steady state, Open quantum systems, Influencefunctional formalism, stochastic density matrix, Langevin equation, noise andfluctuations, Energy flow and power balance relations, Quantum transport,quantum thermodynamics
Contents1 Introduction 3 B and S . . . . . . . . . . . . . . 344.1.2 Energy Flow between S and S . . . . . . . . . . . . . . 374.1.3 Energy Flow between S and B . . . . . . . . . . . . . . 384.2 Condition of Stationarity . . . . . . . . . . . . . . . . . . . . . . 394.3 A Mathematically More Concise Derivation . . . . . . . . . . . . 404.4 Steady State Energy Flow at High and Low Temperatures . . . . 434.5 Heat Conductance . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Nonequilibrium stationary states (NESS) play a uniquely important role inmany-body systems in contact with two or more heat baths at different temper-atures, similar in importance to the equilibrium state of a system in contact withone heat bath which is the arena for the conceptualization and utilization of thecanonical ensemble in statistical thermodynamics. The statistical mechanics[1] and thermodynamics [2] of open systems in NESS have been the focus of Defined in a broader sense (A) an open system is one where some of its informationis difficult or impossible to obtain or retrieve, or is coarse-grained away by design or bynecessity, both in theoretical and practical terms, the latter referring to the limited capabilityof the measuring agent or the precision level of instrumentation. The more specific sense(B) used in nonequilibrium statistical mechanics [3] emphasizing the influence of a system’senvironment on its dynamics goes as follows: Start with a closed system comprising of twosubsystems S and S with some interaction between the two, one can express the dynamicsof S including that of S in terms of an integral differential equation. If one subsystem S contains an overwhelmingly large number of degrees of freedom than the other, we call S anenvironment E of S . The influence of E on S is called the backaction. When the environment closed quantum systems can come to equilibrium and thermalize [11].Equilibration of open quantum systems [12] with strong coupling to a heatbath also shows interesting new features [13]. Transport phenomena in open spin systems has also seen a spur of recent activities [14]. Noteworthy in themathematical properties is the role played by symmetry in the nonequilibriumdynamics of these systems [15]. Our current research program on the nonequilibrium dynamics of quantumopen systems attempts to address four sets of issues with shared common basispertaining to NESS: can be characterized by thermodynamic parameters it is called a heat bath at temperature T or a matter reservoir with chemical potential µ , etc. When a great deal of microscopicinformation of S is discarded or coarse-grained, as is the case when it is described only bya few macroscopic parameters, the effect of the environment can be characterized by noiseand fluctuations [4], and their backaction on the system show up in the “reduced” system’sdynamics as dissipation [5], diffusion (quantum diffusion is responsible for the decoherence [6]of quantum phase information). An open system thus carries the influence or the backactionof its environment. Oftentimes the main task in the treatment of open systems is to find theinfluence of the environment on the subsystem. Defined in this sense (B) it is synonymouswith ”reduced” system – with the burden of explanation now shifted to what “reduced” entailsoperationally. . The approach to NESS . Instead of seeking mathematical proofs for thesebasic issues which are of great importance but not easy to come by it is helpfulto see how these systems evolve in time and find out under what conditionsone or more NESS may exist. For this we seek to derive the quantum stochas-tic equations (master, Langevin, Fokker-Planck) for prototypical quantum opensystems (e.g., for two oscillators in contact with two heat baths and extensionto chains and networks) so one can follow their dynamics explicitly, to examinewhether NESS exist at late times, by checking if energy fluxes reach a steadystate and whether under these conditions a energy flow (power) balance relationexist. This is probably the most explicit demonstration of the NESS possible.In addition, the stochastic equations can be used to calculate the evolution ofkey thermodynamic and quantum quantities such as entropy for equilibration/ thermalization considerations and quantum entanglement for quantum infor-mation inquires. B. Quantum transport : Since the seminal paper of [16], the role of nonlinear-ity and nonintegrability in the violation of Fourier law [17] has been exploredin a wide variety of representative classical systems with different nonlinearinteractions, such as the Fermi-Pasta-Ulam (FPU) models [18] or the Frenkel-Kontorova (FK) model [19] and baths of different natures [17, 20, 21]. For theoriginal papers and current status we refer to two nice reviews [20, 21]. For ap-plications of heat conduction to phononics, see [22]. (Note also the recent workon anomalous heat diffusion [23]). Numerical results are a lot more difficult tocome by for quantum many body systems, thus analytic results, even pertur-bative, for weak nonlinearity, are valuable. Finding solutions to the quantumstochastic equations have been attempted for simple systems like a quantumanharmonic oscillator chain coupled to two heat baths at the ends or harmonicoscillators coupled nonlinearly, each with its own heat bath (namely, with orwithout pinning potentials). The related problem of equilibration of open quan-tum systems with nonlinearity remains an open issue. Even at the classical levelthis is not a straightforward issue. The existence of breather modes [24] and5strange’ behavior [25] have been noted. Nonlinearity in quantum system addsa new dimension bearing some similarity or maybe sharing same origins withthe issue of how to decipher scars of classical chaos in corresponding quantumsystems.
C. Fluctuation Relations : Entropy Production in nonequilibrium system andthe role of large deviations in currents; Fluctuation Theorems both in theGallavoti-Cohen vein [26, 27, 28, 29, 30] and the Jarzynski-Crook relations[31, 32, 33, 34]. Much work in this field is formulated in the context of nonequi-librium thermodynamics. The use of microphysics models such as quantumBrownian motion and open quantum systems techniques, including even deco-herence history concepts (for the definition of trajectories), such as used in [35](see references therein) can provide some new perspective into these powerfulrelations.
D. Quantum entanglement at finite temperature [36, 37] It is generally believedthat at high temperatures thermal fluctuations will overshadow quantum en-tanglement. This problem was explored by Audenaert et al [38] who work outexact solutions for a bisected closed harmonic chain at ground and thermalstates, by Anders [39] for a harmonic lattice in 1-3 dimensions and derived acritical temperature above which the quantum system becomes separable. An-ders and Winters [40] further provided proof of theorems and a phase diagramon this issue. Entanglement of a two particle Gaussian state interacting witha single heat bath is investigated recently in [41]. What makes this issue in-teresting is the suggestion [42] that quantum entanglement can persist at hightemperature in NESS. Recently [43] showed by a coupled oscillator model thatno thermal entanglement is found in the high temperature limit. However, theexistence of quantum entanglement in NESS for driven systems as claimed bythe experiment of Galve et al [44], and the calculations for spin systems [45],remains an open issue. We want to settle this issue theoretically, at least forharmonic oscillator systems, with the help of the formalism set up here.6ur first batch of papers will focus on Issues A and B, which we describebelow. A parallel batch will address Issues C and D in later expositions.
The generic quantum open system we study is a simple 1-dimensional quan-tum oscillator chain, with the two end-oscillators interacting with its own heatbath, each described by a scalar field. The two baths combined make up theenvironment. We begin our analysis with two oscillators linearly coupled andexplore whether a NESS exists for this open system at late times. We do this bysolving for the stochastic effective action and the Langevin equations, which ispossible for a Gaussian system. From this we can derive the expressions for theenergy flow from one bath to another through the system. This is the reasonwhy we begin our study with this model, since in addition to its generic charac-ter and versatility, it provides a nice platform for explaining the methodologywe adopt. For the sake of clarity we will work out everything explicitly, so as tofacilitate easier comparison with other approaches. We name two papers whichare closest to ours, either in the model used or in the concerns expressed: thepaper by Dhar, Saito and Hanggi [46] uses the reduced density matrix approachto treat quantum transport, while that of Ghsquiere, Sinayskiy and Petruccione[43] uses master equations to treat entropy and entanglement dynamics.Similar in spirit is an earlier paper by Chen, Lebowitz and Liverani [47] whichuse the Keldysh techniques in a path integral formalism to consider the dissi-pative dynamics of an anharmonic oscillator in a bosonic heat bath, and recentpapers of Zoli [48], Aron et al [49] for instance. The main tools in nonequilib-rium quantum many-body dynamics such as the closed time path (CTP, in-in,or Schwinger-Keldysh) [50] effective action, the two-particle irreducible (2PI)representation, the large N expansion were introduced for the establishment ofquantum kinetic field theory a quarter of centuries ago [51] and perfected alongthe way [52, 53, 54]. Applications to problems in atomic-optical [55], condensedmatter [56], nuclear-particle [57] and gravitation-cosmology [58] have been onthe rise in the last decade. A description of quantum field theoretic methods7pplied to nonequilibrium processes in a relativistic setting can be found in[59]. By contrast, there is far less applications of these well-developed (power-ful albeit admittedly heavy-duty) methodology for the study of nonequilibriumsteady state in open quantum systems in contact with two or more baths. Wemake such an attempt here for the exploration of fundamental issues of nonequi-librium statistical mechanics for quantum many-body systems and to provide asolid micro-physics foundation for the treatment of problems in quantum ther-modynamics which we see will span an increasingly broader range of applicationsin physics, chemistry and biology. Below we explain our methodology and in-dicate its advantage when appropriate, while leaving the details of how it isrelated to other approaches in the sections proper. The mathematical framework of our methodology is the path-integral influ-ence functional formalism [60, 61, 62], under which the influence action [63], thecoarse-grained effective action [64] and the stochastic effective action [65, 66]are defined. The stochastic equations such as the master equation (see, e.g.,[63]) and the Langevin equations (see, e.g., [67]) can be obtained from takingthe functional variations of these effective actions.There are two main steps in this approach we devised:
1. The derivation of the influence action S IF and coarse-grained effective ac-tions S CG for the reduced system (composed of two linearly interacting os-cillators, then extended to a harmonic chain) obtained by coarse-graining orintegrating over the environmental variables (composed of two baths, coupledto the two end oscillators of a harmonic chain). The baths are here representedby two scalar fields [68, 69]. Noise does not appear until the second stage. Thismaterial is contained in Sec. 2. For Gaussian systems the imaginary part of the influence action can be iden-tified via the Feynman-Vernon integral identity with a classical stochastic force(see, e.g., [70, 71]). Expressing the exponential of the coarse-grained effectiveaction S CG in the form of a functional integral over the noise distribution, thestochastic effective action S SE is identified as the exponent of the integrand.8aking the functional variation of S SE yields a set of Langevin equations for thereduced system . Alternatively one can construct the stochastic reduced densitymatrix.
The averages of dynamical variables in a quantum open system includestaking the expectation values of the canonical variables as quantum operatorsand the distributional averages of stochastic variables as classical noises. Weillustrate how to calculate these quantities with both methods in Sec. 3 and 4.Our methodology includes as subcomponents the so-called reduced densitymatrix approach (e.g., [46, 72]), the nonequilibrium Green function (NEGF)[73, 74, 75]), the quantum master equation and quantum Langevin equationapproaches. It is intimately related to the closed-time-path, Schwinger-Keldyshor in-in effective action method, where one can tap into the many useful fieldtheoretical and diagrammatic methods developed. The stochastic equations ofmotion obtained from taking the functional variation of the stochastic effectiveaction enjoy the desirable features that a) they are real and causal, which guar-antee the positivity of the reduced density matrix, and b) the backaction of theenvironment on the system is incorporated in a self-consistent way. These con-ditions are crucial for the study of nonequilibrium quantum processes includingthe properties of NESS. The physical question we ask is whether a NESS exists at late times. Sincewe have the evolutionary equations and their solutions for this system we canfollow the quantum dynamics (with dissipation and decoherence) of physicalquantities under the influence of the environment (in the form of two noisesources). We describe the behavior of the energy flux and derive the balancerelations in Sec. 4.Paper II [77] will treat the same system but allow for nonlinear interactionbetween the two oscillators. For this we shall develop a functional perturbationtheory for treating weak nonlinearity. Entanglement at high temperatures in This could be any of the three kinds mentioned earlier: see e.g., [63, 76, 67] respectivelyfor derivations of the master, the Fokker-Planck (Wigner) and the Langevin equations.
For the description of the dynamics of an open quantum system obtainingthe time development of the reduced density operator pretty much captures itsessence and evolution. We derive the reduced density operator with the influencefunctional and closed-time-path formalisms (for a ‘no-thrill’ introduction, see,e.g., Chapters 5, 6 of [59]).With this reduced density operator one can compute the time evolution ofthe expectation values of the operators corresponding to physical variables inthe reduced system Here we are interested in the energy flux (heat current)flowing between a chain of n identical coupled harmonic oscillators which to-gether represent the system ( S = (cid:80) nk =1 O k ). Let’s call B the bath which O interacts with, and B the bath oscillator O interacts with. Thus B , B areaffectionately named our oscillators’ ‘private’ baths.Writing the reduced density operator in terms of the stochastic effectiveaction + the probability functional, one can compute the energy current betweeneach oscillator and its private bath in the framework of the reduced densityoperator. This functional method provides a useful platform for the constructionof a perturbation theory, which we shall show in the next paper, in treatingweakly nonlinear cases.Alternatively, from the influence action one can derive the Langevin equationdescribing the dynamics of the reduced system under the influence of a noiseobtained from the influence functional. This is probably a more intuitive andtransparent pathway in visualizing the energy flow between the system and thetwo baths. In the sense described in Footnote 1, a reduced system is an open system whose dynamicsincludes the backaction of its coarse-grained environment. .3.2. FeaturesThe fundamental solutions which together determine the evolutionary oper-ator of the reduced density operator all have an exponentially decaying factor.This has the consequences that(1) the dependence on the system’s initial conditions will quickly become in-significant as the system evolves in time. Because of the exponential decay, onlyduring a short transient period are the effects of initial conditions observable. Atlate time, the behavior of the system is governed by the baths. In other words,for Gaussian initial states, the time evolution of the system is always attractedto the behavior controlled by the bath, independent of the initial conditions ofthe system.(2) the physical variables of interest here tend to relax to – becoming ex-ponentially close to – a fixed value in time. For example, the velocity variancewill asymptotically go to a constant on a time scale longer than the inverse ofthe decay constant. In addition all oscillators O k along the chain have the samerelaxation time scale. The energy currents between B – O , or O k – O k +1 , or O n – B n in general allevolve with time, and will depend on the initial conditions. However, after themotion of the oscillators along the chain is fully relaxed, the energy currentsbetween components approach time-independent values, with the same magni-tude.This time-independence establishes the existence of an equilibrium steadystate . Its insensitivity to the initial conditions of the chain testifies to its unique-ness , the same magnitude ensures there is no energy buildup in any componentof the open system: Heat flows from one bath to another via the intermediaryof the subsystems. To our knowledge, unlike for classical harmonic oscillatorswhere mathematical proofs of the existence and uniqueness of the NESS havebeen provided, there is no such proofs for quantum harmonic systems. It is per-haps tempting to make such an assumption drawing the close correspondencebetween quantum and classical Gaussian systems this is what most authors11acitly assume (e.g. [46]). We have not provided a mathematical proof of theexistence and uniqueness of a NESS for this generic system under study. Whatwe have is an explicit demonstration, drawing our conclusions from solving thedynamics of this system under very general conditions – the full time evolutionof the nonequilibrium open system is perhaps more useful for solving physicalproblems.
1. We have obtained the full nonequilibrium time evolution of the reducedsystem (in particular, energy flow along a harmonic chain between B – O , or O k – O k +1 , or O n – B n in a harmonic chain) at all temperatures andcouplings with arbitrary strength .The formal mathematical expressions of the energy current are given in • Eqs. (4.26), (4.28) and (4.32) as well as (4.53), (4.54) for a two-oscillator chain, and • Eqs. (5.4), (5.5) and (5.9) for an n -oscillator chain,from which we can obtain a profile of energy currents between the com-ponents.2. We have established the steady state value of the energy flux in (5.4), (5.5)and (5.9). Manifest equality and time-independence of these expressionsimplies stationarity. There is no buildup or deficit of energy in any of thecomponents.3. We have demonstrate that the NESS current is independent of the initial(Gaussian) configurations of the chain after the transient period. It thusimplies uniqueness. For comparison, [43] made a weak coupling assumption when working at low temperatures.For Gaussian systems one can solve the full dynamics at least formally in the strong couplingregime – this is well known, see, e.g., [80], and is assumed so in [46]. However, when explicitresults are desired, one often has to make some compromised assumptions, such as weakcoupling between the oscillators and their baths, as done in the last section of [46].
12. We have obtained a Landauer-like formula in • Eq. (4.56) for a two-oscillator chain, and • Eq. (5.10) for an n -oscillator chain5. In particular for the case of two oscillators ( n = 2), we define heat con-ductance (4.61), and have shown that • in the high temperature limit,(a) The steady energy current is proportional to the temperaturedifference between the baths, in (4.60),(b) The heat conductance is independent of the temperature of eitherbath, as is seen in (4.62),(c) The dependence of the conductance on two types of couplingconstants is shown in (4.63) and in Fig. 4.1. • in the low temperature limit(a) the temperature dependence of the steady energy current, (4.64),(4.65) and in Fig. 4.2, and(b) the temperature dependence of the conductivity in (4.67).6. We also plot the general dependence of the NESS energy current on thelength of the chain n in Fig. (5.2), based on our analytical expressions(5.11), (5.12) and (5.13). It shows that • for small n , the NESS current does depend on the length in a non-trivial way; however • for sufficiently large n , the NESS current oscillates but converges toa constant independent of n .
2. Coarse-Grained Effective Action for Open Quantum Systems
Consider a quantum system S = S + S made up of two subsystems S , each consisting of a harmonic oscillator O , interacting with its own bath B ,
13t temperatures T , respectively (assume T > T ). The system by itself isclosed while when brought in contact with heat baths becomes open, owing tothe overwhelming degrees of freedom in the baths which are inaccessible or un-accountable for. The situation of one oscillator interacting with one bath underthe general theme of quantum Brownian motion (QBM) has been studied fordecades and is pretty well-understood, extending to non-Markovian dynamicsin a general environment. Here we wish to extend this study to two such iden-tical configurations, adding a coupling between O and O which is assumedto be linear in this paper and nonlinear in subsequent papers. Let’s call S thecombined system of two coupled quantum Brownian oscillators each interactingwith its private bath. Assume that each oscillator is isolated from the otherthermal bath, thus there is no direct contact between O and B , but there isindirect influence through O ’s coupling to O and its interaction with B . As-sume also that initially the wave functions of the oscillators do not overlap andthat the baths do not occupy the same spacetime region . The physics questionwe are interested in is whether a nonequilibrium steady state exists in S andhow it comes about, in terms of its time evolution. As a useful indicator wewish to describe the energy flow in the three segments: B → O → O → B .It is not clear a priori why energy should flow in a fixed direction (indeedit does not, before each oscillator fully relaxes) and the flow is steady (timetranslational, namely, there is no energy localization or heat accumulation, es-pecially when we extend the two oscillators to a chain). For this purpose weneed to derive the evolution equations for the reduced density operator [62] ofthe system proper, S , after it is rendered open, as a result of tracing over thetwo baths they interact with and including their backaction which shows up asquantum dissipation and diffusion in the equations of motion for the reducedsystem. We do this by way of functional formalisms operating at two levels: 1)at the influence or effective action level, familiar to those with experience of the This can be viewed as an idealization of finite-size bath in the limit that the bath degreesof freedom is sufficiently large and the size of the bath is much larger than the scales associatedwith the oscillator’s motion. z ( i ) (see, e.g., [71]), and its displacement is described by χ ( i ) . The bathsare represented by a massless quantum scalar field φ ( i ) (e.g., [68]) at finite tem-perature. (This is what we refer to as a thermal field, it is a quantum, not aclassical, field, although for Gaussian systems quantum and classical equationsof motions have the same form.) The action of the total system is given by S [ χ, φ ] = (cid:90) t ds (cid:110) (cid:88) i =1 m (cid:104) ˙ χ ( i )2 ( s ) − ω χ ( i )2 ( s ) (cid:105) − mσ χ (1) ( s ) χ (2) ( s ) (cid:111) + (cid:88) i =1 (cid:90) t d x i e i χ ( i ) ( s ) δ ( x i − z ( i ) ( s )) φ ( i ) ( x i )+ (cid:88) i =1 (cid:90) t d x i (cid:2) ∂ µ φ ( i ) ( x i ) (cid:3)(cid:2) ∂ µ φ ( i ) ( x i ) (cid:3) , (2.1)among which we have the actions that describes the oscillators S χ , the bath fields S φ , the interaction between the two oscillators S I and between each oscillatorand its bath S II , respectively, S χ [ χ ( i ) ] = (cid:90) t ds m (cid:104) ˙ χ ( i )2 ( s ) − ω χ ( i )2 ( s ) (cid:105) ,S φ [ φ ( i ) ] = (cid:90) t d x i (cid:2) ∂ µ φ ( i ) ( x i ) (cid:3)(cid:2) ∂ µ φ ( i ) ( x i ) (cid:3) ,S I [ χ (1) , χ (2) ] = (cid:90) t ds (cid:104) − mσχ (1) ( s ) χ (2) ( s ) (cid:105) ,S II [ χ ( i ) , φ ( i ) ] = (cid:90) t d x i e i χ ( i ) ( s ) δ ( x i − z ( i ) ( s )) φ ( i ) ( x i ) . Here we assume that each oscillator is linearly coupled to its own thermal bathwith coupling strength e i , and the oscillators are coupled with each other in the We assume that the coupling strength e i is not so strong as to displace the oscillators. χ (1) − χ (2) ) or χ (1) χ (2) (which are equivalent by a shift in the χ coor-dinate), with an interaction strength denoted by σ . For simplicity without lossof physical contents we let the two oscillators have the same mass m and naturalfrequency ω . We leave the prescribed trajectory z ( i ) ( s ) general here because theposition of the oscillator changes the configuration of the quantum field which inturn affects the other oscillator, aspects which need to be included in quantumentanglement considerations (see, e.g., [81]) and in treating relativistic quantuminformation issues (e.g., [82]). In a later section when we turn to calculating theenergy flow we can safely assume that their external (centers of mass) variablesare fixed in space, and only their internal variables χ enter in the dynamics.Now we assume that the initial state of the total system S at time t = 0 isin a factorizable form . ρ (0) = ρ χ ⊗ ρ β ⊗ ρ β , (2.2)where ρ χ is the initial density operator for the system proper S , consisting oftwo oscillators, with each oscillator described by a Gaussian wavefunction ρ χ ( χ ( i ) a , χ (cid:48) ( i ) a ; 0) = (cid:18) πς (cid:19) / exp (cid:20) − ς (cid:0) χ ( i )2 a + χ (cid:48) ( i )2 a (cid:1)(cid:21) . (2.3)The parameter ς is the width of the wavepacket, and the parameters χ a , χ b arethe shorthand notations for χ evaluated at times t = 0 and t respectively, thatis, χ a = χ (0) and χ b = χ ( t ). This subscript convention will also apply to othervariables. Each bath is initially in its own thermal state at temperature β − i ,so the corresponding initial density matrix is ρ β i ( φ ( i ) a , φ (cid:48) ( i ) a ; 0) = (cid:104) φ ( i ) a | e − β i H φ [ φ ( i ) ] | φ (cid:48) ( i ) a (cid:105) (2.4) H φ [ φ ( i ) ] is the free scalar field Hamiltonian associated with the action S φ [ φ ( i ) ].The density operator of the total system is evolved by the unitary evolutionoperator U ( t, ρ ( t ) = (cid:110) U ( t, ρ (0) U − ( t, (cid:111) . (2.5) For a discussion of the physical consequences of factorizable initial conditions and gener-alizations, see e.g., [63, 65, 83, 84].
16n the path-integral representation the total density matrix at time t is relatedto its values at an earlier moment t = 0 by ρ ( χ ( i ) b , χ (cid:48) ( i ) b ; φ ( i ) b , φ (cid:48) ( i ) b ; t )= (cid:40) (cid:89) i =1 (cid:90) ∞−∞ dχ ( i ) a dχ (cid:48) ( i ) a (cid:90) ∞−∞ dφ ( i ) a dφ (cid:48) ( i ) a (cid:90) χ ( i ) b χ ( i ) a D χ ( i )+ (cid:90) χ (cid:48) ( i ) b χ (cid:48) ( i ) a D χ ( i ) − (cid:90) φ ( i ) b φ ( i ) a D φ ( i )+ (cid:90) φ (cid:48) ( i ) b φ (cid:48) ( i ) a D φ ( i ) − (cid:41) exp (cid:16) (cid:88) i =1 i S χ [ χ ( i )+ ] − i S χ [ χ ( i ) − ] (cid:17) × exp (cid:16) i S I [ χ (1)+ , χ (2)+ ] − i S I [ χ (1) − , χ (2) − ] (cid:17) × exp (cid:16) (cid:88) i =1 i S φ [ φ ( i )+ ] − i S φ [ φ ( i ) − ] (cid:17) × exp (cid:16) (cid:88) i =1 i S II [ χ ( i )+ , φ ( i )+ ] − i S II [ χ ( i ) − , φ ( i ) − ] (cid:17) × ρ χ ( χ ( i ) a , χ (cid:48) ( i ) a ; 0) (cid:89) i =1 ρ β i ( φ ( i ) a , φ (cid:48) ( i ) a ; 0) , (2.6)The subscripts +, − attached to each dynamical variable indicate that thevariable is evaluated along the forward and backward time paths, respectivelyimplied by U and U − in (2.5). When we focus on the dynamics of the oscillators S , accounting for onlythe gross influences of their environments but not the details, we work with thereduced density operator of S obtained by tracing out the microscopic degreesof freedom of its environment, their two baths. We obtain ρ χ ( χ ( i ) b , χ (cid:48) ( i ) b ; t ) = Tr φ (1) Tr φ (2) ρ ( χ ( i ) b , χ (cid:48) ( i ) b ; φ ( i ) b , φ (cid:48) ( i ) b ; t )= (cid:90) ∞−∞ (cid:40) (cid:89) i =1 dχ ( i ) a dχ (cid:48) ( i ) a (cid:41) ρ χ ( χ ( i ) a , χ (cid:48) ( i ) a , t a ) (cid:40) (cid:89) i =1 (cid:90) χ ( i ) b χ ( i ) a D χ ( i )+ (cid:90) χ (cid:48) ( i ) b χ (cid:48) ( i ) a D χ ( i ) − (cid:41) × exp (cid:16) (cid:88) i =1 i S χ [ χ ( i )+ ] − i S χ [ χ ( i ) − ] (cid:17) × exp (cid:16) i S I [ χ (1)+ , χ (2)+ ] − i S I [ χ (1) − , χ (2) − ] (cid:17) × (cid:89) i =1 exp (cid:26) i e i (cid:90) t ds ds (cid:48) (cid:18)(cid:104) χ ( i )+ ( s ) − χ ( i ) − ( s ) (cid:105) G R, β i ( s, s (cid:48) ) (cid:104) χ ( i )+ ( s (cid:48) ) + χ ( i ) − ( s (cid:48) ) (cid:105) + i (cid:104) χ ( i )+ ( s ) − χ ( i ) − ( s ) (cid:105) G H, β i ( s, s (cid:48) ) (cid:104) χ ( i )+ ( s (cid:48) ) − χ ( i ) − ( s (cid:48) ) (cid:105)(cid:19)(cid:27) , (2.7)where the retarded Green’s function G R, β i is defined by G R, β i ( s, s (cid:48) ) = i θ ( s − s (cid:48) ) Tr (cid:16) ρ β i (cid:104) φ ( i ) ( z ( i ) ( s ) , s ) , φ ( i ) ( z ( i ) ( s (cid:48) ) , s (cid:48) ) (cid:105)(cid:17) i θ ( s − s (cid:48) ) (cid:104) φ ( i ) ( z ( i ) ( s ) , s ) , φ ( i ) ( z ( i ) ( s (cid:48) ) , s (cid:48) ) (cid:105) = G R ( s, s (cid:48) ) , (2.8)and the Hadamard function G H, β i by G H, β i ( s, s (cid:48) ) = 12 Tr (cid:16) ρ β i (cid:110) φ ( i ) ( z ( i ) ( s ) , s ) , φ ( i ) ( z ( i ) ( s (cid:48) ) , s (cid:48) ) (cid:111)(cid:17) . (2.9)The Hadamard function is simply the expectation value of the anti-commutatorof the quantum field φ ( i ) , and notice that the retarded Green’s function doesnot have any temperature dependence. The exponential containing G R, β i and G H, β i in (2.7) is the Feynman-Vernon influence functional F , F [ χ + , χ − ] = e i S IF [ χ + ,χ − ] = (cid:89) i =1 exp (cid:26) i e i (cid:90) t ds ds (cid:48) (cid:18)(cid:104) χ ( i )+ ( s ) − χ ( i ) − ( s ) (cid:105) G R, β i ( s, s (cid:48) ) (cid:104) χ ( i )+ ( s (cid:48) ) + χ ( i ) − ( s (cid:48) ) (cid:105) + i (cid:104) χ ( i )+ ( s ) − χ ( i ) − ( s ) (cid:105) G H, β i ( s, s (cid:48) ) (cid:104) χ ( i )+ ( s (cid:48) ) − χ ( i ) − ( s (cid:48) ) (cid:105)(cid:19)(cid:27) , (2.10)where S IF is called the influence action. It captures the influences of the envi-ronment on the system S . The coarse-grained effective action (CG) S CG is made of the influence action S IF from the environment and the actions of the system by S CG [ q ( i ) , r ( i ) ]= (cid:40) (cid:88) i =1 S χ [ χ ( i )+ ] − S χ [ χ ( i ) − ] (cid:41) + S I [ χ (1)+ , χ (2)+ ] − S I [ χ (1) − , χ (2) − ] + S IF [ χ + , χ − ]= (cid:90) t ds (cid:26) (cid:88) i =1 (cid:104) m ˙ q ( i ) ( s ) ˙ r ( i ) ( s ) − mω q ( i ) ( s ) r ( i ) ( s ) (cid:105) (2.11) − mσ (cid:104) q (1) ( s ) r (2) ( s ) + q (2) ( s ) r (1) ( s ) (cid:105)(cid:27) + (cid:88) i =1 e i (cid:90) t ds ds (cid:48) (cid:20) q ( i ) ( s ) G R ( s, s (cid:48) ) r ( i ) ( s (cid:48) ) + i q ( i ) ( s ) G H, β i ( s, s (cid:48) ) q ( i ) ( s (cid:48) ) (cid:21) . (2.12)Here we have introduced the relative coordinate q ( i ) and the centroid coordinate r ( i ) , q ( i ) = χ ( i )+ − χ ( i ) − , r ( i ) = 12 (cid:0) χ ( i )+ + χ ( i ) − (cid:1) . (2.13)18nticipating the oscillator chain treated in a later section, it is convenient tointroduce the vectorial notations by q = q (1) q (2) , r = r (1) r (2) , ΩΩΩ = ω σσ ω , G R ( s, s (cid:48) ) = G R ( s, s (cid:48) ) 00 G R ( s, s (cid:48) ) , G H ( s, s (cid:48) ) = G H, β ( s, s (cid:48) ) 00 G H, β ( s, s (cid:48) ) , and from now on assume the coupling strengths e and e are the same, that is, e = e = e . In so doing the coarse-grained effective action can be written intoa more compact form S CG = (cid:90) t ds (cid:26) m ˙ q T ( s ) · ˙ r ( s ) − m q ( s ) · ΩΩΩ · r ( s ) (2.14)+ e (cid:90) s ds (cid:48) (cid:20) q T ( s ) · G R ( s, s (cid:48) ) · r ( s (cid:48) ) + i q T ( s ) · G H ( s, s (cid:48) ) · q ( s (cid:48) ) (cid:21) . Formally, (2.14) is very general, and can be readily applied to the configurationthat oscillators simultaneously interact with two different thermal baths.Since the coarse-grain effective action S CG governors the dynamics of thesystem S under the influence of the environments, the time evolution of thereduced density matrix can thus be constructed with S CG . We write the reduceddensity matrix as ρ χ ( χ ( i ) b , χ (cid:48) ( i ) b ; t ) = (cid:90) ∞−∞ (cid:40) (cid:89) i =1 dq ( i ) a dr ( i ) a (cid:41) J ( q ( i ) b , r ( i ) b , t ; q ( i ) a , r ( i ) a , ρ χ ( q ( i ) a , r ( i ) a ; 0) , (2.15)where J ( q ( i ) b , r ( i ) b , t ; q ( i ) a , r ( i ) a ,
0) = (cid:40) (cid:89) i =1 (cid:90) q ( i ) b q ( i ) a D q ( i ) (cid:90) r ( i ) b r ( i ) a D r ( i ) (cid:41) exp (cid:110) i S CG [ q ( i ) , r ( i ) ] (cid:111) , (2.16)is the evolutionary operator for the reduced density matrix (from time 0 to t ). The path integral in the evolutionary operator J can be evaluated exactlybecause the coarse-grained effective action (2.14) is quadratic in q and r . Wewon’t pursuit this route in this paper but it will be used later for our study ofnonlinear systems. 19 . Stochastic Effective Action and Langevin Equations We now proceed to derive the stochastic equations and find their solutions
Using the Feynman-Vernon identity for Gaussian integrals we can expressthe imaginary part of the coarse-grained effective action S CG in (2.14) in termsof the distributional integral of a Gaussian noise ξξξ ,exp (cid:20) − e (cid:90) t ds (cid:90) t ds (cid:48) q T ( s ) · G H ( s, s (cid:48) ) · q ( s (cid:48) ) (cid:21) = (cid:90) D ξξξ P [ ξξξ ] exp (cid:20) i (cid:90) t ds q T ( s ) · ξξξ ( s ) (cid:21) , (3.1)with the moments of the noise given by (cid:104) ξξξ ( s ) (cid:105) = 0 , (cid:104) ξξξ ( s ) · ξξξ T ( s (cid:48) ) (cid:105) = e G H ( s, s (cid:48) ) . (3.2)The angular brackets here denote the ensemble average over the probabilitydistribution functional P [ ξξξ ]. Thus we may write the exponential of the coarse-grained effective action S CG in a form of a distributional integral e i S CG [ q , r ] = (cid:90) D ξξξ P [ ξξξ ] exp (cid:20) i (cid:90) t ds (cid:26) m ˙ q T ( s ) · ˙ r ( s ) − m q ( s ) · ΩΩΩ · r ( s )+ q T ( s ) · ξξξ ( s ) + (cid:90) s ds (cid:48) q T ( s ) · G R ( s, s (cid:48) ) · r ( s (cid:48) ) (cid:27)(cid:21) = (cid:90) D ξξξ P [ ξξξ ] e i S SE [ q , r ,ξξξ ] , (3.3)where S SE is the stochastic effective action [65] given by S SE [ q , r , ξξξ ] = (cid:90) t ds (cid:26) m ˙ q T ( s ) · ˙ r ( s ) − m q ( s ) · ΩΩΩ · r ( s ) + q T ( s ) · ξξξ ( s )+ (cid:90) s ds (cid:48) q T ( s ) · G R ( s, s (cid:48) ) · r ( s (cid:48) ) (cid:27) . (3.4)At this point, we may use the stochastic effective action to either derive theLangevin equation, or to construct the stochastic reduced density matrix. Weproceed with the former route below. 20 .2. Langevin Equations Taking the variation of S SE with respect to q and letting q = 0, we arriveat a set of Langevin equation, m ¨ χχχ ( s ) + m ΩΩΩ · χχχ ( s ) − (cid:90) s ds (cid:48) G R ( s, s (cid:48) ) · χχχ ( s (cid:48) ) = ξξξ ( s ) . (3.5)Formally, this equation of motion describes the time evolution of the reducedsystem under the non-Markovian influence of the environment. The influence ismanifested in the form of the local stochastic driving noise ξξξ and the nonlocaldissipative force, (cid:90) s ds (cid:48) G R ( s, s (cid:48) ) · χχχ ( s (cid:48) ) . In general, this nonlocal expression implies the evolution of the reduced system ishistory-dependent. However, in the current configuration, the retarded Green’sfunctions matrix has a very simple form G R ( s, s (cid:48) ) = − e π θ ( s − s (cid:48) ) δ (cid:48) ( s − s (cid:48) ) (cid:32) (cid:33) , (3.6)so the Langevin equation reduces to a purely local form m ¨ χχχ ( s ) + 2 mγ ˙ χχχ ( s ) + m ΩΩΩ R · χχχ ( s ) = ξξξ ( s ) , (3.7)where ΩΩΩ R is obtained by absorbing the divergence of G R ( s, s (cid:48) ) into the diagonalelements of the original ΩΩΩ , and γ = e / πm >
0. We immediately see thateq. (3.7) in fact describes nothing but a bunch of coupled, driven, dampedoscillators. Thus the Langevin equation has a very intuitive interpretation .The general solution to (3.5) or (3.7) can be expanded in terms of fundamen-tal solution matrices D and D . They are simply the homogeneous solutionsof the corresponding equation of motion but satisfy a particular set of initialconditions, D (0) = 1 , ˙ D (0) = 0 , (3.8) Start with two noninteracting oscillators, each interacts with its own thermal bath. Thesystem is described by two decoupled yet almost identical Langevin equations, the only differ-ence is in the noises of the two baths at different temperatures. Now turn on the interactionbetween the two oscillators, then each Langevin equation should acquire an extra force termassociated with the other oscillator’s variable. (0) = 0 , ˙ D (0) = 1 . (3.9)Thus the general solution is given by χχχ ( s ) = D ( s ) · χχχ (0) + D ( s ) · ˙ χχχ (0) + 1 m (cid:90) s ds (cid:48) D ( s − s (cid:48) ) · ξξξ ( s (cid:48) ) . (3.10)This can be the starting point to compute the variance of physical observables,their correlation functions or the variance of the conjugated variables. Forexample, the symmetrized correlation functions of χχχ (where the curly bracketsbelow represent the anti-commutator) are given by12 (cid:104){ χχχ ( t ) · χχχ T ( t (cid:48) ) }(cid:105) = D ( t ) · (cid:104) χχχ (0) · χχχ T (0) (cid:105) · D ( t (cid:48) ) + D ( t ) · (cid:104) ˙ χχχ (0) · ˙ χχχ T (0) (cid:105) · D ( t (cid:48) )+ e m (cid:90) t ds (cid:90) t (cid:48) ds (cid:48) D ( t − s ) · G H ( s, s (cid:48) ) · D ( t (cid:48) − s (cid:48) ) , (3.11)if initially χχχ (0) and ˙ χχχ (0) are not correlated. Notice there that our choice of theparameters m , σ , e and ω renders the fundamental solution matrices symmetric,so we do not explicit show the transposition superscript in the place it is needed.This is a good point to comment on the Langevin equation (3.5) and thederived results such as (3.11). Compared to the equation of motion of a closedsystem the Langevin equation describing the dynamics of a reduced system hastwo additional features, a stochastic forcing term (noise) on the RHS and a dis-sipative term on the LHS. The noise term is a representation of certain measureof coarse-graining of the environment and the backaction of the coarse-grainedenvironment manifests as dissipative dynamics of the reduced system. In theinfluence functional framework the Langevin equation is obtained by taking thefunctional variation of the stochastic effective action S SE about the mean tra-jectory q → ξξξ and thus insensitive to taking the noise distributional average defined by the22robability functional P [ ξξξ ]. If one writes the initial conditions as quantum op-erators of the canonical variables, one may identify the homogeneous part of thecomplete solution as the quantum operators associated with the reduced sys-tem, whose dissipative behavior will in most cases diminish while the reducedsystem relaxes in time. To accommodate both the quantum and the stochasticaspects we only need to extend the meaning of the angular brackets (cid:104)· · · (cid:105) tothat of both taking the expectation value and the distributional average we canproperly incorporate the intrinsic quantum nature of the reduced system andthe noise effects, as demonstrated in (3.11). This approach works very nicelyfor the quantities which take on symmetric ordering, and the result is consis-tent with that computed by the reduced density matrix [85]. It may becomeproblematic if the quantities of interest take on a different ordering from thesymmetric one.Likewise we can find the variances pertinent to the reduced system in theNESS configurations by the Langevin equation. For example, the elements ofthe covariance matrix are given by (cid:104) ∆ χχχ ( l ) b (cid:105) = e m (cid:104)(cid:90) t ds ds (cid:48) D ( s ) · G H ( s − s (cid:48) ) · D ( s (cid:48) ) (cid:105) ll , (3.12) (cid:104) ∆ p ( l ) b (cid:105) = e (cid:104)(cid:90) t ds ds (cid:48) ˙ D ( s ) · G H ( s − s (cid:48) ) · ˙ D ( s (cid:48) ) (cid:105) ll , (3.13)12 (cid:104) { ∆ χχχ ( l ) b , ∆ p ( l ) b }(cid:105) = e m (cid:104)(cid:90) t ds ds (cid:48) D ( s ) · G H ( s − s (cid:48) ) · ˙ D ( s (cid:48) ) (cid:105) ll , (3.14)at late time t (cid:29) γ − . The contributions from the homogenous solutions aretransient, exponentially decaying with time, so they almost vanish on the evo-lution time scale greater than γ − .In the limit t → ∞ , eqs. (3.12)–(3.14) become (cid:104) ∆ χχχ ( l ) b (cid:105) = e m (cid:104)(cid:90) ∞−∞ dω π (cid:101) DDD ( ω ) · (cid:101) G H ( ω ) · (cid:101) DDD ∗ ( ω ) (cid:105) ll , (3.15) (cid:104) ∆ p ( l ) b (cid:105) = e (cid:104)(cid:90) ∞−∞ dω π ω (cid:101) DDD ( ω ) · (cid:101) G H ( ω ) · (cid:101) DDD ∗ ( ω ) (cid:105) ll , (3.16)12 (cid:104) { ∆ χχχ ( l ) b , ∆ p ( l ) b }(cid:105) = 0 , (3.17)23here (cid:101) DDD ( ω ) is given by (cid:101) DDD ( ω ) = (cid:104) − ω I + ΩΩΩ R − i γω I (cid:105) − , ΩΩΩ R = ω R σσ ω R . (3.18)We see that its inverse Fourier transformation DDD ( s ) is the kernel to the Langevinequation, that is, it satisfies (3.7) with an impulse force described by a deltafunction δ ( s ). Furthermore, DDD ( s ) is related to the fundamental solution matrix D ( s ) by DDD ( s ) = θ ( s ) D ( s ).The results in (3.15) and (3.16) cannot be further simplified due to thefact there are two thermal baths with different temperature. The fluctuation-dissipation relation in this case becomes a matrix relation (cid:101) G H ( ω ) = coth β ω β ω · Im (cid:101) G R ( ω ) , (3.19)or in this case (cid:101) G H, β ( ω ) 00 (cid:101) G H, β ( ω ) = coth β ω β ω · Im (cid:101) G R ( ω ) 00 Im (cid:101) G R ( ω ) . Note that each subsystem with its private thermal bath still has its own fluctuation-dissipation relation (cid:101) G H, β i ( ω ) = coth β i ω (cid:101) G R ( ω ) . (3.20)Although the fluctuation-dissipation relation is diagonal, the matrix (cid:101) DDD ( ω ) in(3.15) and (3.16) will blend together the effects of both thermal baths. Forexample, let us examine (cid:104) ∆ χχχ (1) b (cid:105) . In the asymptotic future, it is given by (cid:104) ∆ χχχ (1) b (cid:105) = 1 m (cid:90) ∞−∞ dω π (cid:26)(cid:12)(cid:12) (cid:101) DDD ( ω ) (cid:12)(cid:12) (cid:101) G H ( ω ) + (cid:12)(cid:12) (cid:101) DDD ( ω ) (cid:12)(cid:12) (cid:101) G H ( ω ) (cid:27) (3.21)as seen from (3.15). Physically it is not surprising since each oscillator’s dy-namics needs to reckon with the other oscillator and its bath, albeit indirectly,and thus is determined by both baths.24lternatively, this feature can be seen from the dynamics of the normalmodes of the reduced system that diagonalize ΩΩΩ R . Let v and v be the eigen-vectors of ΩΩΩ R , with eigenvalues ω = ω R + σ and ω − = ω R − σ , respectively.The 2 × U = ( v , v ) can be used to diagonalize the ΩΩΩ R matrix intoΛΛΛ = U T · ΩΩΩ R · U = ω ω − . (3.22)Correspondingly, q and r will be rotated by U to uuu = U T · q , vvv = U T · r . (3.23)If both oscillators are fixed in space, as is the case for our investigation of theNESS, then the Green’s functions looks much simpler because they don’t havespatial dependence. Thus the new Green’s function matrix GGG , transformed by U , is related to the original one G by GGG R ( s, s (cid:48) ) = U T · G R ( s, s (cid:48) ) · U , GGG H ( s, s (cid:48) ) = U T · G H ( s, s (cid:48) ) · U . (3.24)Writing out the matrix GGG explicitly, e.g., taking
GGG H as an example, yields GGG H ( s, s (cid:48) ) = U T · G H ( s, s (cid:48) ) · U = e − · G H, β ( s, s (cid:48) ) 00 G H, β ( s, s (cid:48) ) · − = e G H, β ( s, s (cid:48) ) + G H, β ( s, s (cid:48) ) G H, β ( s, s (cid:48) ) − G H, β ( s, s (cid:48) ) G H, β ( s, s (cid:48) ) − G H, β ( s, s (cid:48) ) G H, β ( s, s (cid:48) ) + G H, β ( s, s (cid:48) ) . (3.25)Again we see that when we decompose the degrees of the freedom of the oscil-lators into their normal modes, the effects from both thermal baths are super-posed.For later reference, we explicitly write down the elements of the matrix (cid:101) DDD ( ω )25s, (cid:101) DDD ( ω ) = (cid:104) − ω I +ΩΩΩ R − i γω I (cid:105) − = − ω + ω R − i γω det (cid:101) DDD − ( ω ) − σ det (cid:101) DDD − ( ω ) − σ det (cid:101) DDD − ( ω ) − ω + ω R − i γω det (cid:101) DDD − ( ω ) , (3.26)with det (cid:101) DDD − ( ω ) = (cid:0) ω − ω R + i γω (cid:1) − σ . Thus the variance? of χ (1) takesthe form (cid:104) ∆ χχχ (1) b (cid:105) (3.27)= 1 m (cid:90) ∞−∞ dω π (cid:0) ω − ω R (cid:1) + 4 γ ω (cid:104) ( ω − ω R + σ ) + 4 γ ω (cid:105) (cid:104) ( ω − ω R − σ ) + 4 γ ω (cid:105) (cid:101) G H ( ω )+ σ (cid:104) ( ω − ω R + σ ) + 4 γ ω (cid:105) (cid:104) ( ω − ω R − σ ) + 4 γ ω (cid:105) (cid:101) G H ( ω ) , and (cid:101) G iiH ( ω ) = ω π coth β i ω . (3.28)Apparently at late time t → ∞ , the displacement variance of O in (3.27)approaches a time-independent constant. In the context of open quantum systems we mention two ways to obtain thedesired physical quantities associated with the dynamics of the reduced system:one is by way of the Langevin equation, which is less formal, more flexible andphysically intuitive. It is particularly convenient if the quantities at hand involvenoise either from the environment or externally introduced. The other is byway of the reduced density operator approach, which is easy to account for theintrinsic quantum dynamics of the reduced system and to enforce the operatorordering. The drawback is that it is less straightforward to use this method tocompute the expectation values of operators corresponding to physical variableswhich contain the environmental noise because the influence functional doesnot have explicit dependence on the noise. Only after invoking the Feynman-Vernon Gaussian integral identity would the noise of the environment, now in26he form of a classical stochastic forcing term, be made explicit. We will showin this section a way to combine the advantages of these two approaches, byincorporating the noise from the environment in the reduced density matrix,whose dynamical equation is obtained by taking the functional variation of thestochastic effective action.Let us rewrite the reduced density matrix (2.15) in terms of the stochasticeffective action S SE in (3.4), ρ χ ( q b , r b ; t ) = (cid:90) ∞−∞ d q a d r a ρ χ ( q a , r a ; 0) (cid:90) q b q a D q (cid:90) r b r a D r exp (cid:110) i S CG [ q , r ] (cid:111) = (cid:90) ∞−∞ d q a d r a ρ χ ( q a , r a ; 0) (cid:90) q b q a D q (cid:90) r b r a D r (cid:90) D ξξξ P [ ξξξ ] e i S SE [ q , r ,ξξξ ] = (cid:90) D ξξξ P [ ξξξ ] ρ χ ( q b , r b , t ; ξξξ ] , (3.29)The term in the integrand ρ χ ( q b , r b , t a ; ξξξ ] is called the stochastic reduced den-sity matrix which is seen to have explicit dependence on the noise ξξξ of theenvironment: ρ χ ( q b , r b , t b ; ξξξ ] = (cid:90) ∞−∞ d q a d r a ρ χ ( q a , r a , t a ) (cid:26)(cid:90) q b q a D q (cid:90) r b r a D r (cid:27) e i S SE [ q , r ,ξξξ ] . (3.30)with the stochastic effective action given by (3.4), S SE [ q , r , ξξξ ] = (cid:90) t ds (cid:26) m ˙ q T ( s ) · ˙ r ( s ) − m q ( s ) · ΩΩΩ · r ( s ) + q T ( s ) · ξξξ ( s )+ (cid:90) s ds (cid:48) q T ( s ) · G R ( s, s (cid:48) ) · r ( s (cid:48) ) (cid:27) . (3.31)In this rendition, the reduced system, now driven by a classical stochastic forceof the environment (by virtue of the Feynman-Vernon transform) as part of theinfluence from the environment, is described by the stochastic density matrix.For each realization of the environmental noise, the reduced system evolves toa state described by the density matrix (3.30). Different realizations make thesystem end up at different final states with probability given by P [ ξξξ ].The reduced system is ostensibly non-conservative with the presence of fric-tion and noise terms, and rightly so, as these effects originate from the inter-action between the system and its environment, which in the case understudy27onsists of two baths. These two processes are, however, constrained by thefluctuation-dissipation relation associated with each bath. This relation plays afundamental role in the energy flow balance between the system and the bath:fluctuations in the environment show up as noise and its backaction on thesystem gives rise to dissipative dynamics. How this relation bears on the prob-lem of equilibration for a quantum system interacting with one bath is prettywell known. We will show explicitly below how this relation underscores theapproach to steady state for systems in nonequilibrium.To compute the quantum and stochastic average of a dynamical variable,say, f ( χχχ ; ξξξ ] at time t , which contains both the stochastic variable ξξξ and thecanonical variables χχχ of the reduced system, we simply evaluate the trace asso-ciated with the system variables and the ensemble average associated with theenvironmental noise, (cid:104) f ( χχχ ; ξξξ ] (cid:105) = (cid:90) D ξξξ P [ ξξξ ] Tr χ ρ χ ( t ; ξξξ ] f ( χχχ ; ξξξ ] . (3.32)The procedure in (3.32) is as follows: for each specific realization of the stochas-tic source, we first calculate the expectation value of the quantum operator f ( χχχ ; ξξξ ] for the state described by the reduced density operator ρ χ ( t ; ξξξ ]. Theobtained result, still dependent on the stochastic variable, will then be averagedover according to the probability distribution P [ ξξξ ] of the noise.As an example, we will compute the power P ξ delivered by the stochasticforce (noise) ξ from Bath 1 to Oscillator 1. The power P ξ is defined by P ξ ( t ) = (cid:104) ξ ( t ) ˙ χ (1) ( t ) (cid:105) , (3.33)and observe that p (1) = m ˙ χ (1) . Thus we have P ξ ( t ) = 1 m (cid:104) ξ ( t ) p (1) ( t ) (cid:105) = − im (cid:90) D ξξξ P [ ξξξ ] (cid:90) ∞−∞ d q b d r b δ ( q b ) ξ ( t ) ∂∂χ (1) b ρ χ ( q b , r b , t ; ξξξ ) , (3.34)where the momentum p ( i ) canonical to the coordinate χ ( i ) is given by p ( i ) = − i ∂∂χ ( i ) , (3.35)28nd the trace over the dynamical variables of the reduced system is defined asTr χ = (cid:90) ∞−∞ d q b d r b δ ( q b ) . (3.36)Since the initial state of the reduced system is a Gaussian state and the stochas-tic effective action is quadratic in the system’s variables, the final state willremain Gaussian and the corresponding reduced density operator thus can beevaluated exactly. To derive the explicit form of the reduced density matrix, wefirst evaluate the path integrals in (3.30), (cid:90) q b q a D q (cid:90) r b r a D r exp (cid:26) i (cid:90) t ds (cid:20) m ˙ q T ( s ) · ˙ r ( s ) − m q ( s ) · ΩΩΩ · r ( s ) + q T ( s ) · ξξξ ( s )+ (cid:90) s ds (cid:48) q T ( s ) · G R ( s, s (cid:48) ) · r ( s (cid:48) ) (cid:21)(cid:27) = N exp (cid:20) i m q Tb · ˙ r b − i m q Ta · ˙ r a (cid:21) , (3.37)where N is the normalization constant, which can be determined by the unitarityrequirement. It is given by N = (cid:16) m π (cid:17) det ˙ µµµ (0) . (3.38)Note that the mean trajectories q , r are solutions to the stochastic Langevinequation (3.7) with the boundary conditions q ( t ) = q b , q (0) = q a and r ( t ) = r b , r (0) = r a . Thus they and their time derivatives are functionals of the stochasticnoise ξξξ . Explicitly, in terms of the boundary values, we can write r ( s ) as r ( s ) = ννν ( s ) · r a + µµµ ( s ) · r b + JJJ r ( s ) , (3.39)for 0 ≤ s ≤ t . The functions µµµ ( s ), ννν ( s ) are defined by µµµ ( s ) = D ( s ) · D − ( t ) , (3.40) ννν ( s ) = D ( s ) − D ( s ) · D − ( t ) · D ( t ) , (3.41)and the current JJJ r ( s ) is given by JJJ r ( s ) = 1 m (cid:90) s ds (cid:48) D ( s − s (cid:48) ) · ξξξ ( s (cid:48) ) − m (cid:90) t ds (cid:48) D ( s ) · D − ( t ) · D ( t − s (cid:48) ) · ξξξ ( s (cid:48) ) . (3.42)29dditionally, we can write the partial derivative ∂/∂χ as ∂∂χ (1) b = ∂∂q (1) b + 12 ∂∂r (1) b . (3.43)Here as an example of the stochastic reduced density matrix approach, we willprovide greater details in the derivation of the power delivered by the stochasticforce ξξξ on O . With these, (3.34) becomes P ξ ( t ) = − im N (cid:90) D ξξξ P [ ξξξ ] (cid:90) ∞−∞ d q b d r b δ ( q b ) (cid:90) ∞−∞ d q a d r a ρ ( q a , r a , ξ ( t ) × (cid:34) ∂∂q (1) b + 12 ∂∂r (1) b (cid:35) exp (cid:20) i m q Tb · ˙ r b − i m q Ta · ˙ r a (cid:21) = − im N (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:90) ∞−∞ d q b d r b δ ( q b ) (cid:90) ∞−∞ d q a d r a ρ ( q a , r a , × (cid:26) i m ˙ r (1) b + i m (cid:20) q Tb · ˙ µµµ ( t ) − q Ta · ˙ µµµ (0) (cid:21) (cid:27) exp (cid:20) i m q Tb · ˙ r b − i m q Ta · ˙ r a (cid:21) = N (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:90) ∞−∞ d r b (cid:90) ∞−∞ d q a d r a ρ ( q a , r a , × exp (cid:26) − i m q Ta · (cid:20) ˙ ννν (0) · r a + ˙ µµµ (0) · r b + ˙ JJJ r (0) (cid:21)(cid:27) × (cid:20) r Ta · ˙ ννν T ( t ) + r Tb · ˙ µµµ T ( t ) + ˙ JJJ Tr ( s ) − q Ta · ˙ µµµ (0) (cid:21) = N (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:90) ∞−∞ d r b (cid:90) ∞−∞ d q a d r a ρ ( q a , r a , × (cid:20) ˙ ννν m ( t ) r ( m ) a + ˙ J (1) r ( t ) −
12 ˙ µµµ T m (0) q ( m ) a + im ˙ µµµ m ( t ) ˙ µµµ − mn (0) ∂∂q ( n ) a (cid:21) × exp (cid:26) − i m q Ta · (cid:20) ˙ ννν (0) · r a + ˙ µµµ (0) · r b + ˙ JJJ r (0) (cid:21)(cid:27) = N (cid:18) πm (cid:19) (cid:18) πς (cid:19) (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:90) ∞−∞ d q a d r a exp (cid:20) − r Ta · r a ς − q Ta · q a ς (cid:21) × exp (cid:26) − i m q Ta · (cid:20) ˙ ννν (0) · r a + ˙ JJJ r (0) (cid:21)(cid:27) × (cid:20) ˙ ννν m ( t ) r ( m ) a + ˙ J (1) r ( t ) −
12 ˙ µµµ T m (0) q ( m ) a + im ˙ µµµ m ( t ) ˙ µµµ − mn (0) ∂∂q ( n ) a (cid:21) δ (2) (cid:2) q Ta · ˙ µµµ (0) (cid:3) = N det ˙ µµµ (0) (cid:18) πm (cid:19) (cid:18) πς (cid:19) (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:90) ∞−∞ d q a d r a δ (2) (cid:0) q a (cid:1) exp (cid:20) − r Ta · r a ς − q Ta · q a ς (cid:21) × exp (cid:26) − i m q Ta · (cid:20) ˙ ννν (0) · r a + ˙ JJJ r (0) (cid:21)(cid:27) (cid:26) ˙ ννν ( t ) · r a + ˙ JJJ r ( t ) −
12 ˙ µµµ T (0) · q a − im ˙ µµµ ( t ) · ˙ µµµ − (0) · (cid:16) − i m (cid:104) ˙ ννν (0) · r a + ˙ JJJ r (0) (cid:105) − ς q a (cid:17)(cid:27) = N det ˙ µµµ (0) (cid:18) πm (cid:19) (cid:18) πς (cid:19) (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:90) ∞−∞ d r a exp (cid:20) − r Ta · r a ς (cid:21) × (cid:104) ˙ JJJ r ( t ) − ˙ µµµ ( t ) · ˙ µµµ − (0) · ˙ JJJ r (0) (cid:105) = N det ˙ µµµ (0) (cid:18) πm (cid:19) (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:104) ˙ JJJ r ( t ) − ˙ µµµ ( t ) · ˙ µµµ − (0) · ˙ JJJ r (0) (cid:105) . (3.44)The expressions in the square brackets can be reduced to˙ JJJ r ( t ) − ˙ µµµ ( t ) · ˙ µµµ − (0) · ˙ JJJ r (0) = 1 m (cid:90) t ds (cid:48) ˙ D ( t − s (cid:48) ) · ξξξ ( s (cid:48) ) . (3.45)Thus the power delivered to Oscillator 1 from Bath 1 is equal to P ξ ( t ) = 1 m (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) (cid:90) t ds (cid:48) ˙ D m ( t − s (cid:48) ) ξ m ( s (cid:48) )= 1 m (cid:90) t ds (cid:48) ˙ D m ( t − s (cid:48) ) (cid:90) D ξξξ P [ ξξξ ] ξ ( t ) ξ m ( s (cid:48) )= e m (cid:90) t ds (cid:48) ˙ D m ( t − s (cid:48) ) G mH ( t − s (cid:48) ) . (3.46)We will see that this is exactly the same as (4.14) except that (4.14) is expressedin terms of the normal modes. Alternatively, we can compare (3.46) with (4.46).They are the same.
4. Energy Transport, Power Balance and Stationarity Condition
As an important application of the formalism developed so far we examinein this section how energy is transported in the combined system S of twooscillators with two baths, in the nature of heat flux, to see whether there isany build-up or localization of energy (the answer is no), or whether there isa balance in the energy flow which signifies the existence of a nonequilibriumsteady state (the answer is yes, with several power balance relations).31 .1. Energy Flow between Components To study the energy transport in the system, as we mentioned in the lastsection, it is physically more transparent to use the Langevin equations (3.7)instead of the matrix form we obtained in the previous section. In this sectionwe will illustrate this approach by deriving the Langevin equations for eachsubsystem, then analyze the heat transfer and energy flux balance relationsfrom them. m ¨ χ (1) ( s ) + 2 mγ ˙ χ (1) ( s ) + mω R χ (1) ( s ) + mσ χ (2) ( s ) = ξ (1) ( s ) , (4.1) m ¨ χ (2) ( s ) + 2 mγ ˙ χ (2) ( s ) + mω R χ (2) ( s ) + mσ χ (1) ( s ) = ξ (2) ( s ) . (4.2)Let’s look at the physics from these equations before delving into the calcula-tions. The motion of any harmonic oscillator in the system, say, O , is alwaysaffected by its neighboring oscillator(s), Oscillator 2 in the present two oscilla-tor case, via their mutual coupling with strength mσ . O is also driven intorandom motion by a stochastic force ξ (1) associated with quantum and thermalfluctuations of Bath 1. Fluctuations in Bath 1, described here by a scalar-field,induces a dissipative force (in general, it is a reactive force) on O . Thus theharmonic oscillator O , other than its own harmonic force, is simultaneouslyacted on by these three seemingly unrelated forces. Certain correlations willbe established among them over time. We will examine the interplay amongthese forces and their roles in energy transport, the power they deliver and thepossible existence of power balance relations, which provide the conditions forthe establishment of a nonequilibrium steady state.First we will find the normal modes of the coupled motion (4.1)–(4.2). Byan appropriate linear combination of the original dynamical variables χ (1) , χ (2) χ + = [ χ (1) + χ (2) ] / , χ − = χ (1) − χ (2) , (4.3)we decouple the motions of O and O and arrive at m ¨ χ + ( s ) + 2 mγ ˙ χ + ( s ) + m ω χ + ( s ) = 12 (cid:2) ξ (1) ( s ) + ξ (2) ( s ) (cid:3) = ξ + ( s ) , (4.4) m ¨ χ − ( s ) + 2 mγ ˙ χ − ( s ) + m ω − χ − ( s ) = ξ (1) ( s ) − ξ (2) ( s ) = ξ − ( s ) (4.5)32here ω ± = ω R ± σ . These are the normal modes of the coupled dynamics, andthey act as two independent driven, damped oscillators. Here we require σ < ω R to avoid any instability in the evolution of the normal modes. Assume the fun-damental solutions to (4.4) and (4.5) are given by d (+) i ( s ), d ( − ) i ( s ) respectivelyand satisfy the initial conditions, d (+)1 (0) = 1 , ˙ d (+)1 (0) = 0 , d (+)2 (0) = 0 , ˙ d (+)2 (0) = 1 , (4.6) d ( − )1 (0) = 1 , ˙ d ( − )1 (0) = 0 , d ( − )2 (0) = 0 , ˙ d ( − )2 (0) = 1 . (4.7)Thus the full solutions to the Langevin equations (4.4) and (4.5) are given by χ + ( s ) = d (+)1 ( s ) χ + (0) + d (+)2 ( s ) ˙ χ + (0) + 1 m (cid:90) s ds (cid:48) d (+)2 ( s − s (cid:48) ) ξ + ( s (cid:48) ) , (4.8) χ − ( s ) = d ( − )1 ( s ) χ − (0) + d ( − )2 ( s ) ˙ χ − (0) + 1 m (cid:90) s ds (cid:48) d ( − )2 ( s − s (cid:48) ) ξ − ( s (cid:48) ) . (4.9)The corresponding solutions to χ (1) ( s ) and χ (2) ( s ) can be obtained by the su-perposition of the normal modes, χ (1) ( s ) = χ + ( s ) + 12 χ − ( s ) , (4.10) χ (2) ( s ) = χ + ( s ) − χ − ( s ) , (4.11)such that χ (1) ( s ) = 12 (cid:104) d (+)1 ( s ) + d ( − )1 ( s ) (cid:105) χ (1) (0) + 12 (cid:104) d (+)1 ( s ) − d ( − )1 ( s ) (cid:105) χ (2) (0)+ 12 (cid:104) d (+)2 ( s ) + d ( − )2 ( s ) (cid:105) ˙ χ (1) (0) + 12 (cid:104) d (+)2 ( s ) − d ( − )2 ( s ) (cid:105) ˙ χ (2) (0)+ 12 m (cid:90) s ds (cid:48) (cid:104) d (+)2 ( s − s (cid:48) ) + d ( − )2 ( s − s (cid:48) ) (cid:105) ξ ( s (cid:48) )+ 12 m (cid:90) s ds (cid:48) (cid:104) d (+)2 ( s − s (cid:48) ) − d ( − )2 ( s − s (cid:48) ) (cid:105) ξ ( s (cid:48) ) , (4.12)and likewise χ (2) ( s ) = 12 (cid:104) d (+)1 ( s ) − d ( − )1 ( s ) (cid:105) χ (1) (0) + 12 (cid:104) d (+)1 ( s ) + d ( − )1 ( s ) (cid:105) χ (2) (0)+ 12 (cid:104) d (+)2 ( s ) − d ( − )2 ( s ) (cid:105) ˙ χ (1) (0) + 12 (cid:104) d (+)2 ( s ) + d ( − )2 ( s ) (cid:105) ˙ χ (2) (0)+ 12 m (cid:90) s ds (cid:48) (cid:104) d (+)2 ( s − s (cid:48) ) − d ( − )2 ( s − s (cid:48) ) (cid:105) ξ ( s (cid:48) )33 12 m (cid:90) s ds (cid:48) (cid:104) d (+)2 ( s − s (cid:48) ) + d ( − )2 ( s − s (cid:48) ) (cid:105) ξ ( s (cid:48) ) . (4.13)These are nothing but the superposition of normal modes in a tethered motion.Now we are ready to compute the power or energy flow or heat transfer betweensubsystems. As we have stressed before, we do not a priori assume the existenceof NESS in this system. The energy flow between the neighboring components ofthe total system is not necessarily time-independent, let alone having the samemagnitude. Rather, we seek to demonstrate the presence of a steady energyflow from one bath to the other after the systems is fully relaxed on a time scale t (cid:29) γ − . B and S The interactions between Subsystem 1 and Bath 1 are summarized in thestochastic force ξ and the dissipative self-force of − mγ ˙ χ (1) , after we coarse-grained the degrees of freedom of B . These two forces mediate the energy flowbetween S and B .The average power delivered to Subsystem 1 by the stochastic force (noise)of Bath 1 is given by P ξ ( t ) = (cid:104) ξ ( t ) ˙ χ (1) ( t ) (cid:105) = 12 m (cid:90) t ds (cid:104) ˙ d (+)2 ( t − s ) + ˙ d ( − )2 ( t − s ) (cid:105) (cid:104) ξ ( t ) ξ ( s ) (cid:105) = e m (cid:90) t ds (cid:104) ˙ d (+)2 ( t − s ) + ˙ d ( − )2 ( t − s ) (cid:105) G β H ( t − s ) y = t − s = e m (cid:90) t dy (cid:104) ˙ d (+)2 ( y ) + ˙ d ( − )2 ( y ) (cid:105) G β H ( y ) , (4.14)with (cid:104) ξ ( t ) ξ ( s ) (cid:105) = e G β H ( t − s ). It tells us the rate at which the energy istransported to Subsystem 1 from Bath 1 by means of the stochastic noise.Since we are particular interested in the existence of NESS, we will payspecial attention to the late-time behavior of the energy transport. In the limit t → ∞ when the motion of Subsystem 1 is fully relaxed and noting that thefundamental solutions d i ( s ) = 0 if s <
0, we write this average power as P ξ ( ∞ ) = e m (cid:90) ∞−∞ dy (cid:104) ˙ d (+)2 ( y ) + ˙ d ( − )2 ( y ) (cid:105) G β H ( y )34 4 πγ (cid:90) ∞−∞ dω π − i ω (cid:104) (cid:101) d (+)2 ( ω ) + (cid:101) d ( − )2 ( ω ) (cid:105) (cid:101) G β H ( ω ) , (4.15)where γ = e / πm , and we have defined the Fourier transformation as f ( t ) = (cid:90) ∞−∞ dω π ˜ f ( ω ) e − i ωt , ⇔ (cid:101) f ( ω ) = (cid:90) ∞−∞ dt f ( t ) e + i ωt , (4.16)so that the convolution integrals are given by (cid:90) ∞−∞ dt f ( t ) g ( t ) = (cid:90) ∞−∞ dω π (cid:101) f ( ω ) (cid:101) g ( − ω ) , (4.17) (cid:90) ∞−∞ dt (cid:90) ∞−∞ dt (cid:48) f ( t ) g ( t (cid:48) ) h ( t − t (cid:48) ) = (cid:90) ∞−∞ dω π (cid:101) f ( ω ) (cid:101) g ( − ω ) (cid:101) h ( ω ) . (4.18)The Fourier transforms (cid:101) d (+)2 ( ω ), (cid:101) d (+)2 ( ω ), and (cid:101) G β H ( ω ) are (cid:101) d ( ± )2 ( ω ) = 1 ω ± − ω − i γω , (4.19) (cid:101) G βH ( ω ) = coth βω (cid:101) G R ( ω ) = ω π coth βω . (4.20)Here and henceforth, we will not make distinction between between DDD i ( t ) and D i ( t ), as in (3.18). They differs only by an unit-step function θ ( t ), that is, DDD i ( t ) = θ ( t ) D i ( t ). Mathematically speaking they are totally different in na-ture; the formal is the inhomogeneous solution to the Langevin equation whilethe latter is the homogeneous solution with a special set of initial conditions.However, in practice, they serve the same purpose to the current case. Thuswhen we refer to the fundamental solution, we use the notations D i ( t ) or d i ( t )for both cases unless mentioned otherwise.The power done by the dissipative force − mγ ˙ χ (1) of Subsystem 1 is P γ ( t ) = − mγ (cid:104) ˙ χ (1) 2 ( t ) (cid:105) = − πγ (cid:90) t ds (cid:90) t ds (cid:48) (4.21) (cid:26)(cid:104) ˙ d (+)2 ( t − s ) + ˙ d ( − )2 ( t − s ) (cid:105)(cid:104) ˙ d (+)2 ( t − s (cid:48) ) + ˙ d ( − )2 ( t − s (cid:48) ) (cid:105) G β H ( s − s (cid:48) )+ (cid:104) ˙ d (+)2 ( t − s ) − ˙ d ( − )2 ( t − s ) (cid:105)(cid:104) ˙ d (+)2 ( t − s (cid:48) ) − ˙ d ( − )2 ( t − s (cid:48) ) (cid:105) G β H ( s − s (cid:48) ) (cid:27) . Here we have ignored the contributions independent of the baths, which willbecome exponentially negligible after the subsystems are fully relaxed. We alsoassumed that both baths are independent of each other and the initial states of35oth subsystems are not correlated with either bath. We see right away that thedissipation power already depends on both reservoirs; the connection of System1 with Bath 2 is established through the coupling mσχ (1) χ (2) between the twosubsystems. Thus in this case that one should not expect that at late times thepower delivered by the stochastic force or noise be exactly in balance with thedissipative power as is the equilibrium case for a single Brownian oscillator incontact with one bath.After the subsystem is fully relaxed, the power done by the self-force becomes P γ ( ∞ ) = − πγ (cid:90) ∞−∞ dω π ω (cid:26)(cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ∗ ˜ G β H ( ω )+ (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105) ∗ ˜ G β H ( ω ) (cid:27) , (4.22)where we have used the convolution integral (4.18) and (4.19). It is instructiveto take a closer look into the expressions in (4.22) that are associated with ˜ G β H of Bath 1. We observe that P (1) γ ( ∞ ) = − πγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ∗ ˜ G β H ( ω )= i πγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) − ˜ d (+)2 ( − ω ) − ˜ d ( − )2 ( − ω ) (cid:105) ˜ G β H ( ω )+ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d ( − )2 ( ω ) ˜ d (+) ∗ ( ω ) + ˜ d (+)2 ( ω ) ˜ d ( − ) ∗ ( ω ) (cid:105) ˜ G β H ( ω )= i πγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ˜ G β H ( ω )+ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d ( − )2 ( ω ) ˜ d (+) ∗ ( ω ) + ˜ d (+)2 ( ω ) ˜ d ( − ) ∗ ( ω ) (cid:105) ˜ G β H ( ω ) , in which we have used˜ d ( ± )2 ( ω ) ˜ d ( ± ) ∗ ( ω ) = − i γω (cid:104) ˜ d ( ± )2 ( ω ) − ˜ d ( ± ) ∗ ( ω ) (cid:105) = − i γω (cid:104) ˜ d ( ± )2 ( ω ) − ˜ d ( ± )2 ( − ω ) (cid:105) , (4.23) (cid:101) G βH ( ω ) = (cid:101) G βH ( − ω ) . (4.24)Now if we combine this contribution P (1) γ ( ∞ ) in the total dissipative power P γ ( ∞ ) with P ξ ( ∞ ), we end up with P ξ ( ∞ ) + P (1) γ ( ∞ ) 36 (cid:26) πγ (cid:90) ∞−∞ dω π − i ω (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ˜ G β H ( ω ) (cid:27) + (cid:26) πγ (cid:90) ∞−∞ dω π i ω (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ˜ G β H ( ω )+ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d ( − )2 ( ω ) ˜ d (+) ∗ ( ω ) + ˜ d (+)2 ( ω ) ˜ d ( − ) ∗ ( ω ) (cid:105) ˜ G β H ( ω ) (cid:27) = − πγ (cid:90) ∞−∞ dω π i ω (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ˜ G β H ( ω )+ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d ( − )2 ( ω ) ˜ d (+) ∗ ( ω ) + ˜ d (+)2 ( ω ) ˜ d ( − ) ∗ ( ω ) (cid:105) ˜ G β H ( ω )= 4 πγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105) ∗ ˜ G β H ( ω ) . (4.25)This implies that after the subsystems are fully relaxed, the net energy transportrate, or power input into Subsystem 1 from Bath 1 is given by P ξ ( ∞ ) + P γ ( ∞ )= P ξ ( ∞ ) + P (1) γ ( ∞ ) + P (2) γ ( ∞ ) (4.26)= − πγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105) ∗ (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) . Note its dependence on ˜ G β H ( ω ) − ˜ G β H ( ω ), which vanishes when there is notemperatures difference between the two baths, β − = β − . S and S Next we consider the energy flow between Subsystems 1 and 2. The powerdelivered by Subsystem 2 to Subsystem 1 is given by P ( t ) = − mσ (cid:104) χ (2) ( t ) ˙ χ (1) ( t ) (cid:105) = − σ m (cid:90) t ds (cid:90) t ds (cid:48) (cid:26)(cid:104) d (+)2 ( t − s ) − d ( − )2 ( t − s ) (cid:105)(cid:104) ˙ d (+)2 ( t − s (cid:48) ) + ˙ d ( − )2 ( t − s (cid:48) ) (cid:105) (cid:104) ξ ( s ) ξ ( s (cid:48) ) (cid:105) + (cid:104) d (+)2 ( t − s ) + d ( − )2 ( t − s ) (cid:105)(cid:104) ˙ d (+)2 ( t − s (cid:48) ) − ˙ d ( − )2 ( t − s (cid:48) ) (cid:105) (cid:104) ξ ( s ) ξ ( s (cid:48) ) (cid:105) (cid:27) + homogeneous terms independent of stochastic forces . (4.27)In the limit t → ∞ , the homogeneous terms vanish and we are left with P ( ∞ ) (4.28)37 i πσγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ∗ (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) . Observe the similarity in form with (4.26) except for the sign difference in thesecond square bracket.Let us now compute the average power P ( t ) = − mσ (cid:104) χ (1) ( t ) ˙ χ (2) ( t ) (cid:105) Sub-system 1 delivers to Subsystem 2 via their mutual interaction. In the samemanners as we arrive at (4.28), we find that at late time t → ∞ , the power P ( ∞ ) is given by P ( ∞ ) = i πσγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ∗ × (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) = − P ( ∞ ) . (4.29)Note this is NOT the consequence of Newton’s third law because P is not thetime rate of work ( = power ) done by the reaction force, namely, the forceSubsystem 2 exerts on Subsystem 1. S and B Finally let us look at the energy flow between Subsystem 2 and Bath 2.First, the average power delivered by the stochastic force ξ from Bath 2 onSubsystem 2 is P ξ ( t ) = (cid:104) ξ ( t ) ˙ χ (2) ( t ) (cid:105) . At late limit t → ∞ , we have P ξ ( ∞ ) = 4 πγ (cid:90) ∞−∞ dω π − i ω (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ˜ G β H ( ω ) . (4.30)Similarly the power delivered by the dissipation force − mγ ˙ χ (2) to this subsys-tem is defined by P γ ( t ) = − mγ (cid:104) ˙ χ (2) 2 ( t ) (cid:105) , and it becomes P γ ( ∞ ) = − πγ (cid:90) ∞−∞ dω π ω (cid:26)(cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105) ∗ ˜ G β H ( ω )+ (cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ∗ ˜ G β H ( ω ) (cid:27) , (4.31)in the limit t → ∞ . Following the procedures that lead to (4.26), we find thatafter the subsystems are fully relaxed, the net power flows into Subsystem 2from Bath 2 is P ξ ( ∞ ) + P γ ( ∞ ) (4.32)38 4 πγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105) ∗ (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) . Compared to the energy flow from B to S , namely, P ξ ( ∞ ) + P γ ( ∞ ) derivedin (4.26) , this carries the opposite sign. This has to be the case for a stationarystate to be established. It says that after the system S = S + S is fully relaxed,the energy which flows into S from Bath 1 at temperature β − is the same asthat out of S into Bath 2 at temperature β − .One last task remains in demonstrating the existence of a NESS: we mustshow that the magnitude of energy flow between the system S and either bath isalso the same as the energy flow between the two subsystems, that is, − P ( ∞ ) = P ξ ( ∞ ) + P γ ( ∞ ). We now show that this is indeed so. Owing to the facts that˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) = − σ ˜ d (+)2 ( ω ) ˜ d ( − )2 ( ω ) , (4.33)˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) = 2 (cid:0) ω R − ω − i γω (cid:1) ˜ d (+)2 ( ω ) ˜ d ( − )2 ( ω ) , (4.34)we can write P ( ∞ ) as P ( t )= i πσγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) + ˜ d ( − )2 ( ω ) (cid:105) ∗ (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) = 16 πγ σ (cid:90) ∞−∞ dω π ω (cid:12)(cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) , (4.35)because the imaginary of the integrand is an odd function of ω .As for P ξ ( ∞ ) + P γ ( ∞ ), we can also show that P ξ ( ∞ ) + P γ ( ∞ )= − πγ (cid:90) ∞−∞ dω π ω (cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105)(cid:104) ˜ d (+)2 ( ω ) − ˜ d ( − )2 ( ω ) (cid:105) ∗ (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) = − πγ σ (cid:90) ∞−∞ dω π ω (cid:12)(cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) , (4.36)with ω ± = ω R ± σ . Eqs. (4.35) and (4.36) indicate that indeed we have − P ( ∞ ) = P ξ ( ∞ ) + P γ ( ∞ ), or P ( ∞ ) + P ξ ( ∞ ) + P γ ( ∞ ) = 0 . (4.37)39his has an interesting consequence. From the Langevin equations (4.1) and(4.2), if we multiply them with ˙ χ (1) and ˙ χ (2) respectively and take their indi-vidual average, we arrive at ddt (cid:104) E (1) k (cid:105) = P + P ξ + P γ , (4.38) ddt (cid:104) E (2) k (cid:105) = P + P ξ + P γ , (4.39)where E ( i ) k is the mechanical energy of each subsystem, E ( i ) k = 12 m ˙ χ ( i ) 2 + 12 mω R χ ( i ) 2 . (4.40)Eq. (4.37) then says that the mechanical energy of each subsystem is conservedwhen the whole system reaches relaxation. In addition, the condition of sta-tionarity P ξ + P γ = − P ξ − P γ (4.41)implies the energy of the whole system will go to a fixed value at late time ddt (cid:104) (cid:88) i =1 , m ˙ χ ( i ) 2 + 12 mω R χ ( i ) 2 + mσ χ (1) χ (2) (cid:105) = 0 , as t → ∞ . (4.42)Eq. (4.41) also says that in the end we must have P + P = 0. This is notobvious when compared with the corresponding closed systems. If there is noreservoir in contact with either subsystem, then although the total energy of thewhole system (internal energy) is a constant value, the mechanical energy in eachsubsystem is not. The energy is transferred back and forth between subsystemsvia their mutual coupling mσ χ (1) χ (2) . Thus in the case of the closed systems, P + P (cid:54) = 0 but oscillate with time. The key difference between an open anda closed system may lie in the fact that, as t → ∞ , the dynamics of an opensystem is determined largely by the reservoirs, at least from the viewpoint ofthe Langevin equation. In the above we sought a balance relation by displaying the energy flowsbetween components explicitly. This has the advantage of seeing the physical40rocesses in great detail and clarity. There is a mathematically more conciseformulation of energy transport which we present here. We adopt a matrixnotation for ease of generalization to the harmonic chain case treated in thefollowing sections.Since Subsystem 2 exerts a force − mσ χ (2) on Subsystem 1, the averagepower delivered by Subsystem 2 to Subsystem 1 is P ( t ) = − mσ (cid:104) χ (2) ( t ) ˙ χ (1) ( t ) (cid:105) = − mσ lim t (cid:48) → t ddτ (cid:48) (cid:104) χ (2) ( t ) χ (1) ( t (cid:48) ) (cid:105) = − mσ (cid:26)(cid:20) ς D ( t ) · ˙ D ( t ) + 12 m ς D ( t ) · ˙ D ( t ) (cid:21) + e m (cid:90) t ds ds (cid:48) (cid:104) D ( t − s ) · G abH ( s − s (cid:48) ) · ˙ D b ( t − s (cid:48) ) (cid:105) (cid:27) , (4.43)where we have used (3.11) and the properties that D i are symmetric. As notedbefore the expressions within the square brackets approach zero at late times,so in the limit τ → ∞ , the power P ( t ) becomes P ( ∞ ) = − e σm (cid:90) ∞−∞ ds ds (cid:48) (cid:104) D ( s ) · G H ( s − s (cid:48) ) · ˙ D ( s (cid:48) ) (cid:105) . (4.44)Recall that since D ij ( s ) = 0 for s <
0, we can extend the lower limit of theintegration to minus infinity. Expressing the integrand by the Fourier transformof each kernel function yields P ( ∞ ) = − e σm (cid:90) ∞−∞ dω π (cid:0) − i ω (cid:1)(cid:104) (cid:101) D ∗ ( ω ) · (cid:101) G H ( ω ) · (cid:101) D ( ω ) (cid:105) . (4.45)It says that the average power delivered by Subsystem 2 ( S ) on Subsystem 1( S ) approaches a constant eventually.Next we examine the corresponding power transfer between the S and itsprivate bath ( B ). The average power delivered by the stochastic force (noise)from B to S is P ξ ( t ) = (cid:104) ξ (1) ( t ) ˙ χ (1) ( t ) (cid:105) = 1 m (cid:90) t ds ˙ D (2)1 a ( t − s ) (cid:104) ξ ( t ) ξ a ( s ) (cid:105) = e m (cid:90) t ds (cid:104) ˙ D ( t − s ) · G H ( t − s ) (cid:105) . (4.46)Hence at late times t → ∞ , the average power P ξ becomes P ξ ( ∞ ) = e m (cid:90) ∞−∞ dω π (cid:0) i ω (cid:1)(cid:104) (cid:101) D ∗ ( ω ) · ˜ G H ( ω ) (cid:105) . (4.47)41ikewise, the average power delivered by the dissipative force in S from thebackaction of B is described by P γ ( t ) = − mγ (cid:104) ˙ χ (1) 2 ( t ) (cid:105) = − e γm (cid:90) t ds ds (cid:48) (cid:104) ˙ D ( s ) · G H ( s − s (cid:48) ) · ˙ D ( s (cid:48) ) (cid:105) . (4.48)Again, we have ignored contributions which are exponentially small at latetimes. The value of P γ in the limit t → ∞ is given by P γ ( ∞ ) = − e γm (cid:90) ∞−∞ dω π ω (cid:104) (cid:101) D ∗ ( ω ) · (cid:101) G H ( ω ) · (cid:101) D ( ω ) (cid:105) . (4.49)Therefore the net energy transfer at late times between B and S is P S = P ξ ( ∞ ) + P γ ( ∞ )= e m (cid:90) ∞−∞ dω π (cid:0) i ω (cid:1) (cid:101) D a ∗ ( ω ) (cid:104) (cid:101) G aH ( ω ) + i γω (cid:101) D b ( ω ) (cid:101) G abH ( ω ) (cid:105) = e m (cid:90) ∞−∞ dω π (cid:0) i ω (cid:1) (cid:101) D a ∗ (cid:104) I b + i γω (cid:101) D b ( ω ) (cid:105) (cid:101) G abH ( ω )= e m (cid:90) ∞−∞ dω π (cid:0) i ω (cid:1) (cid:110)(cid:104) I + i γω (cid:101) D ( ω ) (cid:105) · (cid:101) G H ( ω ) · (cid:101) D † ( ω ) (cid:111) , (4.50)where we have used the symmetric property of G H . We next note that theFourier transform (cid:101) D ( ω ) satisfies (cid:101) D − ( ω ) = ΩΩΩ − ω I − i γω I , ⇒ (cid:104) ΩΩΩ − ω I − i γω I (cid:105) · (cid:101) D ( ω ) = I , ⇒ (cid:104) ΩΩΩ − ω I (cid:105) · (cid:101) D ( ω ) = I + i γω (cid:101) D ( ω ) , (4.51)with I being a 2 × P S = e m (cid:90) ∞−∞ dω π (cid:0) i ω (cid:1) (cid:110)(cid:104) ΩΩΩ − ω I (cid:105) · (cid:101) D ( ω ) · (cid:101) G H ( ω ) · (cid:101) D † ( ω ) (cid:111) . (4.52)From the definition of ΩΩΩ , we see that (cid:2) ΩΩΩ − ω I (cid:3) ab is in fact σ (cid:0) δ a δ b + δ a δ b (cid:1) ,so that (4.52) becomes P S = − e σm (cid:90) ∞−∞ dω π (cid:0) − i ω (cid:1) (cid:104) (cid:101) D ∗ ( ω ) · (cid:101) G H ( ω ) · (cid:101) D ( ω ) (cid:105) , (4.53)42here we have used the symmetry property of (cid:101) D and (cid:101) G H . Compared to (4.45),we see that P S is equal to P ( ∞ ) at late time. It means that if P S >
0, thenthe energy flow from the bath 1 to the subsystem 1 is equal to the energy flowfrom the subsystem 1 to subsystem 2.In addition from (4.45), we can show that at late time P ( ∞ ) = − P ( ∞ )as follows. Since by construction P ( ∞ ) is a real physical quantity, if we takethe complex conjugate of P ( ∞ ) we should return to the very same P ( ∞ ),that is P ( ∞ ) = P ∗ ( ∞ ) = − e σm (cid:90) ∞−∞ dω π (cid:0) i ω (cid:1) (cid:101) D a ( ω ) (cid:101) G abH ( ω ) (cid:101) D b ∗ ( ω ) (4.54)= e σm (cid:90) ∞−∞ dω π (cid:0) − i ω (cid:1) (cid:101) D b ∗ ( ω ) (cid:101) G abH ( ω ) (cid:101) D a ( ω ) = − P ( ∞ ) . Thus we also establish that P S = − P ( ∞ ).To make connection with (4.35), we observe that from (3.26), we can re-late the elements of the fundamental solution matrices with the correspondingfundamental solutions of the normal modes by (cid:101) D ( ω ) = (cid:101) D ( ω ) = 12 (cid:104) (cid:101) d (+)2 ( ω ) + (cid:101) d ( − )2 ( ω ) (cid:105) = (cid:0) ω R − ω − i γω (cid:1) ˜ d (+)2 ( ω ) ˜ d ( − )2 ( ω ) , (cid:101) D ( ω ) = (cid:101) D ( ω ) = 12 (cid:104) (cid:101) d (+)2 ( ω ) − (cid:101) d ( − )2 ( ω ) (cid:105) = − σ ˜ d (+)2 ( ω ) ˜ d ( − )2 ( ω ) , from (4.34)–(4.33). This enable us to write (4.54) as P ( ∞ ) = 16 πγ σ (cid:90) ∞−∞ dω π ω (cid:12)(cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) . (4.55)Thus we recover (4.35). We may define the steady energy flow J by J ≡ P = P ξ ( ∞ ) + P γ ( ∞ )= 16 πγ σ (cid:90) ∞−∞ dω π ω (cid:12)(cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:104) ˜ G β H ( ω ) − ˜ G β H ( ω ) (cid:105) , (4.56)with ˜ G βH ( ω ) = ω π coth βω . (4.57)43n the high temperature limit β i →
0, we have ˜ G β i H ( ω ) ≈ / (2 πβ i ) so that thesteady energy current becomes, when β i ω (cid:28) J (cid:39) γ σ (cid:0) β − − β − (cid:1) (cid:90) ∞−∞ dω π ω (cid:12)(cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12)(cid:12) ∝ ( T − T ) . (4.58)The integral in (4.58) can be exactly carried out, and it is given by (cid:90) ∞−∞ dω π ω (cid:12)(cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12)(cid:12) = 18 γ σ + 4 γ ω R , (4.59)with ω ± = ω R ± σ . Therefore the steady energy current in the high temperaturelimit is given by J = γσ σ + 4 γ ω R ∆ T = γ ∆ T , γω R (cid:28) σ ,σ γω R ∆ T , γω R (cid:29) σ , (4.60)where ∆ T = T − T , for different relative coupling strengths between the sub-systems and the reservoirs. When γ →
0, that is, when the coupling betweenthe subsystems and their baths is turned off, there is no energy flow. Likewise,if there is no coupling between the subsystems σ →
0, the energy flow alsoterminates, as also expected.
We define the thermal conductance K by the ratio of the steady current overthe temperature difference between the reservoirs, K = lim ∆ T → J ∆ T . (4.61)Thus we find that in the high temperature limit βω R (cid:29)
1, the conductance K = γ σ σ + 4 γ ω R , (4.62)becomes independent of temperature but only depends on the parameters σ , γ and ω R . From Fig. 4.1, we see that the conductance monotonically increaseswith the inter-oscillator coupling σ , and gradually approaches the value γ as longas the constraint σ ≤ ω R is still satisfied. On the other hand, when we fix theinter-oscillator coupling, the conductance rises up to a maximum value σ/ ω R .1 0.2 0.3 0.4 0.5 0.60.0050.0100.0150.0200.025 0.5 1.0 1.5 2.00.0050.0100.0150.0200.0250.030 Figure 4.1: variation of the conductance K with respect to the coupling constants σ or γ inthe high temperature limit. at γ = σ/ ω R , and then gradually decreases to zero as the system-environmentcoupling γ increases. From the expression of (cid:12)(cid:12) ˜ d ( ± )2 ( ω ) (cid:12)(cid:12) we note that it tracesout a Breit-Wigner resonance curve with respect to ω . The resonance featureis well-defined only when γ is sufficiently small, that is, γ (cid:28) ω R . The peak islocated at about ω = ( ω R ± σ ) / and the width of the peak is about 2 γ . There-fore for a fixed value of Ω, the inter-oscillator coupling constant σ determinesthe location of the resonance peak, while the system-environment coupling con-stant γ determines the width of the resonance. The integrand (4.56) containsa product of (cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12) (cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12) , which indicates that there are two resonancepeaks at ( ω R ± σ ) / respectively. Hence the distance between these two peaksis (cid:0) ω R + σ (cid:1) − (cid:0) ω R − σ (cid:1) (cid:39) σω R . When the separation of peaks is much greater than the width, it has two distinct,well-defined peaks. If the separation becomes smaller than the characteristicwidth of each peak, σ < γω R , then the two peaks gradually fuse into one peak.This change of the dominant scale is reflected in the behavior of the conductance K in the respective regimes, K = γ , γω R (cid:28) σ ,σ γω R , γω R (cid:29) σ , (4.63)as can be seen from (4.60). 45n the low temperature limit βω ± (cid:29)
1, we may write the Bose-Einsteindistribution factor in (4.56) as˜ G β H ( ω ) − ˜ G β H ( ω ) = ω π ∞ (cid:88) n =1 (cid:2) e − nβ ω − e − nβ ω (cid:3) , and the steady state current becomes J = 8 πγ σ π ∞ (cid:88) n =1 (cid:90) ∞ dω ω (cid:12)(cid:12)(cid:12) ˜ d (+)2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˜ d ( − )2 ( ω ) (cid:12)(cid:12)(cid:12) (cid:104) e − nβ ω − e − nβ ω (cid:105) = 48 πγ σ πω ω − (cid:20) β − β (cid:21) ∞ (cid:88) n =1 n + · · · , = 8 π γ σ (cid:0) ω R − σ (cid:1) (cid:20) β − β (cid:21) . (4.64)The summation over n gives ζ (4) = π /
90. This familiar number comes fromthe higher order expansions of the coth z function in the limit z → ∞ . If wewrite the steady current (4.64) in terms of temperature, we obtain J = 8 π γ σ (cid:0) ω R − σ (cid:1) (cid:0) T − T (cid:1) = 8 π γ σ (cid:0) ω R − σ (cid:1) (cid:0) T ∆ T + T ∆ T (cid:1) , (4.65)where ∆ T = T − T and T = ( T + T ) /
2. We see that the temperaturedependence of the steady current is different from the high temperature limit.In the low temperature limit it is proportional to T − T . However, for fixed T , the current turns out more or less linearly proportional to the temperaturedifference between the reservoirs except for the case ∆ T (cid:39) T , which is equivalentto T < T , where the contributions of the ∆ T terms appreciable. In theregime ∆ T (cid:28) T , the steady state current is J (cid:39) π γ σ (cid:0) ω R − σ (cid:1) T ∆ T , (4.66)and then we may use the definition of the conductance (4.61) to find: K = 32 π γ σ (cid:0) ω R − σ (cid:1) T . (4.67)Apparently in the low temperature limit, the conductance depends on T , whichis the mean temperature of the two reservoirs. In Fig. 4.2, we plot the steady46 .5 1.0 1.5 2.0 2.50.00020.00040.00060.00080.00100.00120.0014 Figure 4.2: variation of the steady state current J with respect to the temperature T of theBath 1 . The temperature difference ∆ T between two reservoirs is fixed. state current as a function of the temperature T of Bath 1, with the temperaturedifference ∆ T fixed, according to (4.56). We see in the high temperature limit,the current approaches a constant, independent of T , T , as long as ∆ T is fixed,which is consistent with (4.60).
5. Harmonic Chain
We now extend the previous results to a one-dimensional chain of n harmonicoscillators. The oscillators at both ends, labelled as O and O n , are attached totheir own pivate baths of respective temperatures T > T n . The remaining n − k = { , , ..., n − } are insulated from these twobaths, and only interact with their nearest neighbors bilinearly with couplingstrength σ .In analogy with the case of two oscillators in the previous sections, here thecolumn matrix χχχ has n entries, so does the row matrix χχχ T = ( χ (1) , χ (2) , · · · , χ ( n ) ).47he matrices ΩΩΩ , I (cid:48) and G are now spanned to n × n matrices,ΩΩΩ = ω σ · · · σ ω σ · · ·
00 . . . ...... . . . 00 · · · σ ω σ · · · σ ω , I (cid:48) = · · ·
00 0 0 0 · · · · · · · · · , and G ( s, s (cid:48) ) = G β ( s, s (cid:48) ) 0 0 0 · · ·
00 0 0 0 · · ·
00 . . . ...... . . . 00 · · · · · · G β n ( s, s (cid:48) ) . (5.1)In addition, we have the same stochastic effective action as in (3.4) except thatnow the stochastic force is a column vector ξξξ T = ( ξ (1) , , · · · , , ξ ( n ) ) and itsmoments satisfy the Gaussian statistics (cid:104) ξξξ ( s ) (cid:105) = 0 , (cid:104) ξξξ ( s ) · ξξξ T ( s (cid:48) ) (cid:105) = G H ( s, s (cid:48) ) . (5.2)Note that although the matrix G in (5.1) is not invertible, it does not prevent usfrom writing down the stochastic effective action. In fact we can do it by com-ponents and then write them back into the tensor notation. Note the stochasticaverage is defined only with respect to the first and the final components of the(vectorial) stochastic force ξξξ .Taking the variation of the stochastic effective action with respect to q andletting q = 0, we arrive at the Langevin equation, m ¨ χχχ ( s ) + 2 mγ I (cid:48) · ˙ χχχ ( s ) + m ΩΩΩ R · χχχ ( s ) = ξξξ ( s ) . (5.3)where ΩΩΩ R is the same as ΩΩΩ in the structure except that we replace ω by ω R due to renormalization. The derivation of the energy currents between the components of the totalsystem is similar to those presented in Sec. 4.3. Thus the net energy flow from48ath 1 ( B ) to Oscillator 1 ( O ) at late times is given by J = P ξ ( ∞ ) + P γ ( ∞ )= 8 γ (cid:90) ∞−∞ dω ω (cid:104) (cid:101) D ( ω ) · I (cid:48) · (cid:101) D ∗ ( ω ) · (cid:101) G TH ( ω ) − (cid:101) D ( ω ) · (cid:101) G H ( ω ) · (cid:101) D † ( ω ) (cid:105) = 8 γ (cid:90) ∞−∞ dω ω (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) (cid:104) (cid:101) G H ( ω ) − (cid:101) G nnH ( ω ) (cid:105) . (5.4)Here we note that the prefactor (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) in the integrand of (5.4) has depen-dence on the length of the chain n . We will take a closer look at its behaviorlater.To demonstrate the existence of a steady state for the present configuration,we have to show that the steady current between the neighboring oscillatorsis the same as J . Let the late-time energy current flow in the intermediateoscillators from O k to O k +1 be J k,k +1 with k = 2 , , . . . , n −
1. From (4.45), weknow it is given by J k,k +1 = − i e σm (cid:90) ∞−∞ dω π ω (cid:104) (cid:101) D ( ω ) · (cid:101) G H ( ω ) · (cid:101) D † ( ω ) (cid:105) k,k +1 (5.5)= − i γσ (cid:90) ∞−∞ dω ω (cid:104) (cid:101) D k, ( ω ) (cid:101) D k +1 , ∗ ( ω ) (cid:101) G H ( ω ) + (cid:101) D k,n ( ω ) (cid:101) D k +1 ,n ∗ ( ω ) (cid:101) G nnH ( ω ) (cid:105) . Eq. (5.5) does not have any reference to time, so the energy current from O k to O k +1 is also a time-independent constant.To show the equality between (5.4) and (5.5), we will relate (cid:101) D k, ( ω ) (cid:101) D k +1 , ∗ ( ω )or (cid:101) D k,n ( ω ) (cid:101) D k +1 ,n ∗ ( ω ) to (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) so that we can factor out the noise kernel (cid:101) G H ( ω ) in (5.5). These relations are provided in Appendix A. Using the resultsin (7.27) and (7.30) (cid:101) D k, (cid:101) D k +1 , ∗ = + i cσ (cid:12)(cid:12) (cid:101) D n (cid:12)(cid:12) + · · · , (5.6) (cid:101) D n − k, ∗ (cid:101) D n − k +1 , = − i cσ (cid:12)(cid:12) (cid:101) D n (cid:12)(cid:12) + · · · , (5.7)where . . . denotes terms which will have vanishing contributions to the integral(5.5), we can rewrite (5.5) as J k,k +1 = 8 γ (cid:90) ∞−∞ dω ω (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) (cid:104) (cid:101) G H ( ω ) − (cid:101) G nnH ( ω ) (cid:105) , (5.8)49ecalling that c = 2 γω . We immediately see that for any neighboring oscillators O k and O k +1 located along the chain, the energy current J k,k +1 between themis exactly the same as the current J transported from Bath 1 to Oscillator O in (5.4).As a final touch, we compute the energy current from Bath B n , located atthe opposite end of the chain, to Oscillator O n . From (5.4), we find this currentis given by J n = 8 γ (cid:90) ∞−∞ dω ω (cid:104) (cid:101) D ( ω ) · I (cid:48) · (cid:101) D ∗ ( ω ) · (cid:101) G TH ( ω ) − (cid:101) D ( ω ) · (cid:101) G H ( ω ) · (cid:101) D † ( ω ) (cid:105) nn = − γ (cid:90) ∞−∞ dω ω (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) (cid:104) (cid:101) G H ( ω ) − (cid:101) G nnH ( ω ) (cid:105) . (5.9)It has the same magnitude as J in (5.4) and J k,k +1 in (5.8), but opposite insign, which says that the current flows from Oscillator O n to Bath B n , as is alsoexpected.In summary, given a quantum harmonic oscillator chain, where each oscil-lator interacts with its nearest neighbors via bilinear coupling, if the two end-oscillators of the chain are placed in contact with two thermal baths of differenttemperatures, while the oscillators in between are kept insulated from thosebaths, we have explicitly shown that after a time when all the oscillators havefully relaxed, the energy flow along the chain becomes independent of time andthe currents between the neighboring oscillators are the same in both magnitudeand direction J NESS = 8 γ (cid:90) ∞−∞ dω ω (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) (cid:104) (cid:101) G H ( ω ) − (cid:101) G nnH ( ω ) (cid:105) . (5.10)This implies that a NESS exists for a quantum harmonic oscillator chain anda steady current flows from the high temperature front along the chain to thelow temperature end. There is no buildup or localization of energy at anysite along the chain. We emphasize that in the transient phase before theconstituent oscillators come to full relaxation, additional contributions fromthe homogeneous solutions of the oscillators’ modes render the current betweenneighboring oscillators unequal, but in the course of the order of the relaxation50ime, the energy along the chain is re-distributed to a final constant value whilethe whole system settles down to a NESS. We have shown that after the motion of the constituents of the chain reachesrelaxation, a steady thermal energy current exists flowing along the harmonicchain across from the hot thermal reservoir to the cold one. However we havenot addressed the scaling behavior of the steady current with the length ofthe chain. To shed some light on this problem, we first analyze the prefactor (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) , which is proportional to the transmission coefficient in the Landauerformula.From (7.21), we have (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) = σ n − | θ n | , (5.11)where | θ n | = f c + 2 f c + f + 2 σ n − c > . (5.12)Once f , f , f is found we can derive the analytic expression of | θ n | . Indeedwe have shown this in the Appendix, where the general expression of f k is givenby f k = µ n − k +11 − µ n − k +12 µ − µ , (5.13)with the roots of the characteristic equation given by µ = a + √ a − σ , µ = a − √ a − σ , (5.14)we immediately see that when a < σ , that is, when ω lies in the interval √ Ω − σ < ω < √ Ω + 2 σ , two roots µ , µ are complex-conjugated. It inturn implies that θ n , as well as (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) , will be highly oscillatory with ω inthis interval. Analytically (cid:12)(cid:12) θ n (cid:12)(cid:12) is described by | θ n | = σ n sin ψ (cid:20) sin ( n + 1) ψ + 2 c σ (cid:16) nψ (cid:17) + c σ sin ( n − ψ (cid:21) , (5.15)with ψ being ψ = tan − √ σ − a a . (5.16)51 ! Figure 5.1: The generic structure of ω (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) . We choose n = 10 for example. We also argue in the Appendix that outside the interval √ Ω − σ < ω < √ Ω + 2 σ , the transmission coefficient falls to a vanishingly small value veryrapidly. Therefore, generically speaking ω (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) will have a comb-likestructure within the band √ Ω − σ < ω < √ Ω + 2 σ , as is shown in Fig 5.1.The lower envelope of the comb structure is traced by sin ψ , while the upper en-velope is determined by the subleading term in (5.15) because the maxima occurapproximately at the locations where sin ( n + 1) ψ vanishes. In addition, recallthat c = 2 γω , so the upper envelope of the comb structure is almost constant.The number of the spikes of the comb structure is equal to the number of theoscillators in the harmonic chain. However, since the bandwidth of the intervalis independent of n , the width of each spike will scale as n − . This implies, asis argued in the Appendix, that the contribution from each spike to the integralover ω in (5.10) also scales as n − . The argument supplied in the Appendix im-proves with growing n , so we see for sufficiently large n , the scaling behavior ofthe contribution of each spike to the steady current will be nicely counteractedby the increasing number of spikes with the band √ Ω − σ < ω < √ Ω + 2 σ .Therefore it implies that the NESS energy current, in the case of the harmonicchain, is independent of the length of the chain. We show in Fig. 5.2 an exam-52
10 15 20 25 3058.558.658.758.858.959.059.1
Figure 5.2: The scaling behavior of the NESS energy current J for a particular set of param-eters. The number n ranges from 3 to 30. ple of the scaling behavior of the NESS energy current along the chain. Thelength of the chain ranges from n = 3 to n = 30. We see that the value quicklyconverges to an almost n -independent constant.While this qualitative behavior for a perfect quantum harmonic chain (withno defect, impurity or nonlinearity) may be well-known or can be reasoned out,we have hereby provided an explicit quantitative proof of it.
6. Concluding Remarks
The setup of an open system interacting with two heat baths serves as thebasis for a wide range of investigations in physics, chemistry and biology. Theexistence of a nonequilibrium steady state in such a system is an issue of fun-damental importance because, to name just one, it is the pre-condition fornonequilibrium thermodynamics, which serves as a powerful springboard forinvestigations in many areas of sciences and engineering. The existence anduniqueness of NESS have been studied for classical systems decades ago withrigorous mathematical proofs. We want to do the same now for quantum opensystems, starting with the simpler case of continuous variable harmonic systems53in contradistinction to discrete variables such as spin chains, a subject whichhas seen many flowering results).The key findings of this investigation have been enumerated in the Introduc-tion so there is no need to repeat them here. A few general concluding remarkswould suffice.The broader value of this work as we see it is twofold: 1) A demonstrationof the existence of a NESS for the system of interest. Rather than constructingmathematical proofs we provide the full dynamics of the system and derive ex-plicit expressions for the energy flow in each component leading to a proof thatan energy flow balance relation exists. 2) Presenting a toolbox whereby onecan derive the stochastic equations and calculate the average values of physicalvariables in open quantum systems – this involves both taking the expectationvalues of quantum operators of the system and the distributional averages ofstochastic variables originating from the environment. The functional methodwe adopt here has the advantage that it is compact and powerful, and it caneasily accommodate perturbative techniques and diagrammatic methods devel-oped in quantum field theory to deal with weakly nonlinear open quantumsystems, as we will show in a sequel paper. The somewhat laborious and expos-itory construction presented in this paper is necessary to build up a platformfor systematic investigations of nonequilibrium open quantum systems, someimportant physical issues therein will be discussed in future communications.To expand the second point somewhat, our approach is characterized by twoessential features; a) we use a microphysics model of generic nature, namelyhere, a chain of harmonic oscillators interacting with two baths described by twoscalar fields at different temperatures. b) this allows us to derive everything from first principles , e.g., starting with an action principle describing the interactionof all the microscopic constituents and components in the model. This wayof doing things has the advantage that one knows the physics which goes intoall the approximations made, in clearly marked stages. For Gaussian systems,namely, bilinear coupling between harmonic oscillators and with baths, one cansolve this problem exactly, providing the fully nonequilibrium evolution of the54pen system with the influences of it environment (the two heat baths here)accounted for in a self-consistent manner.Self-consistency is an absolutely essential requirement which underlies thecelebrated relation of Onsager, for example (one may refer to the balance rela-tions we obtained here as the quantum Onsager relations) and realization of thesymmetries in the open systems which were used for the mathematical proofs ofNESS. This consistency condition is not so well appreciated in the open quantumsystem literature, but we see it as crucial in the exploitation of the symmetryprinciples mentioned above as well as in treating physical processes when mem-ory effects (in non-Markovian processes) and when the effects of backaction areimportant. This includes situations when one wants to a) treat strongly corre-lated systems or systems subjected to colored noises b) design feedback controlof quantum systems c) engineer an environment with sensitive interface withthe open quantum system, to name a few.
Acknowledgment
This work began in the summer of 2013 when both authorsvisited Fudan University’s Center for Theoretical Physics at the invitation ofProf. Y. S. Wu. Earlier that year BLH visited the group of Prof. BaowenLi at the National University of Singapore. Thanks are due to them for theirwarm hospitality. The leitmotiv to understand nonequilibrium energy transportbegan when BLH attended a seminar by Prof. Bambi Hu at Zhejiang Universityin 2010. He also thanks Prof. Dhar for useful discussions. JTH’s research issupported by the National Science Council, R.O.C.
7. Appendix
In this Appendix, we will derive the relations between (cid:101) D k, ( ω ) (cid:101) D k +1 , ∗ ( ω )or (cid:101) D k,n ( ω ) (cid:101) D k +1 ,n ∗ ( ω ) to (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) . They are used in Sec. 5 to set up theequality of the energy current between the neighboring sites along the chain.To this goal, we first establish some useful relations between the elementsof the fundamental solution matrix (cid:101) D . From the definition of the fundamental55olution, (cid:104) − ω I − i γω I (cid:48) + ΩΩΩ (cid:105) · (cid:101) D ( ω ) = I , (7.1)we see that the matrix (cid:101) D − ( ω ) = (cid:104) − ω I − i γω I (cid:48) + ΩΩΩ (cid:105) (7.2)is symmetric with respect to the diagonal and the anti-diagonal, its inverse, thefundamental solution matrix (cid:101) D , also has these properties, that is, (cid:101) D j,k = (cid:101) D k,j , (cid:101) D j,k = (cid:101) D n +1 − k,n +1 − j . (7.3)In particular, it implies (cid:101) D j, = (cid:101) D ,j = (cid:101) D n +1 − j,n , (7.4)and (5.5) becomes J j,j +1 = − i γσ (cid:90) ∞−∞ dω ω (cid:104) (cid:101) D j, ( ω ) (cid:101) D j +1 , ∗ ( ω ) (cid:101) G H ( ω )+ (cid:101) D n − j, ∗ ( ω ) (cid:101) D n − j +1 , ( ω ) (cid:101) G nnH ( ω ) (cid:105) . (7.5)Now the problem reduces to identifying a relation between (cid:101) D j, ( ω ) (cid:101) D j +1 , ∗ ( ω )or (cid:101) D n − j, ∗ ( ω ) (cid:101) D n − j +1 , ( ω ) and (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) .To find the explicit expressions for the elements of the fundamental matrix (cid:101) D , we can invert (7.2). Since (cid:101) D − forms a tridiagonal matrix, its inverse canbe given by the recursion relations (cid:101) D j, = (cid:0) − (cid:1) j +1 σ j − Υ j +1 θ n , j > , (7.6)where Υ n = a n , Υ n +1 = 1 andΥ j +1 = a j +1 Υ j +2 − σ Υ j +2 , Υ j ∈ C , (7.7)with j = 1 , , . . . , n − j ∈ C . The quantity θ n is the determinant ofthe inverse of the fundamental matrix, that is, θ n = det (cid:101) D − , and satisfies therecursion relation θ k = a k θ k − − σ θ k − , (7.8)56ith θ = 1, θ − = 0 and θ k ∈ C .We introduce two shorthand notations a (cid:48) an a by assigning a (cid:48) ≡ a = a n = − ω − i γω + ω R , and a ≡ a = · · · = a n − = − ω + ω R . We note that a (cid:48) is a complex number and its imaginary part is an odd function of ω . It thenproves useful to express Υ j +1 explicitly in terms of a (cid:48) . The motivation behindthis lies in the fact that inside the square brackets of (7.5), only terms that areodd with respect to ω can have nontrivial contributions to the current. On theother hand, the only source that may contribute to the odd power in ω is theimaginary part of a (cid:48) .From the recursion relation (7.7), we see that in general the variable Υ j canbe expanded by f j Υ j = f j a (cid:48) − f j +1 σ , (7.9)where f j is an ( n − j )–order polynomial of a with f n = 1, f n +1 = 0, and itsatisfies the recursion relation, f j = a f j +1 − σ f j +2 , f j ∈ R (7.10)Here we write down the first a few entries in the sequence { f j } , f n +1 = 0 , f n = 1 ,f n − = a , f n − = a − σ ,f n − = a − aσ , f n − = a − a σ + σ ,f n − = a − a σ + 3 aσ , f n − = a − a σ + 6 a σ − σ . (7.11)Although f j follows a similar recursion relation to (7.7), introduction of f j makesit easier to identify the imaginary part of Υ j . Therefore once we find the generalsolution of f j via the recursion relation (7.10), we will have Υ j , which is usefulto construct the general expression for the elements of the fundamental solutionmatrix (cid:101) D .The recursion relation (7.10) can be solved if we substitute f k ∝ µ − k (7.12)57nto (7.10), we find that µ will satisfy a characteristic equation µ − a µ + σ = 0 . (7.13)The two solutions, labeled by µ and µ , µ = a + √ a − σ , µ = a − √ a − σ , (7.14)are assumed distinct, and then the general solution for f k is given by f k = p µ − k + p µ − k . (7.15)We can use the conditions f n = 1, f n +1 = 0 to fix the unknown coefficients p and p , f n = p µ − n + p µ − n = 1 , (7.16) f n +1 = p µ − n − + p µ − n − = 0 , (7.17)so they are p = µ n +11 µ − µ , p = − µ n +12 µ − µ . (7.18)Thus the general solution of f k takes the form f k = µ n − k +11 − µ n − k +12 µ − µ = n − k (cid:88) m =0 µ n − k − m µ m . (7.19)With the help of these results, we can set up relations between (cid:101) D k, ( ω ) (cid:101) D k +1 , ∗ ( ω )or (cid:101) D n − k, ∗ ( ω ) (cid:101) D n − k +1 , ( ω ) and (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) .To gain some insight, we first write down (cid:12)(cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12)(cid:12) by (7.6). Since (cid:101) D n = (cid:0) − (cid:1) n +1 σ n − Υ n +1 θ n = (cid:0) − (cid:1) n +1 σ n − θ n , (7.20)we find (cid:12)(cid:12) (cid:101) D n (cid:12)(cid:12) = σ n − | θ n | . (7.21)Next, (cid:101) D s, (cid:101) D s +1 , ∗ can be given by (cid:101) D s, (cid:101) D s +1 , ∗ = (cid:20)(cid:0) − (cid:1) s +1 σ s − Υ s +1 θ n (cid:21) (cid:20)(cid:0) − (cid:1) s +2 σ s Υ ∗ s +2 θ ∗ n (cid:21) − σ s − Υ s +1 Υ ∗ s +2 | θ n | . (7.22)We can expand the product Υ s +1 Υ ∗ s +2 by (7.9), and getΥ s +1 Υ ∗ s +2 = (cid:104) f s +1 a (cid:48) − f s +2 σ (cid:105)(cid:104) f s +2 a (cid:48)∗ − f s +3 σ (cid:105) = − (cid:16) f s +1 f s +3 a (cid:48) + f s +2 a (cid:48)∗ (cid:17) σ + · · · , (7.23)where . . . are terms that will not contribute to the integral (7.5). Now recallthat the imaginary part of a (cid:48) is odd with respect to ω , so we write a (cid:48) explicitlyas a (cid:48) = a − i c , where a = − ω + ω R has been defined before while c is equalto 2 γω . In so doing, we are able to condense (7.23) further to highlight itsdependence on the imaginary part of a (cid:48) , that is, c ,Υ s +1 Υ ∗ s +2 = i c (cid:16) f s +1 f s +3 − f s +2 (cid:17) σ + · · · , (7.24)The expression in the parentheses can be evaluated by (7.19), and we find f s +1 f s +3 − f s +2 = − σ n − s − , (7.25)where we have used the fact that µ µ = σ in (7.13) at the final step. Thuseq. (7.24) becomes Υ s +1 Υ ∗ s +2 = − i c σ n − s − + · · · . (7.26)We put it back to (7.22) and arrive at (cid:101) D s, (cid:101) D s +1 , ∗ = + i c σ n − | θ n | + · · · = + i cσ (cid:12)(cid:12) (cid:101) D n (cid:12)(cid:12) + · · · , (7.27)where we have compared the result with (7.21). Again the dots represent termsthat do not contribute to the integral (7.5).Next we proceed to evaluate (cid:101) D n − s, ∗ (cid:101) D n − s +1 , . By (7.6), we obtain (cid:101) D n − s, ∗ (cid:101) D n − s +1 , = (cid:20)(cid:0) − (cid:1) n − s +1 σ n − s − Υ ∗ n − s +1 θ ∗ n (cid:21) (cid:20)(cid:0) − (cid:1) n − s +2 σ n − s Υ n − s +2 θ n (cid:21) = − σ n − s ) − Υ ∗ n − s +1 Υ n − s +2 | θ n | . (7.28)59he factor Υ ∗ n − s +1 Υ n − s +2 is then further expanded by f j as shown in (7.9),and we identify the imaginary part of a (cid:48) ,Υ ∗ n − s +1 Υ n − s +2 = (cid:104) f n − s +1 a (cid:48)∗ − f n − s +2 σ (cid:105)(cid:104) f n − s +2 a (cid:48) − f n − s +3 σ (cid:105) = − (cid:16) f n − s +2 a (cid:48) + f n − s +1 f n − s +3 a (cid:48)∗ (cid:17) σ + · · · = − i c (cid:16) f n − s +1 f n − s +3 − f n − s +2 (cid:17) σ + · · · = i c σ s − + · · · , (7.29)where we have used (7.19) and the fact a (cid:48) = a − i c . This implies (cid:101) D n − s, ∗ (cid:101) D n − s +1 , = − i c σ n − | θ n | + · · · = − i cσ (cid:12)(cid:12) (cid:101) D n (cid:12)(cid:12) + · · · , (7.30)where . . . will have vanishing contributions to the integral (7.5). Eqs. (7.27) and(7.30) are the sought-after relations between (cid:101) D s, (cid:101) D s +1 , ∗ or (cid:101) D n − s, ∗ (cid:101) D n − s +1 , and (cid:12)(cid:12) (cid:101) D n (cid:12)(cid:12) .Next we turn to the scaling behavior of (cid:12)(cid:12) (cid:101) D n (cid:12)(cid:12) with n . From (7.21) andfollowing the procedures that lead to the general expression of Υ j , we find | θ n | = f c + 2 f c + f + 2 σ n − c > , (7.31)where we have used (7.19) for the case k = nf − f f = σ n − , (7.32)to simplify | θ n | .To draw further information about | θ n | , we would like to discuss the genericbehavior of f k with respect to ω . We first note that when a − σ < µ , µ of the characteristic equation (7.13) become complex-conjugated. If we write them in terms of polar coordinate, then we have µ = σ e i ψ , µ = σ e − i ψ , ψ = tan − √ σ − a a , (7.33)with 0 ≤ ψ ≤ π . Recall that a = − ω + Ω and c = 2 ωγ . In terms of thefrequency the condition a − σ < (cid:112) Ω − σ < ω < (cid:112) Ω + 2 σ , (7.34)60ithin which the parameter a monotonically decreases from +2 σ to − σ and ψ steadily grows from 0 to π as ω increases. Hence f k can be written as f k = σ n − k sin( n − k + 1) ψ sin ψ , (7.35)which is heavily oscillating in ω . If we substitute (7.35) into (7.31), we find that | θ n | becomes | θ n | = σ n sin ( n + 1) ψ sin ψ + 2 σ n − c sin nψ sin ψ + σ n − sin ( n − ψ sin ψ + 2 σ n − c = σ n sin ψ (cid:20) sin ( n + 1) ψ + 2 c σ (cid:16) nψ (cid:17) + c σ sin ( n − ψ (cid:21) . (7.36)It can be greatly simplified when the coupling with the environment is weaksuch that γ Ω (cid:28) σ . In this case | θ n | reduces to | θ n | ∼ σ n sin ψ (cid:104) sin ( n + 1) ψ + ε (cid:105) . (7.37)Here ε is a small positive number as a reminder that the expression in thesquared brackets is not supposed to totally vanish so that when we put | θ n | back to (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) , it will not introduce artefact poles, (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) = 1 σ sin ψ sin ( n + 1) ψ + ε . (7.38)Owing to the factor sin ( n +1) ψ in the denominator, (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) will have maximaat ψ = kπ/ ( n +1) for k = 1 , , . . . , n . As for k = 0, or n +1, since the numeratorsin ψ cancels with the denominator sin ( n + 1) ψ , there is no maximum of (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) at these two locations. It implies that there are n peaks distributedevenly within the band √ Ω − σ < ω < √ Ω + 2 σ . As n increases, the peaksbecome narrower with the width of the order π/ ( n + 1). in particular when σ (cid:28) Ω because in that limit, ω = (cid:112) Ω − σ cos ψ ∼ Ω − σ Ω cos ψ .
For the neighboring maxima, the separation between them is given by − σ Ω (cid:104) cos( k +1)∆ − cos k ∆ (cid:105) = − σ Ω (cid:104) cos k ∆ cos ∆ − sin k ∆ sin ∆ − cos k ∆ (cid:105) (cid:39) (cid:16) σ Ω ∆ (cid:17) sin k ∆+ O (∆ ) , where ∆ = π/n (cid:28)
1. Thus the separation is independent of k in the neighborhood k ∆ (cid:28) π − k ∆ (cid:28)
61n the other hand, outside the band √ Ω − σ < ω < √ Ω + 2 σ , becausewe have a − σ >
0, both roots of the characteristic equations are real with µ > µ > | µ | > | µ | , a > σ , (7.39) | µ | > | µ | , a < − σ . (7.40)Hence we note that when a > σ , µ k rapidly dominates over µ k as k increases,but when a < − σ , µ k quickly outgrows µ k for large enough k . Thus we have f k ∼ µ n − k +11 µ − µ , a > +2 σ ,µ n − k +12 µ − µ , a < − σ , (7.41)for n − k (cid:29)
1. If we further assume γ Ω (cid:28) σ , then | θ n | will be approximatelygiven by | θ n | ∼ µ n +2 i a − σ (cid:40) c µ i (cid:34) (cid:18) σµ i (cid:19) n (cid:35) + · · · (cid:41) , (7.42)and then (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) (cid:39) a − σ µ i (cid:18) σµ i (cid:19) n − (cid:40) c µ i (cid:34) (cid:18) σµ i (cid:19) n (cid:35) + · · · (cid:41) − , (7.43)with i = 1 for a > σ but i = 2 for a < − σ . Now we note that although thestrong inter-oscillator coupling is allowed, the coupling strength σ is requiredsmaller than Ω / a = − ω + Ω , we have a − σ > ω < √ Ω − σ < <ω < √ Ω + 2 σ . It implies µ − σ = a + √ a − σ − σ > , < σµ < , a > σ , (7.44) − µ − σ = − a + √ a − σ − σ > , − < σµ < , a < − σ . (7.45)We then may conclude that the factor (cid:18) σµ i (cid:19) n − (cid:28) , (7.46)62rops to zero very fast for sufficiently large n , if 0 < ω < √ Ω − σ or ω > √ Ω + 2 σ . Thus, (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) will monotonically and rapidly falls to a relativelysmall value outside the band √ Ω − σ < ω < √ Ω + 2 σ .At this point we may draw some conclusions about (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) . In the fre-quency space, (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) falls monotonically and rapidly to a relatively smallvalue outside the band √ Ω − σ < ω < √ Ω + 2 σ ; on the other hand, withinthe frequency band, it possesses comb-like structure. The number of spikesgrows with the length of the chain, but the width of the spike, on the con-trary, becomes narrower and narrower, inversely proportional to the length ofthe chain.Next how the behavior of the transmission coefficient helps to understandthe dependence of the steady current on the length of the chain? A simplerquestion to ask is how the contribution from each spike in (7.38) will scale with n within the frequency band? We first make an observation for the integral I ( k ) n = (cid:90) k +1 / n +1 π k − / n +1 π dψ sin ψ sin ( n + 1) ψ + (cid:15) (7.47)for the k th spike among n spikes confined within the interval 0 < ψ < π .The parameter (cid:15) is a very small positive number. Change of the variable ψ to (cid:36) = ( n + 1) ψ gives I ( k ) n = 1 n + 1 (cid:90) ( k + ) π ( k − ) π d(cid:36) sin (cid:36)n +1 sin (cid:36) + (cid:15) . (7.48)For large n , the numerator of the integrand is slowly varying compared with thedenominator, so it can be pulled out of the integral and is evaluated for (cid:36) = kπ , I ( k ) n (cid:39) n + 1 sin k πn + 1 (cid:90) ( k + ) π ( k − ) π d(cid:36) (cid:36) + (cid:15) = 1 n + 1 (cid:34) sin k πn + 1 (cid:90) π − π d(cid:36) (cid:36) + (cid:15) (cid:35) . (7.49)The approximation improves for a larger value of n because the denominator of(7.48) becomes more slowly varying with ω . The contribution from the squaredbrackets is approximately the same for the spike at about the same locations63ithin the interval 0 < ψ < π . Thus I ( k ) n will scale with n − for sufficientlylarge n . For example, let us pick one spike, say, at k = n / n = n . Nowsuppose we rescale n from n = n to n = 3 n , and then we see that the threespikes centered at about k = 3 / n . Therefore each spike in the n = 3 n case will contribute one third as much as that in the n = n case tothe integral in I ( k ) n , as we can see from (7.49). When we consider all the spikeswith the range 0 < ψ < π , we expect I n = n (cid:88) k =1 I ( k ) n (7.50)should independent of n for sufficiently large n .Following this argument, we see the steady current J becomes J = 8 γ (cid:90) ∞−∞ dω ω (cid:12)(cid:12) (cid:101) D n ( ω ) (cid:12)(cid:12) (cid:104) (cid:101) G H ( ω ) − (cid:101) G nnH ( ω ) (cid:105) (cid:39) γ σ (cid:90) √ Ω +2 σ √ Ω − σ dω sin ψ sin ( n + 1) ψ + ε (cid:110) ω (cid:104) (cid:101) G H ( ω ) − (cid:101) G nnH ( ω ) (cid:105)(cid:111) = 16 γ σ (cid:90) π dψ sin ψ sin ( n + 1) ψ + ε h ( ω ) , (7.51)with ω = (cid:112) Ω − σ cos ψ , and the function h ( ω ) being given by h ( ω ) = ω (cid:104) (cid:101) G H ( ω ) − (cid:101) G nnH ( ω ) (cid:105) . (7.52)Since only the denominator sin ( n + 1) ψ + ε is vert rapidly oscillating with ψ forlarge n , we can use the previous arguments to support that the steady currentdoes not scale with n for sufficiently large n , that is J (cid:39) O ( n ) , (7.53)when n > N for some large positive integer N .
8. ReferencesReferences [1] See, e.g., H. Spohn and J. L. Lebowitz, “Stationary non-equilibrium statesof infinite harmonic systems”, Commun. Math. Phys. , 97 (1977); G.64allavotti and E. D. G. Cohen , “Dynamical ensembles in stationarystates”, J. Stat. Phys. , 931 (1995); L. Rey-Bellet and L. E. Thomas,“Exponential convergence to non-equilibrium stationary states in classi-cal statistical mechanics”, Commun. Math. Phys. , 305 (2002); J.-P. Eckmann, “Non-equilibrium steady states” Beijing Lecture (2002),arXiv:math-ph/0304043; B. Derrida, “Non-equilibrium steady states:fluctuations and large deviations of the density and of the current”, J.Stat. Mech. P07023 (2007) and references therein.[2] See, e.g., S. R. de Groot and P. Mazur, Non-equilibrium thermodynam-ics (North-Holland, 1962; Dover, 1984); J. Vollmer, Phys. Rep. , 131(2002) and references therein; G. Gallavotti, “Entropy production andthermodynamics of nonequilibrium stationary states: a point of view”,Chaos , 680 (2004); S. Sasa and H. Tasaki, “Steady state thermody-namics”, J. Stat. Phys. , 125 (2006) and references therein.[3] R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford UniversityPress, Oxford, 2001).[4] C. W. Gardiner and P. Zoller,
Quantum Noise, 2nd Enlarge Edition (Springer, Berlin, 2000).[5] U. Weiss,
Quantum Dissipative Systems (World Scientific, Singapore,1993).[6] E. Joos, H. D. Zeh, C. Kiefer, D. J. W. Guilini, J. Kupsch and I.-I. Sta-matescu,