Classical and quantum stochastic thermodynamics

Abstract

The stochastic thermodynamics provides a framework for the description of systems that are out of thermodynamic equilibrium. It is based on the assumption that the elementary constituents are acted by random forces that generate a stochastic dynamics, which is here represented by a Fokker-Planck-Kramers equation. We emphasize the role of the irreversible probability current, the vanishing of which characterizes the thermodynamic equilibrium and yields a special relation between fluctuation and dissipation. The connection to thermodynamics is obtained by the definition of the energy function and the entropy as well as the rate at which entropy is generated. The extension to quantum systems is provided by a quantum evolution equation which is a canonical quantization of the Fokker-Planck-Kramers equation. An example of an irreversible systems is presented which shows a nonequilibrium stationary state with an unceasing production of entropy. A relationship between the fluxes and the path integral is also presented.

Keywords:
Stochastic thermodynamics; Quantum thermodynamics; Entropy production

Introduction

Thermodynamics was conceived as a discipline based on principles and laws that refer to macroscopic quantities such as the principles of energy conservation and of entropy increase, which are the first and second laws of thermodynamics. Although these two principles are valid for systems in equilibrium as well as for systems out of equilibrium, the initial development of the discipline lead to the establishing of a theory of thermodynamics of system in equilibrium. This was possible because the energy of a system in thermodynamic equilibrium is related functionally to the entropy, which allows the definition of temperature.

The derivation of the thermodynamics from the microscopic laws of motion was the aim of the kinetic theory. One of its consequences was the development of the equilibrium statistical mechanics advanced by Gibbs. The statistical mechanics is based on the description of system by the probability distribution that bears the name of Gibbs, which for a system in contact with a heat reservoir is proportional to $e^{-E/kT}$, where E is the energy function and T, the temperature. The crucial property of the equilibrium distribution is that the probability depends on the states of the system only through the energy function. This property, along with the Gibbs expression for the entropy, leads to the relation between energy and entropy, mentioned above, which characterizes the thermodynamic equilibrium.

The entropy of an isolated system in equilibrium remains invariant. But if it is not in equilibrium, its entropy increases and the increase, in this case, is not due to the flux of entropy because the system is isolated. Entropy is being created spontaneously and in this sense it differs from the energy which is a conserved quantity. If a system is not isolated then the variation of the entropy S with time is the algebraic sum of two terms,

(1) $$\frac{dS}{dt}=\mathrm{\Pi }-\mathrm{\Phi },$$

where $\mathrm{\Pi }$ is the rate in which entropy is being created, the rate of entropy production, and $\mathrm{\Phi }$ is the flux of entropy.

The production of entropy is related to irreversible processes occurring inside the system which are understood as processes that are more likely to occur than their time-reverse counterparts. Thermodynamic equilibrium is thus characterized as the state where a process and its time reversal are equally probable. This characterization of equilibrium, embodied in the stochastic thermodynamics, is a dynamical definition, being more comprehensible than the static definition given above in terms of the Gibbs distribution. The stochastic thermodynamics [1 –32] provides an approach to out of equilibrium an equilibrium thermodynamics that takes into account the dynamical characterization of the irreversible processes by assuming that a system evolves in time according to a microscopic stochastic dynamics.

The elementary constituents of a system are assumed to be acted by random forces in addition to the usual deterministic forces. As a consequence the trajectory followed by a particle is not in general deterministic. There are many possible trajectories that a particle may follow from a given point, each one with a certain probability of occurrence. The approach we follow here uses a representation of the dynamics in terms of the probability of the occurrence of a state at a certain instant of time, which is assumed to be governed by an evolution equation.

The main features of the approach that we follow here are: (1) a stochastic dynamics which is here represented by a Fokker-Planck-Kramers equation [33 –36]; (2) the assignment of an energy function; (3) the definition of entropy as having the same form as the Gibbs entropy; (4) a proper definition of the rate of entropy production.

The first part of this text is dedicated to the classical case. In the second part we extend the results obtained in the first part to the quantum case. In particular, we use as the evolution equation a quantum version of the Fokker-Planck-Kramers equation. In this case, the probability density is replaced by the probability density operator, usually called density matrix. In a third part we generalize the results for the case of many degrees of freedom and present an example of a system that displays a nonequilibrium stationary state with an unceasing production of entropy.

2. Evolution equation

Our object of study is a system of particles that interact among themselves and may also be subject to external forces. In addition each particle is acted by random forces, the origin of which may be internal or external to system. The system evolves in time according to the Newton equation of motion. Due to random forces, the trajectory is not uniquely defined. There are many possible trajectories starting from a given state, each one with a certain probability.

If the system is at a given state at the initial time, one may ask for the probability that it is at a given state at a later time. An answer to this question is provided by the Fokker-Planck-Kramers (FPK) equation which gives the time evolution of the probability density $\rho (x,p,t)$. We will focus initially on a system with just one degree of freedom in which case x is the position and p is the momentum of the particle, and both quantities constitute the state of the system. The probability that at time t the state of the system is inside $dxdp$ around $(x,p)$ is $\rho (x,p,t)dxdp$.

In contrast to the equilibrium statistical mechanics for which the probability density is constant in time, here it depends on time. If we wish that the system reaches thermal equilibrium for long times, usually called thermalization, then the solution of the evolution equation for long times must be a Gibbs equilibrium distribution, which is characterized by depending on $(x,p)$ only through the energy function associated to the system.

The FPK equation is given by

(2) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-\frac{p}{m}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }x}-\frac{\mathrm{\partial }f\rho }{\mathrm{\partial }p}+\frac{\mathrm{\Gamma }}{2}\frac{{\mathrm{\partial }}^2\rho }{\mathrm{\partial }{p}^2},$$

where m is the mass of the particle, f is the ordinary force acting on the particle and $\mathrm{\Gamma }$ is a constant that is associated to the stochastic forces. The force f is understood to be a sum of an internal force $F^c$, considered to be a conservative force, and a dissipative force D,

(3) $$f=F^c+D.$$

The property of the dissipative force D that distinguishes it from the other forces is that it is an odd function of the momentum.

A derivation of the FPK equation (2) can be obtained from a Langevin equation and can be found in reference [36]. The Langevin equation leading to (2) is the equation of motion for a particle of mass m moving along a straight line which in addition to the ordinary force f is also under the action of a stochastic force with zero mean and variance proportional to $\mathrm{\Gamma }$.

An essential property of the FPK equation, and for that matter of any equation that governs the time evolution of a probability distribution, is that it preserves the normalization of $\rho$, that is,

(4) $$\int \rho (x,p,t)dxdp=1,$$

for any instant of time, where the integration is performed on the whole space of states. If $\rho$ is normalized at the initial time, it remains normalized forever. The demonstration of this fundamental property of the FPK equation is given in the appendix 16.

As $F^c$ is a conservative force, we may write $F^c=-(\mathrm{\partial }H/\mathrm{\partial }x)$, $p/m=(\mathrm{\partial }H/\mathrm{\partial }p)$, where

(5) $$H=\frac{p^2}{2m}+V(x)$$

is the energy function, which is the sum of the kinetic energy and the potential energy V. The first term and the one involving $F^c$ of the FPK equation become

(6) $$-\frac{p}{m}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }x}-\frac{\mathrm{\partial }{F}^c\rho }{\mathrm{\partial }p}=-\frac{\mathrm{\partial }H}{\mathrm{\partial }p}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }x}+\frac{\mathrm{\partial }H}{\mathrm{\partial }x}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}.$$

The right-hand side of this equation is written in the abbreviated form as

(7) $$\frac{\mathrm{\partial }H}{\mathrm{\partial }x}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}-\frac{\mathrm{\partial }H}{\mathrm{\partial }p}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }x}=\{H,\rho \},$$

which is called the Poisson brackets. Replacing this result in the equation (2), the FPK equation acquires the form

(8) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=\{H,\rho \}-\frac{\mathrm{\partial }J}{\mathrm{\partial }p},$$

where

(9) $$J=D\rho -\frac{\mathrm{\Gamma }}{2}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}.$$

The FPK equation (2) can also be written in the form

(10) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-\frac{\mathrm{\partial }{J}_x}{\mathrm{\partial }x}-\frac{\mathrm{\partial }{J}_p}{\mathrm{\partial }p},$$

where $J_x=p/m$ and $J_p=F^c\rho +J$. In this form, the FKP equation is understood as a continuity equation and $J_x$ and $J_p$ are understood as the components of the probability current. This is the form that allowed us to show the property (4), as presented in the appendix 16. The part J of the component $J_p$ is the irreversible probability current, which plays a fundamental role in the present approach. We remark that without J the FPK equation reduces to the Liouville equation of classical statistical mechanics [37],

(11) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=\{H,\rho \}.$$

3. Thermodynamic equilibrium

The FPK equation as it stands may or may not describe a system that for long times will be in thermodynamic equilibrium. For long times the solution of the FPK equation is its stationary solution, that is, the solution obtained by setting to zero the right-hand side of the equation, in which case J may or may not vanish.

A system in thermodynamic equilibrium may be said to be the one described by a Gibbs distribution. However, this is a static definition. We need here a dynamic characterization of thermodynamic equilibrium. This is provided by characterizing the thermodynamic equilibrium as the stationary state such that the irreversible current vanishes. Therefore, the equilibrium distribution ${\rho }_e$ obeys the condition

(12) $$D{\rho }_e-\frac{\mathrm{\Gamma }}{2}\frac{\mathrm{\partial }{\rho }_e}{\mathrm{\partial }p}=0,$$

and in addition

(13) $$\{H,{\rho }_e\}=0.$$

The solution of this last condition leads to the result that ${\rho }_e$ is a function of H, that is, ${\rho }_e$ depends on x and p through the energy function $H(x,p)$. This characterizes any Gibbs equilibrium distribution.

The condition (12) is understood as the relation between dissipation, described by D, and noise or fluctuations, described by $\mathrm{\Gamma }$. That is, in equilibrium there must be a relation between dissipation and fluctuations. We wish to describe a system in contact with a heat reservoir at a temperature T, in which case the appropriate Gibbs equilibrium distribution is [37]

(14) $${\rho }_e=\frac{1}{Z}{e}^{-\beta H},$$

where $\beta =1/kT$ and k is the Boltzmann constant. Replacing (14) in equation (12), we find

(15) $$D=-\gamma p,$$

where $\gamma$ is the constant that connects $\mathrm{\Gamma }$ with temperature,

(16) $$\mathrm{\Gamma }=2m\gamma kT.$$

Notice that D is the usual type of dissipation force proportional to the velocity.

4. Energy, heat and entropy

A thermodynamic system may have its energy changing in time. The energy variation is due to the exchange of heat or work with the environment. Here we will treat the case where the external forces are absent so that the variation in energy is only due to the exchange of heat. The energy U of a thermodynamic system, sometimes called internal energy, is the average of the energy function introduced above, that is,

(17) $$U=\int H\rho dxdp,$$

and may depend on time as $\rho$ depends on time. The variation of U with time is

(18) $$\frac{dU}{dt}=\int H\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}dxdp.$$

Replacing the time derivative of $\rho$ by using the FPK equation in the form (8), we find

(19) $$\frac{dU}{dt}=-\int H\frac{\mathrm{\partial }J}{\mathrm{\partial }p}dxdp.$$

The term involving the Poisson brackets, vanishes after an integration by parts has been performed. Here and in the following, whenever we do an integration by parts, the integrated term is assumed to vanish. As shown in the appendix 16, this is so because we are considering that at the limits of integration the probability density vanishes rapidly. Performing the integral in (19) by parts, we find

(20) $$\frac{dU}{dt}=\int J\frac{\mathrm{\partial }H}{\mathrm{\partial }p}dxdp.$$

The right-hand side of this equation is interpreted as the rate at which heat is introduced into the system, or the flux of heat ${\mathrm{\Phi }}_q$,

(21) $${\mathrm{\Phi }}_q=\int J\frac{\mathrm{\partial }H}{\mathrm{\partial }p}dxdp.$$

Thus we may write

(22) $$\frac{dU}{dt}={\mathrm{\Phi }}_q,$$

an equation that may be understood as the conservation of energy. Notice that the flux of heat ${\mathrm{\Phi }}_q$ involves the irreversible part of the probability current.

From thermodynamics, we know that heat is related to entropy through the Clausius equation $dQ=TdS$, valid for systems in equilibrium, where $dQ$ is the infinitesimal heat exchanged with the system and $dS$ is the infinitesimal increase of the entropy S. Clausius equation is equivalent to the equation ${\mathrm{\Phi }}_q=T(dS/dt)$ which could be used to define entropy. However, this equation is of no use here because it is valid only for system in equilibrium. The appropriate way to define entropy of a system is to use the Gibbs form

(23) $$S=-k\int \rho \ln \rho dxdp,$$

which is a generalization of the Boltzmann entropy $S=k\ln W$ where W is the number of accessible states. Although, S given by (23) is usually used for system in equilibrium, here we are assuming that this form is also appropriate for systems out of equilibrium.

5. Entropy production

If the probability $\rho$ is found as a function of time by solving the FPK equation, then S is obtained as a function of time. Deriving (23) with respect to time, we find

(24) $$\frac{dS}{dt}=-k\int \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}\ln \rho dxdp.$$

There is another part involving the time derivative of $\ln \rho$ but it vanishes if we take into account that $\rho$ is normalized at any time, a result given by equation (4).

Replacing the time derivative of $\rho$ in (24), by using the FPK equation in the form (8), we find

(25) $$\frac{dS}{dt}=k\int \frac{\mathrm{\partial }J}{\mathrm{\partial }p}\ln \rho dxdp.$$

Again the part involving the Poisson brackets vanishes by an integration by parts. The entropy is not a conserved quantity like the energy and as a consequence its variation with time is not equal to the flux of entropy. In other terms, the right-hand side of (14) cannot be identified as the flux of entropy. In addition to the flux of entropy, there is another contribution related to the creation of entropy. This contribution is the rate of how entropy is being generated or created, which is called the rate of entropy production, denoted by $\mathrm{\Pi }$. Thus we variation of the entropy of a system with time is written as

(26) $$\frac{dS}{dt}=\mathrm{\Pi }-\mathrm{\Phi },$$

where $\mathrm{\Phi }$ is the flux of entropy from the system to the outside. The next step is to define or postulate the expression for one of the two quantities, $\mathrm{\Pi }$ or $\mathrm{\Phi }$. Once one of them is given, the other is obtained by observing that their difference should be equal to the right-hand side of (25).

The rate of entropy production $\mathrm{\Pi }$ is a quantity that vanishes when the thermodynamic equilibrium sets in and gives a measure of the deviation from equilibrium. As we have seen above, the vanishing of the irreversible probability current J is a condition for the the thermodynamic equilibrium. Since the entropy production is a nonnegative quantity and vanishes when J vanishes, it should be related to $J^2$. The expression for the $\mathrm{\Pi }$ that we are about to introduce meets this two conditions.

Let ${\rho }_0$ be the probability distribution that makes the irreversible probability current J to vanish. Writing J in the form

(27) $$\frac{J}{\rho }=D-\frac{\mathrm{\Gamma }}{2}\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }p},$$

this condition is equivalent to

(28) $$D=\frac{\mathrm{\Gamma }}{2}\frac{\mathrm{\partial }\ln {\rho }_0}{\mathrm{\partial }p}.$$

The probability ${\rho }_0$ does not need to be necessarily the equilibrium probability distribution ${\rho }_e$ because we are not requiring that $\{H,{\rho }_0\}$ vanishes as it occurs with ${\rho }_e$. In analogy with the right-hand side of equation (25), we define the rate of entropy production by

(29) $$\mathrm{\Pi }=k\int \frac{\mathrm{\partial }J}{\mathrm{\partial }p}(\ln \rho -\ln {\rho }_0)dxdp.$$

If we integrate by parts and use the relation

(30) $$\frac{J}{\rho }=\frac{\mathrm{\Gamma }}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }p}(\ln {\rho }_0-\ln \rho ),$$

that follows from (27) and (28), we reach the expression

(31) $$\mathrm{\Pi }=\frac{2k}{\mathrm{\Gamma }}\int \frac{J^2}{\rho }dxdp.$$

We see that the rate of entropy production is the integral of an expression proportional to $J^2$, as desired. It is nonnegative and vanishes in equilibrium, when $J=0$, that is,

(32) $$\mathrm{\Pi }\ge 0,$$

which is a brief statement of the second law of thermodynamics.

The flux of entropy $\mathrm{\Phi }$ is obtained from $dS/dt=\mathrm{\Pi }-\mathrm{\Phi }$ by using the expressions (25) and (29),

(33) $$\mathrm{\Phi }=-k\int \frac{\mathrm{\partial }J}{\mathrm{\partial }p}\ln {\rho }_{0}dxdp.$$

Performing an integration by parts and using (28), we find

(34) $$\mathrm{\Phi }=\frac{2k}{\mathrm{\Gamma }}\int JDdxdp.$$

The flux of entropy can also be written as

(35) $$\mathrm{\Phi }=\frac{2k}{\mathrm{\Gamma }}⟨D^2⟩+k⟨\frac{\mathrm{\partial }D}{\mathrm{\partial }p}⟩,$$

after the replacement of J, given by (8) in (34) and performing an integration by parts in the second term. This is an interesting form for the flux of entropy because it can be understood as an average over the probability distribution $\rho$, which is not the case of the rate of entropy production.

From the expressions for the flux of entropy $\mathrm{\Phi }$ and the rate of entropy production $\mathrm{\Pi }$, we draw an important conclusion concerning the Liouville equation (11). Since this equation can be understood a the FPK equation without the irreversible probability current J, and since $\mathrm{\Phi }$ and $\mathrm{\Pi }$ vanishes if $J=0$, it follows that the Liouville equation predicts no entropy production nor flux of entropy, and the entropy S is constant in time. If the Liouville equation is employed to describe a closed system that approaches equilibrium, but initially is out or equilibrium, then this prediction of the Liouville equation is in contradiction with thermodynamics which predicts an increase of entropy with time.

Let us consider in the following that D and $\mathrm{\Gamma }$ are related by (12) where ${\rho }_e$ is the canonical Gibbs distribution (14) so that the FPK equation describes a system in contact with a thermal reservoir and at equilibrium. Replacing the results (15) and (16) in (9), the expression for the irreversible probability current becomes

(36) $$J=-\gamma (p\rho +mkT\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}).$$

Replacing the expression for D given by (12) in the equation (34), we find

(37) $$\mathrm{\Phi }=k\int J\frac{\mathrm{\partial }\ln {\rho }_e}{\mathrm{\partial }p}dxdp=-\frac{1}{T}\int J\frac{\mathrm{\partial }H}{\mathrm{\partial }p}dxdp.$$

Comparing with (21), we see that the entropy flux and the heat flux are related by

(38) $$\mathrm{\Phi }=-\frac{1}{T}{\mathrm{\Phi }}_q.$$

Using equations (22) and (26) we get the following relation

(39) $$\frac{dU}{dt}-T\frac{dS}{dt}=-T\mathrm{\Pi },$$

valid at any time. Near equilibrium, $\mathrm{\Pi }$ vanishes faster then the other two time derivatives and

(40) $$\frac{dU}{dt}=T\frac{dS}{dt},$$

or $dU=TdS$ which is the Clausius relation, valid at thermodynamic equilibrium. We remark that T here is the temperature of the heat reservoir. The temperature of the system is $(\mathrm{\partial }U/\mathrm{\partial }S)=T^{\ast }$ if U could be written as a function of S. In out of equilibrium, when $\mathrm{\Pi }\ne 0$, this is not possible, but in equilibrium, in view of the relation $dU=TdS$, then U becomes a function of S. The relation $dU=TdS$ is translated into $T=(\mathrm{\partial }U/\mathrm{\partial }S)$, which implies $T=T^{\ast }$, and T becomes also the temperature of the system.

It is worth mentioning that the variation in time of the free energy F, defined by $F=U-TS$, at T constant, is

(41) $$\frac{dF}{dt}=-T\mathrm{\Pi }.$$

which follows form (39). Since $\mathrm{\Pi }\ge 0$, then $dF/dt\le 0$ and F decreases monotonically towards its equilibrium value. This inequality can also be viewed as the H theorem of Boltzmann. Defining the Boltzmann H by

(42) $$H=\int \rho \ln \frac{\rho }{{\rho }_e}dxdp,$$

we see that it is equal to $-\beta F$ plus a constant. Then, it follows from $dF/dt\le 0$ that $dH/dt\ge 0$, which is the H theorem of Boltzmann.

6. Work

The specific systems that we have consider so far are those that exchange only heat with the environment and are described by the FPK equation (8). Now we wish to consider the case where the systems are also subject to external forces. The appropriate way to treat this case is to add to the ordinary force appearing in the FPK equation (2) an external force $F^e$, so that f now reads

(43) $$f=F^c+F^e+D.$$

Repeating the reasoning leading to (8) from (2), we reach the evolution equation

(44) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=\{H,\rho \}-\frac{\mathrm{\partial }{F}^e\rho }{\mathrm{\partial }p}-\frac{\mathrm{\partial }J}{\mathrm{\partial }p},$$

where J is the irreversible probability current, given by (36), which is the one appropriate for the contact with a heat reservoir at a temperature T.

Due to the presence of the external force, the variation in time of the energy has another contribution in addition to the flux of heat,

(45) $$\frac{dU}{dt}={\mathrm{\Phi }}_q-{\mathrm{\Phi }}_w,$$

where ${\mathrm{\Phi }}_q$ is the heat flux into the system and ${\mathrm{\Phi }}_w$ is the work performed by the system per unit time against the the external forces, or power.

To determine the variation of energy with time, we proceed in the same way as we did to derive (22) from the evolution equation (2), but now we use the evolution equation (44). The result is the equation (45) where ${\mathrm{\Phi }}_q$ is the expression (21) and the power ${\mathrm{\Phi }}_w$ is

(46) $${\mathrm{\Phi }}_w=-\int F^e\rho \frac{\mathrm{\partial }H}{\mathrm{\partial }p}dxdp.$$

The equation (25) for the variation of entropy with time remains unchanged by the addition of the external force. To see this we replace the expression (44) into (24). The term involving the Poisson brackets vanishes as we have already seen. The term involving the external force is

(47) $$k\int \frac{\mathrm{\partial }{F}^e\rho }{\mathrm{\partial }p}\ln \rho dxdp=-k\int F^e\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}dxdp,$$

where we have performed an integration by parts. But this integral also vanishes if we assume that $F^e$ does not depend on p.

The rate of entropy production is defined by (31), and considering that the expression for $dS/dT$ remains unchanged, as we have just seen, so does the expression (35) for the entropy flux. As these relations are not modified by the presence of the external force, then the relation $\mathrm{\Phi }=-{\mathrm{\Phi }}_q/T$, as expressed by equation (38), between the entropy flux and the heat flux, valid for a system is in contact we a heat reservoir at a temperature T, also remains unchanged.

Taking into account that $dS/dt=\mathrm{\Pi }-\mathrm{\Phi }$ and that $dU/dt={\mathrm{\Phi }}_q-{\mathrm{\Phi }}_w$, we reach the following relation

(48) $$\frac{dU}{dt}-T\frac{dS}{dt}=-{\mathrm{\Phi }}_w-T\mathrm{\Pi }.$$

Considering a process in which T is constant, then the left hand side is $dF/dt$ where $F=U-TS$ is the free energy, that is,

(49) $$\frac{dF}{dt}=-{\mathrm{\Phi }}_w-T\mathrm{\Pi }.$$

Integrating in time, between $t_1$ and $t_2$, we find

(50) $$\mathrm{\Delta }F=-W-T{\int }_{t_1}^{t_2}\mathrm{\Pi }dt.$$

where W is the work performed by the external force,

(51) $$W={\int }_{t_1}^{t_2}{\mathrm{\Phi }}_{w}dt,$$

which is the time integration of the power ${\mathrm{\Phi }}_w$. Since $\mathrm{\Pi }\ge 0$, it follows that

(52) $$\mathrm{\Delta }F\le -W,$$

that is, the variation of the free energy is smaller than the work done on the system. The equality holds in an equilibrium process, when $\mathrm{\Pi }=0$.

The following remark is in order. When the system is in contact with a heat reservoir at a certain temperature T, it does not mean that T is the temperature of the system, as we have pointed out in the remark just below equation (40). If the rate of entropy production is nonzero, no temperature could be assigned to the system and $F=U-TS$ could not be, strictly speaking, the free energy of the system although U and S are the energy and entropy of the system. But we may suppose that for $t\le t_1$ and $t\ge t_2$, the system is in equilibrium, in which case $\mathrm{\Pi }$ is nonzero only for $t_1<t<t_2$. Within this scenario, T can be considered the temperature of the system at time $t_1$ and at time $t_2$, and F at these two instants of time will be the free energy of the system, and the relation $\mathrm{\Delta }F$ in (50) will represent the difference in the free energies of the system.

The derivations that we have just carried out, such as that of the inequality (52), made no restriction on the type of external force. It could be a nonconserved force or a time dependent force. This latter type of external force happens, for instance, when the system is driven at our will.

7. Harmonic oscillator

Let us apply the results we have found so far to the case of a harmonic oscillator for which $F^c=-Kx$ and

(53) $$H=\frac{p^2}{2m}+\frac{1}{2}K{x}^2.$$

It is in contact with a heat reservoir at a temperature T so that the FPK equation is

(54) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-\frac{p}{m}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }x}+Kx\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}+\gamma \frac{\mathrm{\partial }p\rho }{\mathrm{\partial }p}+\frac{m\gamma }{\beta }\frac{{\mathrm{\partial }}^2\rho }{\mathrm{\partial }{p}^2}.$$

The FPK equation can be solved exactly by assuming the following Gaussian form for the probability distribution

(55) $$\rho =\frac{1}{Z}{e}^{-(a{x}^2+b{p}^2+2cxp)/2},$$

where the parameters a, b, c depend on time and

(56) $$Z=\frac{2\pi }{\sqrt{ab-c^2}}.$$

The solution is given in the appendix 16, where we find the parameters a, b, and c as functions of time.

From the probability distribution (55), we determine the covariances,

(57) $$⟨x^2⟩=\frac{b}{ab-c^2},$$

(58) $$⟨p^2⟩=\frac{a}{ab-c^2},$$

(59) $$⟨xp⟩=\frac{-c}{ab-c^2},$$

and other properties. The energy is

(60) $$U=\frac{1}{2m}⟨p^2⟩+\frac{K}{2}⟨x^2⟩,$$

and the entropy is found by using its definition and is

(61) $$S=k+\ln 2\pi -\frac{1}{2}\ln (ab-c^2).$$

To determine $dS/dt$, $\mathrm{\Pi }$ and $\mathrm{\Phi }$, we need to find J, which is defined by (36). In the present case it reads

(62) $$J=\frac{m\gamma }{\beta }(cx+bp-\frac{\beta }{m}p)\rho .$$

From J we find

(63) $$\frac{dS}{dt}=\frac{km\gamma }{\beta }b-k\gamma ,$$

(64) $$\mathrm{\Phi }=\frac{k\beta \gamma }{m}⟨p^2⟩-k\gamma ,$$

(65) $$\mathrm{\Pi }=\frac{k\beta \gamma }{m}⟨p^2⟩+\frac{km\gamma }{\beta }b-2k\gamma ,$$

and it becomes clear that $dS/dt=\mathrm{\Pi }-\mathrm{\Phi }$.

From the asymptotic values of the parameter a, b, and c, given in the appendix 16, which are $a=K\beta$, $b=\beta /m$, and $c=0$, we find the equilibrium values of the various quantities which are $⟨xp⟩=0$,

(66) $$\frac{1}{2m}⟨p^2⟩=\frac{1}{2}K⟨x^2⟩=\frac{1}{2\beta }=\frac{1}{2}kT,$$

(67) $$U=\frac{1}{\beta }=kT,$$

(68) $$S=k+\ln 2\pi -\frac{1}{2}\ln \frac{K{\beta }^2}{m},$$

$dS/dt=0$, $\mathrm{\Phi }=0$, and $\mathrm{\Pi }=0$. As the parameters a, b, and c decay exponentially to their asymptotic values, so do the properties obtained above.

We remark that the probability distribution approaches the equilibrium probability distribution

(69) $${\rho }_e=\frac{1}{Z}{e}^{-\beta H},$$

where H is the energy function (53), and the system thermalizes properly.

8. Quantum evolution equation

To extend the stochastic approach developed above to the quantum case we need to provide a quantum version of the evolution equation. One way of setting up the quantum version is to use a procedure known as canonical quantization, which amounts to replace the Poisson brackets of classical mechanics by a quantum commutator. For the case of just one degree of freedom that we are considering here, the Poisson brackets between A and B are given by

(70) $$\{A,B\}=\frac{\mathrm{\partial }A}{\mathrm{\partial }x}\frac{\mathrm{\partial }B}{\mathrm{\partial }p}-\frac{\mathrm{\partial }B}{\mathrm{\partial }x}\frac{\mathrm{\partial }A}{\mathrm{\partial }p}.$$

The canonical quantization is obtained by performing the replacement

(71) $$\{A,B\}\to \frac{1}{i\hslash }[\hat{A},\hat{B}],$$

where $\hslash$ is the Planck constant and the quantities $\hat{A}$ and $\hat{B}$ on the right are understood as quantum operators, and $[\hat{A},\hat{B}]=\hat{A}\hat{B}-\hat{B}\hat{A}$.

From the quantization rule above, we obtain two useful rules. The first is obtained by setting $A=x$ in the Poisson brackets, which gives

(72) $$\{x,B\}=\frac{\mathrm{\partial }B}{\mathrm{\partial }p}.$$

Using the quantization rule, we obtain

(73) $$\frac{\mathrm{\partial }B}{\mathrm{\partial }p}\to \frac{1}{i\hslash }[\hat{x},\hat{B}].$$

In an analogous way, if we set $B=p$, we find

(74) $$\{A,p\}=\frac{\mathrm{\partial }A}{\mathrm{\partial }x},$$

and using the quantization rule, we obtain

(75) $$\frac{\mathrm{\partial }A}{\mathrm{\partial }x}\to \frac{1}{i\hslash }[\hat{A},\hat{p}].$$

We remark that $\hat{x}$ and $\hat{p}$ on the right-hand sides of (73) and (75) are quantum operators representing the position and the momentum of a particle, and we recall that according to quantum mechanics $[\hat{x},\hat{p}]=i\hslash$.

The last two rules are useful in the transformation of a differential equation such as the FPK equation into a quantum equation. It should be remarked however that the equation obtained by this procedure is not a mathematical derivation of the quantum equation from the classical equation. In fact, the opposite is true. From the quantum equation one reaches the classical equation by taking the classical limit. Thus, the quantization rules should be used as a guide to find a quantum equation which at the end should be introduced as a postulate.

It is usual to use the hat symbol to denote a quantum operator as we have done above. But from now on we will drop the hat symbol and denote an operator by a letter without the hat. Thus the position and momentum operator, for instance, will be denoted by x and p.

Let us consider the FPK equation in the form (8). According to the quantization rules, the quantum evolution equation is

(76) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=\frac{1}{i\hslash }[H,\rho ]-\frac{1}{i\hslash }[x,J],$$

where now, $\rho$, H, and J are quantum operators. Since the quantum operators can be represented by matrices most properties of quantum operators are better understood if stated in terms of matrices. For instance the density operator $\rho$, which is the quantum version of the probability density distribution, holds the following property: its diagonal elements are nonnegative and the sum of the diagonal elements, its trace, equals the unity,

(77) $$\mathrm{T}\mathrm{r}\rho =1.$$

The operator H is the quantum energy function, or the Hamiltonian, given by

(78) $$H=\frac{1}{m}{p}^2+V,$$

where V is a function of x. We remark that without J the equation reduces to the quantum Liouville equation of quantum statistical mechanics.

The property (77) is the analogue of the normalization of the probability density that we used before and thus it should be preserved in time. To see that this is the case let us take the trace of the right-hand side of equation (76). There are two terms to be considered and both vanish because each one is the trace of a commutator, and the trace of a commutator vanishes. Therefore, the left hand side, which is the time derivative of the trace of $\rho$, vanishes and $\rho$ should be constant in time.

The trace of a commutator vanishes because $\mathrm{T}\mathrm{r}(AB)=\mathrm{T}\mathrm{r}(BA)$. This is a particular case of the cyclic property of the trace $\mathrm{T}\mathrm{r}(ABC)=\mathrm{T}\mathrm{r}(BCA)$. This cyclic property allows us to write the following property

(79) $$\mathrm{T}\mathrm{r}([A,B]C)=\mathrm{T}\mathrm{r}(A[B,C]),$$

that we will employ further on.

9. Energy and entropy

The average $U=⟨H⟩$ of the energy function H with respect to the density operator $\rho$ is given by

(80) $$U=\mathrm{T}\mathrm{r}(H\rho ),$$

and is the quantum analog of the integral in (17). Deriving it with respect to time,

(81) $$\frac{dU}{dt}=\mathrm{T}\mathrm{r}\left(H\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}\right),$$

and using the evolution equation (76), we find

(82) $$\frac{dU}{dt}=-\frac{1}{i\hslash }\mathrm{T}\mathrm{r}(H[x,J]),$$

where the term involving the commutator $[H,\rho ]$ vanishes owing to the property (79). Using again this same property we get

(83) $$\frac{dU}{dt}=\frac{1}{i\hslash }\mathrm{T}\mathrm{r}([x,H]J).$$

The right-hand side is interpreted as the heat flux into the system

(84) $${\mathrm{\Phi }}_q=\frac{1}{i\hslash }\mathrm{T}\mathrm{r}([x,H]J),$$

and

(85) $$\frac{dU}{dt}={\mathrm{\Phi }}_q.$$

The definition of entropy for the quantum case is that introduced by von Neumman,

(86) $$S=-k\mathrm{T}\mathrm{r}(\rho \ln \rho ),$$

and corresponds to the extension of the Gibbs entropy to the quantum case. Deriving this equation with respect to time we find

(87) $$\frac{dS}{dt}=-k\mathrm{T}\mathrm{r}(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}\ln \rho ),$$

where the term involving the derivative of $\ln \rho$ vanishes in view of the normalization property (77). Using the evolution equation, we find

(88) $$\frac{dS}{dt}=\frac{k}{i\hslash }\mathrm{T}\mathrm{r}([x,J]\ln \rho ),$$

where again the term involving the commutator $[H,\rho ]$ vanishes owing to the property (79).

10. Rate of entropy production

The right-hand side of the equation (88) is not equal to to the flux of entropy because the entropy is not a conserved quantity. There is another source of entropy which comes from dissipation inside the system. Thus as before the time variation of the entropy has two term,

(89) $$\frac{dS}{dt}=\mathrm{\Pi }-\mathrm{\Phi },$$

where $\mathrm{\Pi }$ is the rate of entropy production and $\mathrm{\Phi }$ is the flux of entropy from the system to the outside.

Guided by the classical version, the production of entropy is defined as follows. The quantum irreversible current J has not been defined yet but it is expressed in terms of the density operator, that is, $J(\rho )$. Let us denote by ${\rho }_0$ the quantity such that $J({\rho }_0)=0$. If the commutator of ${\rho }_0$ with H also vanishes then ${\rho }_0$ is identified as the density at thermodynamic equilibrium. However, here we do not demand that this commutator vanishes. The rate of entropy production is defined as

(90) $$\mathrm{\Pi }=\frac{k}{i\hslash }\mathrm{T}\mathrm{r}\{[x,J](\ln \rho -\ln {\rho }_0\},$$

and it becomes clear that $\mathrm{\Pi }$ vanishes whenever J vanishes. Taking into account (88) and (89), the expression for the flux of entropy $\mathrm{\Phi }$ is

(91) $$\mathrm{\Phi }=-\frac{k}{i\hslash }\mathrm{T}\mathrm{r}\{[x,J]\ln {\rho }_0\}.$$

Let us assume that the irreversible current J, which has not yet been specified, is so defined that the evolution equation describes a system in contact with a heat reservoir at temperature T. In thermodynamic equilibrium J vanishes and ${\rho }_0$ should be identified as corresponding to the Gibbs canonical distribution which in the quantum case reads

(92) $${\rho }_0=\frac{1}{Z}{e}^{-\beta H},$$

where

(93) $$Z=\mathrm{T}\mathrm{r}(e^{-\beta H}).$$

Replacing ${\rho }_0$ in the expressions for $\mathrm{\Pi }$ and $\mathrm{\Phi }$, we find

(94) $$\mathrm{\Pi }=\frac{k}{i\hslash }\mathrm{T}\mathrm{r}\{[x,J](\ln \rho +\beta H\},$$

(95) $$\mathrm{\Phi }=\frac{k\beta }{i\hslash }\mathrm{T}\mathrm{r}([x,J]H).$$

Now let us compare equations (95) and (84). We see that $\mathrm{\Phi }$ and ${\mathrm{\Phi }}_q$ are related by

(96) $$\mathrm{\Phi }=-\frac{1}{T}{\mathrm{\Phi }}_q,$$

which is the thermodynamic relation that should exist between the flux of entropy $\mathrm{\Phi }$ and the heat flux ${\mathrm{\Phi }}_q$ when a system is in contact with a heat reservoir at a temperature T. Other thermodynamic relations that we have obtained for the classical case, such those given by equations (39) and (40), can also be shown to be valid in the quantum case.

11. Irreversible current

The quantum irreversible current J has not yet been specified. Here we will choose it by applying the rules of the canonical quantization to the expression (9). The form chosen for the irreversible current is

(97) $$J=\frac{1}{2}(D^{†}\rho +\rho D)-\frac{\mathrm{\Gamma }}{2}\frac{1}{i\hslash }[x,\rho ],$$

where D is a quantum operator representing the dissipative force and $\mathrm{\Gamma }$ is a real constant as in the classical case. Notice that we have written a symmetrized form for the product of the dissipative force and the density operator. In one of the products, we have used the Hermitian conjugate of D, denoted $D^{†}$, so that the whole expression is Hermitian.

The matrix $A^{†}$ that represents the Hermitian conjugate of an operator is obtained from the matrix A that represent the operator by transposing the elements of the matrix and taking the complex conjugate of each element. If $A_{ij}$ and $(A^{†}{)}_{ij}$ represent an element of these two matrices then $(A^{†}{)}_{ij}=(A_{ji}{)}^{\ast }$. A Hermitian operator is the operator which is equal to its Hermitian conjugate. An important property of such an operator is that its eigenvalue are real.

If we wish to describe a system that at long times reaches the thermodynamic equilibrium then D and $\mathrm{\Gamma }$ should have a relation such that J vanishes when $\rho$ is the equilibrium density operator ${\rho }_e$. Imposing the vanish of J when $\rho$ is replaced by ${\rho }_e$ we find the condition

(98) $$\frac{1}{2}(D^{†}{\rho }_e+{\rho }_{e}D)=\frac{\mathrm{\Gamma }}{2}\frac{1}{i\hslash }[x,{\rho }_e].$$

The solution for D is

(99) $$D=\frac{\mathrm{\Gamma }}{2}\frac{1}{i\hslash }{\rho }_e^{-1}[x,{\rho }_e]=\frac{\mathrm{\Gamma }}{2}\frac{1}{i\hslash }({\rho }_e^{-1}x{\rho }_e-x),$$

which is understood as the relation between dissipation, described by D, and noise or fluctuations, described by $\mathrm{\Gamma }$. In equilibrium there must be a relation between dissipation and fluctuations.

That (99) is a solution can be verified by substitution and using the property that the Hermitian conjugate of a product of two operators equals the product of the conjugate of each operator in the inverse order. In the present case, $({\rho }_{e}D{)}^{†}=D^{†}{\rho }_e$ because ${\rho }_e$ is Hermitian.

Next we seek for J that could describe the contact of the system with a heat reservoir. In this case

(100) $${\rho }_e=\frac{1}{Z}{e}^{-\beta H},$$

which replace in (99) gives

(101) $$D=\frac{\mathrm{\Gamma }}{2}\frac{1}{i\hslash }(e^{\beta H}x{e}^{-\beta H}-x).$$

This form of dissipation is certainly not the form of the classical dissipation found above which is proportional to the momentum. However, at high temperatures this is the case. If we expand the terms between parentheses in the right-hand side of (101) up to terms of order $\beta$, we find $D=-\gamma p$ where $\gamma =\mathrm{\Gamma }\beta /2m$, or

(102) $$\mathrm{\Gamma }=\frac{2\gamma m}{\beta }.$$

For an arbitrary temperature, the dissipation, according to the present approach, is not proportional to the momentum and is given by (101), which we write, by using (102), as

(103) $$D=-\gamma g,$$

where

(104) $$g=-\frac{m}{i\hslash \beta }(e^{\beta H}x{e}^{-\beta H}-x).$$

Replacing the irreversible current (97) into the evolution equation (76) we may write it in the more explicit form

(105) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=\frac{1}{i\hslash }[H,\rho ]+\frac{\gamma }{2i\hslash }[x,g^{†}\rho +\rho g]-\frac{\gamma m}{\beta {\hslash }^2}[x,[x,\rho ]],$$

which we call the quantum FPK equation.

12. Quantum harmonic oscillator

For a quantum Harmonic oscillator the energy function is

(106) $$H=\frac{1}{2m}{p}^2+\frac{1}{2}m{\omega }^2{x}^2,$$

where $\omega$ is the frequency of oscillations. To find the solution of the quantum FPK equation from which we may determine the thermodynamic properties it is necessary to know the explicit expression of the dissipation force $D=-\gamma g$, that is, we need to know g as a function of x and p. For the harmonic oscillator we show in the appendix 16 that

(107) $$g=ap+ibx,$$

where a and b are real numbers given by

(108) $$a=\frac{1}{\beta \hslash \omega }\sinh \beta \hslash \omega ,b=\frac{m}{\beta \hslash }(\cosh \beta \hslash \omega -1).$$

A solution of the quantum FPK equation (105) can be obtained by a method similar to that used in the classical case, which is to assume a solution of the form

(109) $$\rho =\frac{1}{Z}{e}^{-(a{x}^2+b{p}^2+cxp+cpx)/2},$$

where a, b, and c are real constant that depends on time. Here we will limit ourselves to write down the equations for the covariances and determine their asymptotic values, which are the values at thermodynamic equilibrium. The time dependent solution can be found in the reference [30]. Multiplying the quantum FPK equation successively by $x^2$, $p^2$, and $xp$, and taking the trace, we obtain the following equations for the covariances

(110) $$\begin{array}{rcl}& & \frac{d}{dt}⟨p^2⟩=-m{\omega }^2(⟨px⟩+⟨xp⟩)\\ & & +\hslash b\gamma -2a\gamma ⟨p^2⟩+\frac{2\gamma m}{\beta },\end{array}$$

(111) $$\frac{d}{dt}⟨x^2⟩=\frac{1}{m}(⟨px⟩+⟨xp⟩),$$

(112) $$\frac{d}{dt}⟨xp⟩=\frac{1}{m}⟨p^2⟩-m{\omega }^2⟨x^2⟩-\frac{a\gamma }{2}(⟨px⟩+⟨xp⟩).$$

The equation for $⟨px⟩$ is not needed because $⟨xp⟩-⟨px⟩=i\hslash$.

At the stationary state, which is a state of thermodynamic equilibrium we find

(113) $$⟨xp⟩=-⟨xp⟩=\frac{i\hslash }{2},$$

(114) $$\frac{1}{2m}⟨p^2⟩=\frac{1}{2}m{\omega }^2⟨x^2⟩=\frac{1}{2}\hslash \omega (\frac{1}{e^{\beta \hslash \omega }-1}+\frac{1}{2}).$$

From these results one reaches the expected expression for the average energy of a quantum oscillator,

(115) $$⟨H⟩=\hslash \omega (\frac{1}{e^{\beta \hslash \omega }-1}+\frac{1}{2}).$$

We remark that the probability distribution approaches the equilibrium probability distribution

(116) $${\rho }_e=\frac{1}{Z}{e}^{-\beta H},$$

where H is the quantum energy function (106), and the system thermalizes properly.

13. Multiple degrees of freedom

Up to this point, we have considered systems with just one degree of freedom. Here we wish to consider the case of a system with multiple degrees of freedom. The derivations of the results for the present case parallel those obtained for one degree of freedom and will not be shown in full detail. We restrict ourselves to the classical case but the quantum case can be obtained in a way similar to the case of one degree of freedom and can be found in reference [30].

An appropriate treatment of a system with many degrees of freedom begins with the generalization of the FKP equation (44) for this case. As we have seen, this equation describes a system in contact with a heat reservoir and is subject to an external force. The generalization that we give below is the one appropriate to describe a system in contact with multiple reservoirs at distinct temperatures and subject to several forces. We denote by $x_i$ a Cartesian component of the position and by $p_i$ the respective component of the momentum related to a certain degree of freedom. The energy function is a sum of kinetic energy and a potential energy,

(117) $$H=\sum _i\frac{p_i^2}{2m}+V.$$

The FPK equation, which governs the time evolution of the probability density $\rho$, now reads

(118) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=\{H,\rho \}-\sum _i\frac{\mathrm{\partial }{F}_i^e\rho }{\mathrm{\partial }{p}_i}-\sum _i\frac{\mathrm{\partial }{J}_i}{\mathrm{\partial }{p}_i},$$

where $F_i^e$ are the Cartesian components of the external force, which may be nonconservative and time dependent, and

(119) $$J_i=D_i\rho -\frac{{\mathrm{\Gamma }}_i}{2}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }{p}_i}$$

are the components of the dissipative probability current. We choose the dissipation force $D_i$ to be of the usual form $D_i=-\gamma p_i$ and ${\mathrm{\Gamma }}_i=2mk{T}_i$ so that we may interpret the FPK equation of describing a system in contact with several heat reservoirs at temperatures $T_i$. Therefore,

(120) $$J_i=-\gamma (p_i\rho +mk{T}_i\frac{\mathrm{\partial }\rho }{\mathrm{\partial }{p}_i}).$$

In the absence of external forces and if all heat baths have the same temperature $T_i=T$, then the stationary state is a state of thermodynamic equilibrium because in this case $J_i$ vanishes for all i. Indeed, if we replace the Gibbs probability distribution

(121) $${\rho }_e=\frac{1}{Z}{e}^{-\beta H},$$

where $\beta =1/kT$, in the expression for $J_i$ we see that it vanishes. We remark that the Poisson brackets also vanish.

The time variation of the energy $U=⟨H⟩$ has the same form of equation (45),

(122) $$\frac{dU}{dt}={\mathrm{\Phi }}_q-{\mathrm{\Phi }}_w,$$

but now ${\mathrm{\Phi }}_q$ is a sum of the heat fluxes coming from each reservoir, that is,

(123) $${\mathrm{\Phi }}_q=\sum _i{\mathrm{\Phi }}_{qi},$$

where each heat flux has the form (21), with J replaced by $J_i$,

(124) $${\mathrm{\Phi }}_{qi}=\int J_i\frac{\mathrm{\partial }H}{\mathrm{\partial }{p}_i}dxdp.$$

Using the relation $\mathrm{\partial }H/\mathrm{\partial }{p}_i=p_i/m$ and performing an integration by parts, the heat flux can be written as

(125) $${\mathrm{\Phi }}_{qi}=-\gamma (\frac{1}{m}⟨p_i^2⟩-k{T}_i).$$

It becomes clear that if $k{T}_i/2$ is larger than the average kinetic energy $⟨p_i^2⟩/2m$ then ${\mathrm{\Phi }}_{qi}>0$, and heat flows into the system, otherwise, heat flows from the system to the heat reservoir.

The expression for the power ${\mathrm{\Phi }}_w$ is similar to that of equation (46) but now there is a sum over all components of the force

(126) $${\mathrm{\Phi }}_w=-\sum _i\int F_i^e\rho \frac{\mathrm{\partial }H}{\mathrm{\partial }{p}_i}dxdp.$$

Using again the relation $\mathrm{\partial }H/\mathrm{\partial }{p}_i=p_i/m$ and performing an integration by parts, the power can be written as

(127) $${\mathrm{\Phi }}_w=-\frac{1}{m}\sum _i⟨F_i^e{p}_i⟩.$$

It is useful to understand that the forces $F_i^e$ are being acted by an external agent which is a power device. The role played by the power device in relation to the transfer of mechanical work is analogous to the role played by the heat reservoir in relation to the transfer of heat. We see from (127) that the power has the usual form of a force multiplied by the velocity. If $F_i^e$ and $p_i$ have the same sign, work is performed by the power device onto the system, otherwise, the system performs work on the power device.

The heat flux and the power in (122) might be understood as functions of the time so that we may write ${\mathrm{\Phi }}_q=dQ/dt$ and ${\mathrm{\Phi }}_w=dW/dt$ in which case equation (122) reduces to the form

(128) $$dU=dQ-dW,$$

which is the usual way of writing the conservation of energy. However, it should be remarked that both $dQ$ and $dW$ are not in general exact differential, although $dU$ is. The concepts of exact and inexact differentials are commented in the appendix 16.

The variation of entropy with time is

(129) $$\frac{dS}{dt}=\mathrm{\Pi }-\mathrm{\Phi },$$

where the rate of entropy production $\mathrm{\Pi }$ and the entropy flux $\mathrm{\Phi }$ are generalization of the equations (31) and (35) for the present case,

(130) $$\mathrm{\Pi }=\frac{1}{m\gamma }\sum _i\frac{1}{T_i}\int \frac{J_i^2}{\rho }d\xi ,$$

where the integration extends over the space of states, and

(131) $$\mathrm{\Phi }=\sum _i\frac{\gamma }{T_i}(\frac{1}{m}⟨p_i^2⟩-k{T}_i).$$

In view of the equation (125), we see that the flux of entropy is related to the heat fluxes by

(132) $$\mathrm{\Phi }=-\sum _i\frac{{\mathrm{\Phi }}_{qi}}{T_i}.$$

Using this relation, we may derive again all the results involving the free energy obtained in section 6, including the inequality 52.

14. Nonequilibrium steady state

Here we wish to show by an example that the present approach is indeed capable of describing systems displaying nonequilibrium steady states. That is, for long times the systems reaches a stationary state in which the irreversible currents are nonzero and entropy is permanently being created. One way of setting a system in a nonequilibrium steady state is to place the system in contact with heat reservoirs with different temperatures. Another possibility is to place the system under the action of a power device. The two possibilities are embodied in the development made in the previous section.

We examine a system with two degrees of freedom. The potential energy is harmonic and the energy function is

(133) $$H=\frac{1}{2m}(p_1^2+p_2^2)+\frac{1}{2}K(x_1^2+x_2^2)-L{x}_1{x}_2,$$

so that the conservative forces are linear, $F_1=-K{x}_1+L{x}_2$ and $F_2=-K{x}_2+L{x}_1$. In addition to the conservative forces, the system is under the action of a nonconservative external force given by

(134) $$F_1^e=-c{x}_2,F_2^e=c{x}_1.$$

According to equation (127),

(135) $${\mathrm{\Phi }}_w=\frac{c}{m}⟨x_2{p}_1⟩-\frac{c}{m}⟨x_1{p}_2⟩,$$

which is minus the power performed by the power device. Each degree of freedom is understood as being coupled to one of the two heat reservoirs at the temperatures $T_1$ and $T_2$. The fluxes of heat associated to each reservoir are given by (125) and are

(136) $${\mathrm{\Phi }}_{q1}=-\frac{\gamma }{m}(⟨p_1^2⟩-mk{T}_1),$$

(137) $${\mathrm{\Phi }}_{q2}=-\frac{\gamma }{m}(⟨p_2^2⟩-mk{T}_2).$$

From the two heat fluxes we determine the entropy flux by means of the relation (132)

(138) $$\mathrm{\Phi }=-\frac{{\mathrm{\Phi }}_{q1}}{T_1}-\frac{{\mathrm{\Phi }}_{q2}}{T_2}.$$

From now on we wish to determine the above quantities in the steady state. In the steady state $\mathrm{\Pi }=\mathrm{\Phi }$ and in view of the relation (138), it suffices to determine the heat fluxes ${\mathrm{\Phi }}_{qi}$ and the power ${\mathrm{\Phi }}_w$. To find these quantities we need the covariances $⟨x_i{p}_j⟩$ and $⟨p_i^2⟩$ at the steady state. Therefore, we should solve the FPK equation (118), with the energy function H given by (133) and external forces given by (134).

In view of the fact that the forces, conservative and nonconservative, are linear, it is possible to solve the FPK equation exactly. The solution is carried out in the appendix 16 where the covariances at the stationary state are determined. Replacing the covariances in the expressions (135), (136), and (137), we find

(139) $${\mathrm{\Phi }}_w=-\frac{2c}{m}\mathcal{C},$$

(140) $${\mathrm{\Phi }}_{q1}=-\frac{c-L}{m}\mathcal{C},{\mathrm{\Phi }}_{q2}=-\frac{c+L}{m}\mathcal{C},$$

where

(141) $$\mathcal{C}={\mathcal{C}}_0(L\mathrm{\Delta }T+2cT),$$

$T=(T_1+T_2)/2$, $\mathrm{\Delta }T=T_1-T_2$, and

(142) $${\mathcal{C}}_0=\frac{m\gamma k}{2(m{\gamma }^{2}K+L^2-c^2)}.$$

The range of values of c are such that the denominator be positive.

In the steady state the production of entropy equals the entropy flux which, from (138) is given by

(143) $$\mathrm{\Pi }=\mathrm{\Phi }=\frac{{\mathcal{C}}_0}{T_1{T}_{2}m}(2cT+L\mathrm{\Delta }T{)}^2,$$

which is clearly nonnegative.

Let us analyze the results we have just found for $L>0$. The various possibilities for the heat fluxes and power are shown in table 1. In all cases where $c\ne 0$, except one, energy is being dissipated, that is, work is performed onto the system $({\mathrm{\Phi }}_w<0)$ which in turn releases it in the form of heat to one or to both heat reservoirs. The exception is the case in which the system perform works $({\mathrm{\Phi }}_w>0)$ in which case heat flows from the hotter reservoir to the colder reservoir through the system, and the whole system functions as a heat engine.

Table 1
Heat fluxes and power for the case of the system defined by the energy function (133) and by the external force (134), where $\mathrm{\Delta }T=T_1-T_2$. The convention for the fluxes are: if ${\mathrm{\Phi }}_{q1}$ is positive, heat flows from the reservoir 1 to the system, if ${\mathrm{\Phi }}_{q2}$ is positive, heat flows from the reservoir 2 to the system, if ${\mathrm{\Phi }}_w$ is positive, the system performs work on to the external device. Notice that ${\mathrm{\Phi }}_w={\mathrm{\Phi }}_{q1}+{\mathrm{\Phi }}_{q2}$. We are considering here $L>0$ and $|c|<L$, and that $T_1\ge T_2$.

$\mathrm{\Delta }T$	c	$\mathcal{C}$	${\mathrm{\Phi }}_w$	${\mathrm{\Phi }}_{q1}$	${\mathrm{\Phi }}_{q2}$
0	+	+	-	+	-
0	-	-	-	-	+
+	0	+	0	+	-
+	+	+	-	+	-
+	-	+	+	+	-
+	-	-	-	-	+

15. Path integral

The approach to the stochastic thermodynamics that we have developed here is based on the FPK equation which is understood as an equation that governs the time evolution of the probability distribution $\rho$. The solution of the evolution equation gives $\rho$ at any instant of time. For convenience here we denote a state by $\xi$ which is understood as the collection of positions and momenta of the particles of the system.

If we solve the FPK equation considering that at the initial time $t_0$ it was in a certain state ${\xi }_0$, then $\rho (\xi ,t)$ is understood as the conditional probability of finding the system in state $\xi$ at time t, given that it was in the state ${\xi }_0$ at time $t_0$, and we denote it by $P(\xi ,t|{\xi }_0,t_0)$.

Let us consider now a discretized trajectory, that is, a trajectory for which the system is at state ${\xi }_0$ at an initial time $t_0$, in ${\xi }_1$ at time $t_1$, in ${\xi }_2$ at time $t_2$, …, and in ${\xi }_n$ at the final time $t_n$. The probability of the occurrence of this discretized trajectory is a successive product of

(144) $$P({\xi }_{\ell },t_{\ell }|{\xi }_{\ell -1},t_{\ell -1}),$$

from $\ell =1$ until $\ell =n$, multiplied by the probability $P({\xi }_0,t_0)$. Omitting the reference to the instants of time, the trajectory probability is

(145) $$P({\xi }_n|{\xi }_{n-1})\dots P({\xi }_2|{\xi }_1)P({\xi }_1|{\xi }_0)P({\xi }_0).$$

From now on we consider that the time intervals between two successive instants of time are the same and equal to $\tau$. It is understood that $\tau$ is small enough so that the trajectory approaches a continuous trajectory.

Generally speaking the probability of a trajectory is a joint probability, which we denote by

(146) $$P({\xi }_n,{\xi }_{n-1},\dots ,{\xi }_2,{\xi }_1,{\xi }_0).$$

The identification of (146) with (145) defines a type of stochastic dynamics associated with the name of Markov and the approach we are using here, with the FKP equation as the evolution equation, is thus a Markovian stochastic dynamics.

The probability distribution (146) is a joint probability distribution. If we integrate in all variables except one one them, we find the probability of this variable. For instance, if we integrate in ${\xi }_0,{\xi }_1,\dots ,{\xi }_{n-1}$, we find the probability of ${\xi }_n=\xi$ at time $t_n=t$, which is $\rho (\xi ,t)$, that is,

(147) $$\rho (\xi )=\int \cdots \int P(\xi ,{\xi }_{n-1},\dots ,{\xi }_0)d{\xi }_{n-1}\dots d{\xi }_0.$$

In which circumstance should we use the path probability (146)? If we wish to find the average U of the energy function $H(\xi )$ at time t, for instance, it suffices to use the probability distribution $\rho (\xi ,t)$. However, if we wish to find the average of the mechanical work, we should use the path probability because the work is a path integral.

Let L be the work along a certain trajectory of the force with components $f_i$,

(148) $$L=\sum _i{\int }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}{f}_{i}d{x}_i,$$

where the index ‘path’ serves to remember that the integral is an integral along a certain trajectory. In a discretized form, the path integral is written the sum

(149) $${\int }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}{f}_{i}d{x}_i=f_i^0{a}_i^0+f_i^1{a}_i^1+\dots +f_i^n{a}_i^n,$$

where $f_i^{\ell }$ is the value of $f_i$ at the $\ell$-th step, and $a_i^{\ell }$ is the increment in $x_i$ at the $\ell$-th step. In this discretized form we see clearly that the path integral depends on ${\xi }_0,{\xi }_1,\dots ,{\xi }_n$, and to find its average we should use the path probability (146). This amounts to multiply the right-hand side of (149) by the right-hand side of (146) and integrate in all variables, ${\xi }_0,{\xi }_1,\dots ,{\xi }_n$. The result of this procedure is indicated by an index ‘path’ in the signs of the average. Therefore, the average of L which we call W is written as

(150) $$W=⟨L{⟩}_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}.$$

Although we call work both L and W, it should be understood that L is the actual work and W its average.

The path integral (148) can be written as a time integral of the power as is well known. A trajectory may be defined parametrically by given ${\xi }_i$, and thus $x_i$ and $p_i$, as functions of this parameter which we take to be the time t. In terms of this parameter, the integral (148) becomes a time integral,

(151) $$L=\sum _i{\int }_{t_1}^{t_2}{f}_i\frac{p_i}{m}dt,$$

where we have replaced $d{x}_i/dt$ by the velocity $p_i/m$. This expression is written as

(152) $$L={\int }_{t_1}^{t_2}\varphi dt,$$

where $\varphi$ is the power at time t, and given by

(153) $$\varphi =\sum _i{f}_i\frac{p_i}{m}.$$

Taking the average of the expression (152) we find

(154) $$⟨L{⟩}_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}={\int }_{t_1}^{t_2}⟨\varphi ⟩dt,$$

where in the right-hand side the average is the usual average taken by the use of the probability density $\rho (\xi ,t)$ at time t because $\varphi$ depends only on $\xi$ at time t. Thus the average over a path integral is transformed into an average over the ordinary probability density. The integrand on the right-hand side of (154) is understood as the average power,

(155) $${\mathrm{\Phi }}_w=⟨\varphi ⟩=\int \sum _i{f}_i\frac{p_i}{m}\rho d\xi ,$$

which coincides with the power of the external force given by the equivalent forms (127) or (126), if we recall that $f_i=-F_i^e$. Since W is the average of L, we may write

(156) $$W={\int }_{t_1}^{t_2}{\mathrm{\Phi }}_{w}dt.$$

The result (156) or its equivalent form (154) where L and $\varphi$ are given by (148) and (153), respectively, can immediately be generalized by replacing $f_i(\xi )$ by any other function of $\xi$. Let us suppose that it is replaced by $J_i/\rho$ where $J_i$ is the irreversible current given by (119)

(157) $$J_i=D_i\rho -\frac{{\mathrm{\Gamma }}_i}{2}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }{p}_i}.$$

The path integral of this quantity is

(158) $$\mathrm{\Psi }=\sum _i{\int }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}\frac{J_i}{\rho }d{x}_i,$$

and the associated flux according to the result above is

(159) $${\mathrm{\Phi }}_q=\frac{1}{m}\sum _i⟨\frac{1}{\rho }{J}_i{p}_i⟩=\frac{1}{m}\sum _i\int J_i{p}_{i}d\xi ,$$

which is the heat flux as given by (124). The quantity

(160) $$Q={\int }_{t_1}^{t_2}{\mathrm{\Phi }}_{q}dt$$

is thus the heat exchanged between the two instants of time, and according to the results above may be written as the path average

(161) $$Q=⟨\mathrm{\Psi }{⟩}_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}.$$

Let us consider that a system in contact with a heat bath at a temperature T is acted by external forces during a certain interval of time. The following equality relating the free energy to the work performed by the system during this interval of time has been shown to be valid [6],

(162) $$e^{-\beta \mathrm{\Delta }F}=⟨e^{+\beta L}{⟩}_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}.$$

Therefore, if the right-hand side of equation (162) is measured, we may obtain $\mathrm{\Delta }F$. Using the inequality $⟨e^{\eta }⟩\ge e^{⟨\eta ⟩}$, valid for any random variable $\eta$, we find

(163) $$\mathrm{\Delta }F\le -⟨L{⟩}_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{h}}=-W,$$

which should be compared with the equation (52).

Before we end this section it is appropriate to place a discussion on the experimental measurement of the several quantities presented in the theory. The quantities that are mensurable are those that we call state functions, that is, quantities that are functions of the random variables, in the present case the positions and momenta, and themselves random variables. A experimental result obtained for a state function E, be it the value of one trial or the arithmetic average of several trials, should be compared with the average predicted by the theory which is written as

(164) $$⟨E⟩=\int E(\xi )\rho (\xi )d\xi ,$$

such as the energy, or as a path average as is the case of the mechanical work.

Let us consider the case of the entropy which is

(165) $$S=-k\int \rho (\xi )\ln \rho (\xi )d\xi ,$$

which sometimes is written as

(166) $$S=⟨-k\ln \rho ⟩.$$

Although one may write in this form, this expression is merely an abbreviation of the expression on the right-hand side of (165) and cannot be understood as the average of a state function merely because $-k\ln \rho$ is not a state function. Although, sometimes $-k\ln \rho$ is called an instantaneous entropy, it is not a mensurable quantity. This point can be better understood if we try to calculate the entropy from a Monte Carlo simulation. One immediately realizes that it is impossible to determine $\ln \rho$ along a Monte Carlo run and an alternative should be used.

The observation made above with respect to the average (166) is extended to the average in (161) because $J_i/\rho$ is not a state function and it is not a mensurable quantity. This is paradoxical because no one denies that heat Q is mensurable. However, a moment of reflection will reveal that heat is measured through the work dissipated and not as the quantity $\mathrm{\Psi }$ above.

16. Discussion and conclusion

We have developed an approach to the stochastic thermodynamics based on the use of the FPK equation, which governs the time evolution of the probability distribution. The main feature of the approach in addition to the evolution equation is the assignment of an energy function, the definition of entropy and the introduction of an expression for the rate of production of entropy. According to the approach, these quantities are well defined quantities in equilibrium or out of equilibrium. This is in contrast to other quantities of thermodynamics such as the temperature, which is defined only when the system is in thermodynamic equilibrium.

The evolution equation contains the mechanism of dissipation and stochastic fluctuation or noise which leads the system toward equilibrium, if an appropriate relation exists between dissipation and noise. The mechanism is included in the irreversible current by the term containing the dissipative force and the term containing the quantity $\mathrm{\Gamma }$, which is a measure of the noise. If the thermodynamic equilibrium sets in, the irreversible current vanishes. In out of equilibrium the irreversible current is nonzero and the production of entropy, which is related to the square of the irreversible currents is greater than zero. The rate of entropy production is thus a measure of the deviation of a system from thermodynamic equilibrium and of the irreversibility.

We have considered systems with a continuous space of states in which case the appropriate evolution equation is the FPK equation. However, the present approach can be extended to a discrete space of states in which case the evolution equation is called a master equation [29, 31]. It can be applied to systems of interacting particles with different species including reactions among them [31]. We did not treat the case where the parameters taking place in the evolution equation depends on time but the present approach can also be applied for instance to the case where the temperature oscillate periodically in time [39, 40]. In this case, for long times the system may not properly reaches a stationary state in the sense of being independent of time but may reach a state with a probability that oscillates in time.

Stochastic thermodynamics is sometimes called a discipline whose quantities are defined at the level of single trajectories. This denomination emphasizes the fluctuation aspect of the theory, which is a relevant feature in systems with few degrees of freedom, the main application of the theory. In this respect the theory looks like statistical mechanics, which incorporates fluctuations and may also be applied to small systems. Thus, an alternative name to discipline would be stochastic mechanics avoiding the term thermodynamics which is usually associated to macroscopic systems.

The emphasis on trajectories and path integralsl is a distinguish feature as is used for instance in the Jarzynski equality (162) The present approach, on the other hand, the emphasis rests on the fluxes and currents of various types but a relationship between fluxes and path integrals exists as shown in section 15, revealing the equivalence between the two approaches. The present approach also emphasizes the connection with the laws of thermodynamics, particularly the second law expressed by the nonnegativity of the rate of entropy production.

The present approach to the quantum stochastic thermodynamics is based on the quantum evolution equation which is a canonical quantization of the FPK equation. It differs from other approaches such as those based on the Lindblad operators [41]. However, the present quantum FPK equation has similarity with the quantum master equation derived by Dekker [42] and by Caldeira and Leggett [43, 44]. The main features of the quantum FPK equation that distinguishes from other approaches is that it is centered on the irreversible density current operator, the analog of the classical irreversible probability current, which plays a fundamental role in defining the fluxes of various type as well as the rate of entropy production. The similarity of the quantum evolution equation to the classical counterpart is useful because the generalization of the concepts of classical stochastic thermodynamics, such as those associated to the current of probability, to the quantum case become easier. A final distinguishing feature is that a system described by the quantum FPK equation thermalizes properly. That is, for long times the system approaches the equilibrium state, if, of course, the relationship (99) between noise and fluctuation is obeyed.

Appendix A

The FPK equation (2) can be written in the form

(167) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-\frac{\mathrm{\partial }{J}_x}{\mathrm{\partial }x}-\frac{\mathrm{\partial }{J}_p}{\mathrm{\partial }p},$$

where $J_x$ and $J_p$ are the components of the probability current, and given by

(168) $$J_x=\frac{p}{m},J_p=F\rho -\frac{\mathrm{\Gamma }}{2}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p},$$

Let us integrate both sides of the equation (167) in a region R of the space $(x,p)$ delimited by a boundary line,

(169) $$\frac{d}{dt}{\int }_R\rho dxdp=-{\int }_R\frac{\mathrm{\partial }{J}_x}{\mathrm{\partial }x}dxdp-{\int }_R\frac{\mathrm{\partial }{J}_p}{\mathrm{\partial }p}dxdp.$$

The first integral can be integrated in x,

(170) $${\int }_R\frac{\mathrm{\partial }{J}_x}{\mathrm{\partial }x}dxdp=\int [J_x(x_2,p)-J_x(x_1,p)]dp.$$

For simplicity we are considering that R is a convex region so that there are two values of x at the boundary of R for a given p, which we are denoting by $x_1(p)$ and $x_2(p)$. In an analogous way we write the second integral as

(171) $${\int }_R\frac{\mathrm{\partial }{J}_p}{\mathrm{\partial }p}dxdp=\int [J_p(x,p_1)-J_p(x,p_2)]dx.$$

If $J_x$ and $J_p$ vanish at the boundary then both integrals vanish and

(172) $$\frac{d}{dt}{\int }_R\rho dxdp=0,$$

from which follows that the integral is a constant that we set equal to unity,

(173) $${\int }_R\rho dxdp=1.$$

This result is extended to the case where the region R is the whole space of states, in which case we demand that $J_x$ and $J_p$ vanish at infinity, a requirement that is provided by demanding that $\rho$ vanishes rapidly at infinity.

Let us consider now the case of an integral of the type

(174) $$\int A\frac{\mathrm{\partial }B}{\mathrm{\partial }x}dxdp,$$

where the integral is over the whole space of states. If we perform and integration by parts the result is

(175) $$\int \frac{\mathrm{\partial }AB}{\mathrm{\partial }x}dxdp-\int \frac{\mathrm{\partial }A}{\mathrm{\partial }x}Bdxdp.$$

Assuming that $AB$ vanishes rapidly at the limits of integration as we did above, the first integral vanishes and we are left with the result

(176) $$\int A\frac{\mathrm{\partial }B}{\mathrm{\partial }x}dxdp=-\int \frac{\mathrm{\partial }A}{\mathrm{\partial }x}Bdxdp.$$

Appendix B

Here we solve the FPK equation for the case of a harmonic force $F^c=-Kx$. The equation is

(177) $$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-\frac{p}{m}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }x}+Kx\frac{\mathrm{\partial }\rho }{\mathrm{\partial }p}+\gamma \frac{\mathrm{\partial }p\rho }{\mathrm{\partial }p}+\frac{m\gamma }{\beta }\frac{{\mathrm{\partial }}^2\rho }{\mathrm{\partial }{p}^2},$$

and it can be solved exactly by assuming the following Gaussian form for the probability distribution

(178) $$\rho =\frac{1}{Z}{e}^{-(a{x}^2+b{p}^2+2cxp)/2},$$

where the parameters a, b, c and Z depend on time. Replacing this form in the equation (177), we see that the left and right-hand sides will only have terms of the types $x^2$, $p^2$ and $xp$. Equating the respective coefficients of these terms we find equations for the parameters a, b, and c. There is no need to seek an equation for Z because this quantity can be obtained from the three parameters a, b, and c. This follows from the normalization of (178), which gives

(179) $$Z=\int e^{-(a{x}^2+b{p}^2+2cxp)/2}dxdp=\frac{2\pi }{\sqrt{ab-c^2}}.$$

Replacing the Gaussian distribution (178) in the FPK equation we may find the equations for the three parameters. However, the equations are too complicated and we will instead seek for equations that determine the covariances $⟨x^2⟩$, $⟨p^2⟩$, and $⟨xp⟩$. Before that we should write down the relations between the covariances and the three parameters, which are are obtained from the probability distribution (178), and are

(180) $$⟨x^2⟩=\frac{b}{ab-c^2},$$

(181) $$⟨p^2⟩=\frac{a}{ab-c^2},$$

(182) $$⟨xp⟩=\frac{-c}{ab-c^2}.$$

Inverting these relations we find

(183) $$a=\frac{⟨p^2⟩}{⟨x^2⟩⟨p^2⟩-⟨xp{⟩}^2},$$

(184) $$b=\frac{⟨x^2⟩}{⟨x^2⟩⟨p^2⟩-⟨xp{⟩}^2},$$

(185) $$c=\frac{-⟨xp⟩}{⟨x^2⟩⟨p^2⟩-⟨xp{⟩}^2}.$$

It remains now to determine the covariances as functions of time. To find the equations for the covariances we proceed as follows. We multiply both sides of the FPK equation successively by $x^2$, $p^2$ and $xp$, and integrate in x and p. Performing appropriate integration by parts, we find

(186) $$\frac{d}{dt}⟨x^2⟩=\frac{2}{m}⟨xp⟩,$$

(187) $$\frac{d}{dt}⟨p^2⟩=-2K⟨xp⟩-2\gamma ⟨p^2⟩+\frac{2m\gamma }{\beta },$$

(188) $$\frac{d}{dt}⟨xp⟩=\frac{1}{m}⟨p^2⟩-K⟨x^2⟩-\gamma ⟨xp⟩.$$

The stationary solution of this equation is $⟨x^2⟩=1/K\beta$, $⟨p^2⟩=m/\beta$, and $⟨xp⟩=0$. Taking these result into account, we define variables that are deviations of the covariances from their stationary values as follows, $A=⟨x^2⟩-1/K\beta$, $B=⟨p^2⟩-m/\beta$, and $C=⟨rp⟩$. These variables obey the set of linear equations

(189) $$\frac{dA}{dt}=\frac{2}{m}C,$$

(190) $$\frac{dB}{dt}=-2KC-2\gamma B,$$

(191) $$\frac{dC}{dt}=\frac{1}{m}B-KA-\gamma C,$$

which we write in matrix form

(192) $$\frac{d}{dt}\left(\begin{array}{c}A\\ B\\ C\end{array}\right)=\left(\begin{array}{rrr}0& 0& 2/m\\ 0& -2\gamma & -2K\\ -K& 1/m& -\gamma \end{array}\right)\left(\begin{array}{c}A\\ B\\ C\end{array}\right).$$

The solution for each variable is of the type $e^{\lambda t}$ where $\lambda$ is an eigenvalue of the square matrix above. They are

(193) $${\lambda }_1=-\gamma +\sqrt{{\gamma }^2-4K/m},$$

(194) $${\lambda }_2=-\gamma ,$$

(195) $${\lambda }_3=-\gamma -\sqrt{{\gamma }^2-4K/m},$$

and are all negative. The general solution is

(196) $$A=A_1{e}^{{\lambda }_{1}t}+A_2{e}^{{\lambda }_{2}t}+A_3{e}^{{\lambda }_{3}t},$$

(197) $$B=B_1{e}^{{\lambda }_{1}t}+B_2{e}^{{\lambda }_{2}t}+B_3{e}^{{\lambda }_{3}t},$$

(198) $$C=C_1{e}^{{\lambda }_{1}t}+C_2{e}^{{\lambda }_{2}t}+C_3{e}^{{\lambda }_{3}t},$$

and the coefficients are not all independent, but are related by

(199) $${\lambda }_i{A}_i=\frac{2}{m}{C}_i,$$

(200) $${\lambda }_i{B}_i=-2K{C}_i-2\gamma B_i,$$

(201) $${\lambda }_i{C}_i=\frac{1}{m}{B}_i-K{A}_i-\gamma C_i.$$

Thus only three, say $A_1$, $A_2$, and $A_3$ can be chosen to be independent, and they are determined by the initial conditions.

It is worth determining the solution for long times. In this case the solution is dominated by the largest eigenvalue, which is ${\lambda }_1$. The covariances are

(202) $$⟨x^2⟩=\frac{1}{K\beta }+A_1{e}^{{\lambda }_{1}t},$$

(203) $$⟨p^2⟩=\frac{m}{\beta }+B_1{e}^{{\lambda }_{1}t},$$

(204) $$⟨xp⟩=C_1{e}^{{\lambda }_{1}t},$$

from which we find the three parameters

(205) $$a=K\beta -a_1{e}^{{\lambda }_{1}t},$$

(206) $$b=\frac{\beta }{m}-b_1{e}^{{\lambda }_{1}t},$$

(207) $$c=-c_1{e}^{{\lambda }_{1}t},$$

where $a_1=K^2{\beta }^2{A}_1$, $b_1={\beta }^2{B}_1/m^2$, $c_1=K{\beta }^2{C}_1/m$.

Appendix C

We wish to determine here in an explicit form the dissipative force $D=-\gamma$ for the quantum harmonic oscillator were

(208) $$g=-\frac{m}{i\hslash \beta }(e^{\beta H}x{e}^{-\beta H}-x),$$

and

(209) $$H=\frac{1}{2m}{p}^2+\frac{1}{2}m{\omega }^2{x}^2,$$

To this end we start with the following identity [38]

(210) $$\begin{gathered} e^{\beta H}x{e}^{-\beta H}=x+\beta [H,x]+\frac{{\beta }^2}{2}[H,[H,x]]+ \\ +\frac{{\beta }^3}{3!}[H,[H,[H,x]]]+\dots \end{gathered}$$

Using the notation

(211) $$A_n=[H,[H,\dots [H,x]\dots ]],$$

where the numbers of commutator is equal to n, the identity above is written as

(212) $$e^{\beta H}x{e}^{-\beta H}=\sum _{n=0}^{\mathrm{\infty }}\frac{{\beta }^n}{n!}{A}_n,$$

where $A_0=x$. The quantities $A_n$ obeys the recursive relations

(213) $$A_{n+1}=[H,A_n].$$

To determine $A_n$ is easier if we use the relations

(214) $$[H,x]=-\frac{i\hslash }{m}p,$$

(215) $$[H,p]=im{\omega }^2\hslash x,$$

which are obtained by using the commutation relation $[x,p]=i\hslash$. From these relations we get the two useful rules,

(216) $$[H,[H,x]={\hslash }^2{\omega }^{2}x,$$

(217) $$[H,[H,p]={\hslash }^2{\omega }^{2}p.$$

The first two coefficients of the expansion are

(218) $$A_0=x,A_1=-\frac{i\hslash }{m}p.$$

Next, with the two rules above in mind, we observe that $A_2$ will be proportional to x and $A_3$ will be proportional to p, and, in general, $A_n$ will be proportional to x if n is even, and it will be proportional to p if n is odd.

Let us consider the case n even and write $A_n=a_{n}x$. Then using the two rules above,

(219) $$A_{n+2}=[H,[H,A_n]]={\hslash }^2{\omega }^2{a}_{n}x,$$

so that

(220) $$a_{n+2}={\hslash }^2{\omega }^2{a}_n,$$

from which we find

(221) $$a_n=(\hslash \omega {)}^n,$$

because $a_0=1$. The part of the expansion (212) corresponding to n even is

(222) $$\sum _{n(\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n})}\frac{{\beta }^n}{n!}{A}_n=x\sum _{n(\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n})}\frac{(\beta \hslash \omega {)}^n}{n!}=x\cosh \beta \hslash \omega .$$

Now we consider the case n odd and write $A_n=b_{n}p$. Then using the two rules above

(223) $$A_{n+2}=[H,[H,A_n]]={\hslash }^2{\omega }^2{b}_{n}p,$$

so that

(224) $$b_{n+2}={\hslash }^2{\omega }^2{b}_n,$$

from which we find

(225) $$b_n=-\frac{i}{m\omega }(\hslash \omega {)}^n,$$

because $b_1=-i\hslash /m$. The part of the expansion (212) corresponding to n odd is

(226) $$\begin{array}{rcl}& & \sum _{n(\mathrm{o}\mathrm{d}\mathrm{d})}\frac{{\beta }^n}{n!}{A}_n=-\frac{ip}{m\omega }\sum _{n(\mathrm{o}\mathrm{d}\mathrm{d})}\frac{(\beta \hslash \omega {)}^n}{n!}\\ & & =-\frac{ip}{m\omega }\sinh \beta \hslash \omega .\end{array}$$

Collecting the results above we find

(227) $$e^{\beta H}x{e}^{-\beta H}=-\frac{ip}{m\omega }\sinh \beta \hslash \omega +x\cosh \beta \hslash \omega ,$$

and the quantity (208) becomes

(228) $$g=p\frac{1}{\beta \hslash \omega }\sinh \beta \hslash \omega +ix\frac{m}{\beta \hslash }(\cosh \beta \hslash \omega -1).$$

Appendix D

Let suppose that f is a function of several variables that we denote by $x_i$, and that these variables depend on time. The derivative of f with respect to time is

(229) $$\frac{df}{dt}=\sum _i{f}_i\frac{d{x}_i}{dt},$$

where $f_i$ are functions of $x_i$ given by

(230) $$f_i=\frac{\mathrm{\partial }f}{\mathrm{\partial }{x}_i}.$$

Equation (229) can be written in simplified form

(231) $$df=\sum _i{f}_{i}d{x}_i,$$

where $d{x}_i$ are the differentials of the variable $x_i$ and $df$ is the differential of f. Since f is a function of $x_i$ then the following relation is valid

(232) $$\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }{f}_j}{\mathrm{\partial }{x}_i}.$$

Now we raise the following question. Let g be a function of t and let $x_i$ depend on time as before and let us assume that

(233) $$\frac{dg}{dt}=\sum _i{g}_i\frac{d{x}_i}{dt},$$

where $g_i$ are given function of the variables $x_j$. The question now arises whether g could depend on time only through the variables $x_i$, that is, whether

(234) $$g(t)=g(x_1(t),x_2(t),\dots ).$$

If that is possible then according to our reasoning above the given functions $g_i$ of the variables $x_j$ must fulfill the condition

(235) $$\frac{\mathrm{\partial }{g}_i}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }{g}_j}{\mathrm{\partial }{x}_i},$$

for all pairs $i,j$. If this condition is not satisfied it is not possible to write g as a function of $x_i$. In this case if we write (233) in the simplified form

(236) $$dg=\sum _i{g}_{i}d{x}_i.$$

we say that $dg$ is not an exact differential.

Appendix E

We determine here the covariances $⟨x_i{x}_j$, $⟨x_i{p}_j⟩$ and $⟨p_i{p}_j⟩$ for the system described by the FPK equation (118) where $F_1^c=-K{x}_1+L{x}_2$, $F_2^c=-K{x}_2+L{x}_1$, and $F_1^e=-c{x}_2$ and $F_2^e=c{x}_1$, which we reproduce here in the following form

$$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=-\frac{p_1}{m}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }{x}_1}-\frac{p_2}{m}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }{x}_2}$$

$$+\frac{\mathrm{\partial }}{\mathrm{\partial }{p}_1}(K{x}_1+b{x}_2+\gamma p_1)\rho +\frac{\mathrm{\partial }}{\mathrm{\partial }{p}_2}(K{x}_2+a{x}_1+\gamma p_2)\rho$$

(237) $$+\gamma mk{T}_1\frac{{\mathrm{\partial }}^2\rho }{\mathrm{\partial }{p}_1^2}+\gamma mk{T}_2\frac{{\mathrm{\partial }}^2\rho }{\mathrm{\partial }{p}_2^2},$$

where $b=-L+c$ and $a=-L-c$.

Multiplying successively equation (237) by $x_i{x}_j$, $x_i{p}_j$, and $p_i{p}_j$, and performing the integration we find the following equations, after appropriate integration by parts,

(238) $$\frac{d}{dt}⟨x_1^2⟩=\frac{2}{m}⟨x_1{p}_1⟩,$$

(239) $$\frac{d}{dt}⟨x_2^2⟩=\frac{2}{m}⟨x_2{p}_2⟩,$$

(240) $$\frac{d}{dt}⟨x_1{x}_2⟩=\frac{1}{m}⟨x_1{p}_2⟩+\frac{1}{m}⟨x_2{p}_1⟩,$$

(241) $$\begin{array}{rcl}& & \frac{d}{dt}⟨x_1{p}_2⟩=\frac{1}{m}⟨p_1{p}_2⟩-K⟨x_1{x}_2⟩\\ & & -a⟨x_1^2⟩-\gamma ⟨x_1{p}_2⟩,\end{array}$$

(242) $$\frac{d}{dt}⟨x_2{p}_1⟩=\frac{1}{m}⟨p_1{p}_2⟩-K⟨x_1{x}_2⟩-b⟨x_2^2⟩-\gamma ⟨x_2{p}_1⟩,$$

(243) $$\frac{d}{dt}⟨x_1{p}_1⟩=\frac{1}{m}⟨p_1^2⟩-K⟨x_1^2⟩-b⟨x_1{x}_2⟩-\gamma ⟨x_1{p}_1⟩,$$

(244) $$\frac{d}{dt}⟨x_2{p}_2⟩=\frac{1}{m}⟨p_2^2⟩-K⟨x_2^2⟩-a⟨x_1{x}_2⟩-\gamma ⟨x_2{p}_2⟩,$$

(245) $$\begin{array}{rcl}& & \frac{d}{dt}⟨p_1^2⟩=-2K⟨x_1{p}_1⟩-2b⟨x_2{p}_1⟩\\ & & -2\gamma ⟨p_1^2⟩+2\gamma mk{T}_1,\end{array}$$

(246) $$\begin{array}{rcl}& & \frac{d}{dt}⟨p_2^2⟩=-2K⟨x_2{p}_2⟩-2a⟨x_1{p}_2⟩\\ & & -2\gamma ⟨p_2^2⟩+2\gamma mk{T}_2,\end{array}$$

$$\frac{d}{dt}⟨p_1{p}_2⟩=-K⟨x_1{p}_2⟩-K⟨x_2{p}_1⟩$$

(247) $$-b⟨x_2{p}_2⟩-a⟨x_1{p}_1⟩-2\gamma ⟨p_1{p}_2⟩.$$

Now we look for the stationary solution. Setting the above equation to zero, we find that the following covariances vanish, $⟨x_1{p}_1⟩=0$, $⟨x_2{p}_2⟩=0$, $⟨p_1{p}_2⟩=0$. The other covariances are the solution of the set of linear equations

(248) $$⟨x_1{p}_2⟩+⟨x_2{p}_1⟩=0,$$

(249) $$K⟨x_1{x}_2⟩+a⟨x_1^2⟩+\gamma ⟨x_1{p}_2⟩=0,$$

(250) $$K⟨x_1{x}_2⟩+b⟨x_2^2⟩+\gamma ⟨x_2{p}_1⟩=0,$$

(251) $$⟨p_1^2⟩-mK⟨x_1^2⟩-bm⟨x_1{x}_2⟩=0,$$

(252) $$⟨p_2^2⟩-mK⟨x_2^2⟩-am⟨x_1{x}_2⟩=0,$$

(253) $$b⟨x_2{p}_1⟩+\gamma ⟨p_1^2⟩=\gamma mk{T}_1,$$

(254) $$a⟨x_1{p}_2⟩+\gamma ⟨p_2^2⟩=\gamma mk{T}_2.$$

A straightforward calculation leads us to the result

(255) $$⟨x_1{p}_2⟩=-⟨x_2{p}_1⟩=\frac{m\gamma k(b{T}_2-a{T}_1)}{2(m{\gamma }^{2}K+ab)},$$

(256) $$⟨x_1{x}_2⟩=-\frac{k(a{T}_1+b{T}_2)}{2(K^2-ab)},$$

(257) $$⟨p_1^2⟩=mk{T}_1+\frac{b}{\gamma }⟨x_1{p}_2⟩,$$

(258) $$⟨p_2^2⟩=mk{T}_2-\frac{a}{\gamma }⟨x_1{p}_2⟩,$$

(259) $$⟨x_1^2⟩=-\frac{\gamma }{a}⟨x_1{p}_2⟩-\frac{K}{a}⟨x_1{x}_2⟩,$$

(260) $$⟨x_2^2⟩=\frac{\gamma }{b}⟨x_1{p}_2⟩-\frac{K}{b}⟨x_1{x}_2⟩.$$

We remark that, as $⟨x_1^2⟩$, $⟨x_2^2⟩$, $⟨p_1^2⟩$, and $⟨p_2^2⟩$ must be nonnegative, the following conditions should be fulfilled

(261) $$mK{\gamma }^2+ab\ge 0,K^2-ab\ge 0.$$

It is worth mentioning that the probability density can also be determined. On account of the linearity of the FPK equation in relation to the variable $x_i$ and $p_i$, the solution is a multivariate Gaussian distribution, which we write as

(262) $$\rho =\frac{1}{Z}\exp \{-\frac{1}{2}\sum _{i,j=1}^4{L}_{ij}{\xi }_i{\xi }_j\},$$

where we are using the abbreviations ${\xi }_1=x_1$, ${\xi }_2=x_2$, ${\xi }_3=p_1$, and ${\xi }_4=p_2$. The matrix L whose elements are $L_{ij}$ is the inverse of the covariance matrix C, whose elements we have just determined, and are $C_{11}=⟨x_1^2⟩$, $C_{22}=⟨x_2^2⟩$, $C_{12}=⟨x_1{x}_2⟩$, $C_{33}=⟨p_1^2⟩$, $C_{44}=⟨p_2^2⟩$, $C_{34}=⟨p_1{p}_2⟩$, $C_{13}=⟨x_1{p}_1⟩$, $C_{14}=⟨x_1{p}_2⟩$, $C_{23}=⟨x_2{p}_1⟩$, and $C_{24}=⟨x_2{p}_2⟩$. The expression given by equation (262) is the probability distribution describing the nonequilibrium stationary state of the present problem.

References

• [1] J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976).
• [2] L. Jiu-Li, C. Van den Broeck and G. Nicolis, Z. Phys. B 56,165 (1984).
• [3] C.Y. Mou, J.L. Luo and G. Nicolis, J. Chem. Phys. 84, 7011 (1986).
• [4] A. Pérez-Madrid, J.R. Rubí and P. Mazur, Physica A 212, 231 (1994).
• [5] T. Tomé and M.J. de Oliveira, Braz. J. Phys. 27, 525 (1997).
• [6] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
• [7] K. Sekimoto, Prog. Theor. Phys. Suppl. 130, 17 (1998).
• [8] J.L. Lebowitz and H. Spohn, J. Stat. Phys. 95, 333 (1999).
• [9] P. Mazur, Physica A 274, 491 (1999).
• [10] C. Maes and K. Netočný, J. Stat. Phys. 110, 269 (2003).
• [11] L. Crochik and T. Tomé, Phys. Rev. E 72, 057103 (2005).
• [12] T. Tomé, Braz. J. Phys. 36, 1285 (2006).
• [13] R.K.P. Zia and B Schmittmann, J. Phys. A: Math. Gen. 39, L407 (2006).
• [14] D. Andrieux and P. Gaspard, Phys. Rev. E 74, 011906 (2006).
• [15] T. Schmiedl and U. Seifert, J. Chem. Phys. 126, 044101 (2007).
• [16] R.J. Harris and G. M. Schütz, J. Stat. Mech. 2007,P07020 (2007).
• [17] U. Seifert, Eur. Phys. J. B 64, 423 (2008).
• [18] R.A. Blythe, Phys. Rev. Lett. 100, 1010060 (2008).
• [19] M. Esposito, K. Lindenberg and C. Van den Broeck, Phys. Rev. Lett. 102, 130602 (2009).
• [20] M. Esposito, U. Harbola and S. Mukamel, Rev. Mod. Phys. 81, 1665 (2009).
• [21] T. Tomé and M.J. de Oliveira, Phys. Rev. E 82, 021120 (2010).
• [22] C. Van de Broeck and M. Esposito, Phys. Rev. E 82, 011144 (2010).
• [23] C. Jarzynski, Annual Review of Condensed Matter Physics 2, 329 (2011).
• [24] T. Tomé and M.J. de Oliveira, Phys. Rev. Lett. 108, 020601 (2012).
• [25] R.E. Spinney and I.J. Ford, Phys. Rev. E 85, 051113 (2012).
• [26] U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
• [27] M. Santillan and H. Qian, Physica A 392, 123 (2013).
• [28] D. Luposchainsky and H. Hinrichsen, J. Stat. Phys. 153, 828 (2013).
• [29] T. Tomé and M.J. de Oliveira, Phys. Rev. E 91, 042140 (2015).
• [30] M.J. de Oliveira, Phys. Rev. E 94, 012128 (2016).
• [31] T. Tomé and M.J. de Oliveira, J. Chem. Phys. 148, 224104 (2018).
• [32] M.J. de Oliveira, Phys. Rev. E 99, 052138 (2019).
• [33] N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1981).
• [34] C.W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and Natural Sciences (Springer, Berlin, 1983).
• [35] H. Risken, The Fokker-Planck Equation, Methods of Solution and Applications (Springer, Berlin, 1984).
• [36] T. Tomé and M.J. de Oliveira, Stochastic Dynamics and Irreversibility (Springer, Heidelberg, 2015).
• [37] S.R.A. Salinas, Introduction to Statistical Physics (Springer, New York, 2001).
• [38] E. Merzbacher, Quantum Mechanics (Wiley, New York, 1970). 2ª ed.
• [39] M.J. Oliveira, J. Stat. Mech., 073204 (2019).
• [40] C.E. Fiore, M.J. de Oliveira, Phys. Rev. E 99, 052131 (2019).
• [41] G. Lindblad, Comm. Math. Phys. 48, 19 (1976).
• [42] H. Dekker, Phys. Rev. A 16, 2116 (1977).
• [43] A.O. Caldeira and A. Leggett, Physica 121, 587 (1983).
• [44] A.O. Caldeira, Introduction to Macroscopic Quantum Phenomena and Quantum Dissipation (Cambridge University Press, Cambridge, 2014).

Publication Dates

• Publication in this collection
  14 Sept 2020
• Date of issue
  2020

History

• Received
  21 May 2020
• rev-received
  11 June 2020
• Accepted
  18 June 2020