Non-Equilibrium Thermodynamics

Chapter 2 Non-equilibrium Thermodynamics 2.1 Response, Relaxation and Correlation At the beginning of the 21st century, the thermodynamics of systems far from equilibrium remains poorly understood. However, it turns out that many nonequilibrium phenomena can be described rather well in terms of equilibrium fluctuations; this is especially true of systems near equilibrium [1, 2]. By designating a system as “near equilibrium”, we mean that the system is perturbed from its equilibrium state by some time-dependent external force f (t). The external force is deterministic, not random; typical examples include mechanical forces and forces due to an applied electric or magnetic field. This force drives the expectation values of some of the system’s observables away from their equilibrium values. For example, a typical observable A affected by the external force might be the system’s velocity or its magnetic moment. If the response of the observable A to the external force f (t) satisfies the linearity property δA(λf (t), t) = λδA(f (t), t) (2.1) where δA = A − ⟨A⟩eq and λ describes the strength of the force, then we call the time-dependent behavior of A the linear response of A to the external force f (t). The linearity property Eq.(2.1) implies that the shape of the response curve A vs. t is independent of the value of λ in the case of linear response. After achieving a short-lived nonequilibrium steady state (between t2 and t3 in Figure 2.1), the system is allowed to relax back to equilibrium. This process is also known as regression. Linear response and regression of a system driven from equilibrium are both described in terms of the time correlation function of the observable A, and so we turn first to the definition and properties of the time correlation function [3, 4]. The time correlation function CAA (t, t′ ) of the observable A is defined by CAA (t, t′ ) = ⟨A(t)A(t′ )⟩ = Tr[A(t)A(t′ )ρeq ] Tr[ρeq ] (2.2) Here, ρeq denotes the equilibrium density matrix of the system; hence the average denoted by ⟨⟩ is the ensemble average. This function describes how the value of A at time t is correlated to its 2 Figure 2.1: Response, nonequilibrium steady state, and relaxation value at time t′ ; it is sometimes referred to as the autocorrelation function of A to distinguish it from correlation functions between A and other observables. For a system which is time-translational invariant, we often choose for convenience to set t′ = 0 and to drop the subscript on CAA , so that the time correlation function becomes simply C(t) = ⟨A(t)A(0)⟩ (2.3) The correlation function may in general take on complex values. This result is in keeping with our phenomenological understanding of quantum mechanics in the following way. In order to mea­ sure the correlation function of an observable A, the quantity A must be measured twice (first at time zero, then again at time t). However, the first measurement at t = 0 collapses the system wavepacket, and the state that would have been exhibited by the unperturbed system at time t becomes irrecoverable. We now identify some important features and properties of correlation functions. 1. All inner products ⟨X|Y ⟩ satisfy the Schwarz inequality |⟨X|Y ⟩|2 ≤ ⟨X 2 ⟩⟨Y 2 ⟩ (2.4) Thus the correlation function for any relaxation process has the property C 2 (t) = |⟨A(t)|A(0)⟩|2 ≤ ⟨A2 (t)A2 (0)⟩ ≤ ⟨A2 (0)⟩2 = C 2 (0) (2.5) The second inequality above arises from the fact that A2 (t) < A2 (0) for relaxation processes when t > 0. More concisely, the Schwarz inequality implies that |C(t)| ≤ C(0) (2.6) 2. Correlation functions are time-invariant, that is, their value depends only on the time interval between the two measurements of the observable: ⟨A(t)A(0)⟩ = ⟨A(t − t0 )A(t0 )⟩ = ⟨A(0)A(−t)⟩ (2.7) 3 3. Time-invariance imparts the following identities on the time derivative of a time correlation function: Ċ(t) = ⟨Ȧ(t)A(0)⟩ = −⟨A(0)Ȧ(−t)⟩ = −⟨A(t)Ȧ(0)⟩ (2.8) 4. If the equilibrium value of A is ⟨A⟩eq = 0, then the long-time limit of the correlation function is zero, (2.9) lim ⟨A(t)A(0)⟩ = ⟨A⟩eq ⟨A(0)⟩ = 0 t→∞ 5. For quantum systems, the time-invariance properties imply that C(−t) = C ∗ (t). In the classical limit, the correlation function is always real-valued, so this relation becomes C(−t) = C(t) and C(t) is thus even. The fact that classical correlation functions are real-valued should seem sensible because we can (and do) measure correlation functions every day for classical systems, for example, when we try to steady a cord dangling from the ceiling. In this case, we determine the appropriate time and place to apply an external steadying force by looking ˙ for time correlations between the various motions the cord undergoes. Note that C(t) is odd ˙ with C(0) = 0 for classical time correlation functions. 6. For ergodic systems, the time correlation function can be calculated as a time average instead of an ensemble average: ∫ 1 T ⟨A(t)A(0)⟩ = lim A(t + τ ′ )A(τ ′ ) dτ ′ (2.10) T →∞ τ 0 Since most systems amenable to analysis by the methods of statistical mechanics are inher­ ently ergodic, we are generally free to choose whichever formulation is easier to work with. The time average is often easier to implement experimentally because it only requires inte­ gration along a trajectory rather than a simultaneous sampling of every state accessible to the system. Example: The classical linear harmonic oscillator with mass m and frequency ω obeys the equation of motion ẍ + ω 2 x = 0 If we provide initial conditions x(0) and ẋ(0) = v(0), then this equation of motion has the closed-form solution v(0) x(t) = x(0) cos ωt + sin ωt ω Taking the inner product of x(t) with the initial value x(0), we find ⟨x(t)x(0)⟩ = ⟨x2 (0)⟩ cos ωt + ⟨v(0)x(0)⟩ sin ωt ω The second term is zero because ⟨x(0)⟩ = 0, so the time correlation function is just C(t) = ⟨x2 (0)⟩ cos ωt kT Finally, invoking the equipartition result ⟨x2 (0)⟩ = m� 2 , where k is the Boltzmann constant, the correlation function for the classical linear harmonic oscillator is C(t) = kT cos ωt mω 2 4 2.2 Onsager Regression Theory At first glance, the relaxation of macroscopic non-equilibrium disturbances in a system might seem completely unrelated to the regression of microscopic fluctuations in the corresponding equilibrium system. However, they are intimately related by so-called fluctuation-dissipation theorems. The existence of this link between microscopic fluctuations and macroscopic relaxation was conjectured by Lars Onsager in 1931, some twenty years before it was finally proven to be true; hence it is often referred to as the Onsager regression hypothesis. To formulate the hypothesis, we consider an observable A with ⟨A⟩eq = 0 that takes on a nonequilibrium average value ∆A due to an applied external force f which acts during the time interval t ≤ 0 but becomes identically zero for t > 0. For t ≤ 0, the ensemble average of ∆A can be expressed as ∆A = [ ] ⟨Ae−)(H−f A) ⟩ 2 � βf ⟨A(0)A(0)⟩ − ⟨A(0)⟩ = βf C(0) ⟨e−)(H−f A) ⟩ (2.11) where the approximation being made is truncation of the Taylor series for each exponential to first order. For t > 0, the system evolves according to H instead of H − f A, so ∆A is no longer stationary, but acquires a time-dependence: ∆A = [ ] ⟨A(t)e−)(H−f A) ⟩ � βf ⟨A(t)A(0)⟩ − ⟨A(0)⟩2 = βf C(t) −)(H−f A) ⟨e ⟩ (2.12) Onsager’s hypothesis states that the relaxation of the non-equilibrium value of ∆A is related to its value at t = 0 in the same way that the time correlation function for a spontaneous fluctuation is related to its value at t = 0: ∆A(t) C(t) = (2.13) ∆A(0) C(0) Example: The transition state theory of chemical kinetics can be formulated through the Onsager relation we’ve just presented. Consider a chemical equilibrium established between two species A and B, kf A�B kb with forward rate constant kf and backward rate constant kb . We can describe the population dynamics of A and B deterministically in the macroscopic limit through a pair of coupled differential equations, { ṖA = −kf PA + kb PB 1. Equilibrium popuations ṖB = k f P A − kb P B The equilibrium state of this system satisfies the detailed balance condition kf ⟨PA ⟩ = kb ⟨PB ⟩ 5 where the angle brackets denote the equilibrium values of the populations. Taking the populations to be normalized to unity, ⟨PA ⟩ + ⟨PB ⟩ = 1, we can express ⟨PA ⟩ in terms of the rate constants: ⟨PA ⟩ = ⟨PA ⟩ kb = kf + kb ⟨PA ⟩ + ⟨PB ⟩ For notational simplicity, we introduce k = kf + kb and refer to the equilibrium populations ⟨PA ⟩ and ⟨PB ⟩ by qA and qB , respectively. With this new notation, we can express the equilibrium populations of A and B as  kb   qA = k k  q = f B k If the initial state is all species A, the solution to the coupled differential equations indicates a decay to equilibrium with rate constant k, which we can write in terms of ∆PA (t) = PA (t) − qA as ∆PA (t) = ∆PA (0)e−kt Setting this result aside for a moment, note that if we consider the energies of species A and B to be potential wells connected along a reaction coordinate x, then we can write down an expression for the fluctuation in occupation number n for each species as a function of x. The barrier between the A and B potential wells is a maximum at x = xb ; see Figure 2.2. Figure 2.2: Projection of the potential energy surface connecting species A and B along reaction coordinate x 2. Application of Onsager Regression hypothesis To reflect the fact that a particle to the left of the barrier is species A and a particle to the right is species B, we write the occupation numbers in terms of the Heaviside step function, { nA = θ(xb − x) nB = θ(x − xb ) where ⟨nA ⟩ = qA and ⟨nB ⟩ = qB . Applying Onsager’s regression hypothesis to this example, we can relate the dissipation of PA to the fluctuations in occupation number as follows: C(t) ⟨δnA (t)δnA (0)⟩ ∆PA (t) = = = e−kt 2 ∆PA (0) C(0) ⟨δnA (0)⟩ 6 The second equality arises from our integrated rate equation for the dissipation of PA . Also note that 2 2 ⟨δnA ⟩ = ⟨n2A ⟩ − ⟨nA ⟩2 = qA − qA = qA − qA (1 − qB ) = qA qB Differentiating the fluctuation-dissipation relation above with respect to t and invoking the identity just shown, we find ⟨δṅA (t)δnA (0)⟩ ⟨nA (t)ṅA (0)⟩ ke−kt = − = qA qB ⟨δn2A (0)⟩ Recasting this equation in terms of the reaction coordinate x, we arrive at an expression for the time dependence of the forward rate constant kf (t), kf (t) = kf e−kt = ⟨θ (x (t) − xb ) δ (xb − x (0)) v⟩ ⟨θ (xb − x (t))⟩ where v = ṅA (0) is the initial rate of reaction. Finally, to determine the transition state theory (TST) rate constant, we consider our time-dependent expression for kf in the short-time limit, since transition states typically only survive a few molecular vibrations. In this limit, 3. Expression for the TST rate constant lim kf (t) = t→0+ ⟨θ(x(0+ ) − xb )δ(xb − x(0))v⟩ ⟨θ(v)δ(xb − x(0))v⟩ = ⟨nA ⟩ ⟨nA ⟩ From the kinetic theory of gases, we recognize that √ kB T ⟨θ(v)v⟩ = = (2πmβ)−1/2 2πm If we stipulate now that the height of the barrier is Eb , some rearrangement of the preceding formulas reveals that √ ⟨δ(xb − x)⟩ mω 2 β −(Eb = e 2π ⟨θ(xb − x)⟩ where ω is the fundamental frequency of the left potential well. It follows that the TST rate constant takes on the simple form ω −(Eb kT ST = e 2π To conclude our excursion into TST kinetics, note that the ratio k(t) ⟨θ(x(t) − xb )δ(x(0) − xb )v⟩ = kT ST ⟨θ(x(0+ ) − xb )δ(x(0) − xb )v⟩ is always less than or equal to one. This result indicates that the TST flux is partially trapped in the product well while part of the TST flux recrosses back to the reactant state. This result is in keeping with our intuition of chemical dynamics in that every macroscopic reaction is, to some degree, a process of establishing equilibrium rather than a perfect flow from all reactants to all products. 2.3 Linear Response Theory and Causality The concept of linear response was introduced in section 2.1. Here, we explore further how the linear response of a system is quantified by considering the important relations regularly invoked by practitioners of linear response theory. 7 2.3.1 Response Functions The motivating idea behind linear response is that the response of a system to an external force depends on the strength of that force at all times during which the force acts on the system. That is, the response at time t depends on the history of the force’s action on the system. An appropriately weighted sum of the strength of the external force at each moment during the interaction will describe the overall response. Mathematically, therefore, we express the response as an integral over the history of the interaction, ∫ ∞ ∆A(t) = K(t, τ )f (τ ) dτ (2.14) −∞ The kernel K(t, τ ) in this expression, which provides the weight for the strength of the external force at each time, is called the response function. The response function has two very important properties: • Time invariance: K depends only on the time interval between τ and t, not on the two times independently. More succinctly, K(t, τ ) = K(t − τ ) (2.15) • Causality: The system cannot respond until the force has been applied. This places an upper limit of t on the integration over the history of the external force. With these observations in place, we arrive at the standard formula describing the linear response of an observable A to an external force f (t), ∫ t ∆A(t) = K(t − τ )f (τ ) dτ (2.16) −∞ Linear response – described by the response function K(t) – and linear regression – described by the time correlation function C(t) – are directly related to one another. To see the connection, consider a force f (t) which is constant with strength f for t ≤ 0 and is zero for t > 0. We have established two ways to describe the response of an observable A to this force: • Linear regression: ∆A(t) = βf C(t) ∫0 • Linear response: ∆A(t) = −∞ K(t − τ )f (τ ) dτ From this information, we conclude that the correlation function and response function are related by K(t) = −βĊ(t)θ(t) (2.17) where θ(t) is the Heaviside function. Sometimes the linear response function is more conveniently expressed in the frequency domain, in which case it is called the frequency-dependent response function. In many physical situations, it plays the role of a susceptibility to a force and consequently is denoted by χ(ω), ∫ ∞ χ(ω) = ei t K(t) dt (2.18) 0 This response function is often partitioned into real and imaginary parts, which can also be thought of as even and odd parts, respectively, χ(ω) = χ′ (ω) + iχ′′ (ω) (2.19) 8 Example: The response function for the classical linear harmonic oscillator can be quickly deduced from its time correlation function. Recall from the first example in this chapter that the time correlation function for the classical linear harmonic oscillator is kB T C(t) = cos ωt mω 2 Applying Eq.(2.17), we differentiate with respect to t and multiply by β = kB1T to determine that 1 sin (ωt)θ(t) mω This is the response function for the classical linear harmonic oscillator. K(t) = 2.3.2 Absorption Power Spectra The frequency-dependent response function is directly related to the absorption spectrum: in fact, knowledge of χ(ω) and the time-dependent external force f (t) is sufficient to fully describe the absorption spectrum. The rate at which work is done on a system by a generalized external force f (t) is f (t)Ȧ(t), where A is the observable corresponding to the generalized force f . This quantity has units of power, so we can calculate the total absorption energy by integrating this power over time, ∫ ∫ P (t) dt = f (t)Ȧ(t) dt (2.20) To recast this result in terms of χ(ω), we first consider the Fourier transform of the time-dependent observable A, ∫ Ã(ω) = Applying Eq.(2.16), we have ∫ ˜ A(ω) = ei t A(t) dt (2.21) ∫ ei t K(t − τ )f (τ ) dτ dt (2.22) ˜ The following rearrangements allow us to express A(ω) entirely in terms of frequency-dependent functions: ∫ ∫ Ã(ω) = ei (t−T ) ei T K(t − τ )f (τ ) dτ dt (2.23) ∫ ∫ = ei (t−T ) K(t − τ ) d(t − τ ) ei T f (τ ) dt (2.24) = χ(ω)f˜(ω) (2.25) where f˜(ω) is the Fourier transform of f (t). Returning to Eq.(2.20) for the absorption energy, ∫ ∫ ∫ 1 P (t) dt = f (t)Ȧ(t) dt = f˜(−ω)A˜˙ (ω) dω (2.26) 2π Rearranging the expression once again and evaluating the integral over time yields ∫ ∫ 1 P (t) dt = (−iω)f˜(−ω)Ã(ω) dω (2.27) 2π ∫ 1 (2.28) = (−iω)χ(ω)|f˜(ω)|2 dω 2π ∫ 1 (2.29) = ωχ′′ (ω)|f˜(ω)|2 dω 2π 9 Hence the absorption power spectrum P (ω) has the form P (ω) = ωχ′′ (ω)|f˜(ω)|2 (2.30) Example: For a monochromatic force F (t) = F� cos ω0 t, the Fourier transform of F (t) is given by F̃ (ω) = F� π[δ(ω − ω0 ) + δ(ω + ω0 )] Hence our expression for the absorption power spectrum Eq.(2.30) tells us that the absorption rate (i.e. the time average of P (t)) for such a system is ∫ 1 T ω lim P (t) dt = χ′′ (ω)|F� |2 T →∞ τ 0 2 Furthermore, using Eq.(2.17) and the fact that χ′′ (ω) is odd, we find that ∫ ∞ χ′′ (ω) = sin ωtK(t) dt ∫0 ∞ ( ) = sin ωt −βĊ(t) dt 0 ∫ ∞ = βω C(t) cos ωt dt ∫0 ∞ βω = ei�t C(t) dt 2 −∞ βω ˜ = C(ω) 2 This simple relationship illustrates the close connection between frequency-dependent response functions (or susceptibilities) and time correlation functions. 2.3.3 Causality and the Kramers-Kronig Relations Our final consideration in this chapter is the relationship between the real and imaginary parts of the frequency-dependent response function as defined in Eq.(2.19). The equations relating these two functions are known as the Kramers-Kronig relations. In their most general form, they govern the response function as a function of the complex frequency z = ω +iϵ, though under most physical circumstances of interest they can be expressed in terms of real-valued frequencies alone. The relations arise from the causality requirement, which we originally expressed by∫ requiring ∞ K(t) = 0, ∀t < 0. It turns out that this requirement, along with the assumption that 0 K(t) dt converges, implies that the response function χ(z) is analytic on the upper half of the complex plane. Consider integrating the function χ(z) z − ω0 around the contour in the complex plane shown in Figure 2.3. The contour has been explicitly constructed to avoid the pole at z = ω0 , so Cauchy’s integral theorem implies that ∮ χ(z) dz = 0 (2.31) z − ω0 10 Figure 2.3: Contour in the complex frequency plane used to motivate the Kramers-Kronig relations We can also, however, integrate piecewise over each part of the contour; some manipulation with the residue theorem is required, but the final result is ∮ ∫ ∞ χ(z) χ(z) (2.32) dz = P dz − iπχ(ω0 ) z z − ω0 −∞ − ω0 where P denotes the Cauchy principal value of the integral. Setting the two results above equal and solving for the response function, ∫ ∞ 1 χ(z) χ(ω) = P dz (2.33) iπ −∞ z − ω The decomposition of Eq.(2.33) into real and imaginary parts yields the Kramers-Kronig relations, ∫ ∞ ′′ 1 χ (z) ′ dz (2.34) χ (ω) = P π −∞ z − ω ∫ ∞ ′ 1 χ (z) ′′ dz (2.35) χ (ω) = − P π −∞ z − ω As promised, this pair of equations provides a concise relationship between the real and imaginary parts of any response function χ(ω). References [1] N. Kubo, R.; Saito and N. Hashitsume. Statistical Physics II: Nonequilibrium Statistical Me­ chanics. Springer, 2003. [2] Robert Zwanzig. Nonequilibrium Statistical Mechanics. New York: Oxford University Press, 2001. [3] David Chandler. Introduction to Modern Statistical Mechanics. New York: Oxford University Press, 1987. [4] L.E. Reichl. A Modern Course in Statistical Physics. New York: Wiley-Interscience, 1998.