Kinetic Energy of a Free Quantum Brownian Particle

We consider a paradigmatic model of a quantum Brownian particle coupled to a thermostat consisting of harmonic oscillators. In the framework of a generalized Langevin equation, the memory (damping) kernel is assumed to be in the form of exponentially-decaying oscillations. We discuss a quantum counterpart of the equipartition energy theorem for a free Brownian particle in a thermal equilibrium state. We conclude that the average kinetic energy of the Brownian particle is equal to thermally-averaged kinetic energy per one degree of freedom of oscillators of the environment, additionally averaged over all possible oscillators’ frequencies distributed according to some probability density in which details of the particle-environment interaction are present via the parameters of the damping kernel.


Introduction
One of the enduring milestones of classical statistical physics is the theorem of the equipartition of energy [1], which states that energy is shared equally amongst all energetically-accessible degrees of freedom of a system and relates average energy to the temperature T of the system. In particular, for each degree of freedom, the average kinetic energy is equal to E k = k B T/2, where k B is the Boltzmann constant. This relation is exploited in various aspects of many areas of physics, chemistry and biology. However, in many cases, it is applied in an unjustified way forgetting about assumptions used in proving this theorem. One can notice confusion and mess, in particular in the case of quantum systems. In the standard course of classical statistical physics, the equality E k = k B T/2 is derived under the following conditions: 1.
The system is at thermal equilibrium of temperature T.

2.
The state of the system is described by the Gibbs canonical probability distribution. 3.
The Gibbs probability distribution describes an equilibrium state of the system in the limit of weak coupling with the thermostat. 4.
The Gibbs probability distribution does not depend on the system-thermostat coupling constant.
It should be stressed that in classical statistical physics, the equality E k = k B T/2 is universal: it does not depend either on the number of particles of the system or on the single-particle potential U( r i ) in which the i-th particle dwells, as well as it does not depend on the form of mutual interaction U( r i , r j ) between i-th and j-th particles of the system.
In quantum physics, the problem is more complicated. In the weak coupling limit, the density operator (or the density matrix) ρ for the canonical ensemble describes a thermal equilibrium state and has the form: where H is a Hamiltonian of the system. In the book by Feynman [2], one can find Expression (2.88) for average kinetic energy E k (ω 0 ) of a quantum harmonic oscillator of the eigenfrequency ω 0 . It has the form: where p is momentum and m is the mass of the oscillator. The Hamiltonian of the oscillator has the well-known form: We can perform the limit ω 0 → 0, which corresponds to the Hamiltonian of a free particle. In this limit, Equation (2) assumes the form: It is the same expression as for the classical free particle. However, Equations (2) and (4) are different. This means that for quantum systems, in contrast to a classical case, the average kinetic energy depends on the potential U(x), even in the weak coupling limit. If the weak coupling limit does not hold, the problem is even much more complicated. It is the aim of this paper to discuss the question of equipartition energy for an arbitrary system-thermostat coupling. We study only one specific and as simple as possible model of a quantum open system to present basic concepts and ideas. Therefore, we consider a free quantum particle coupled to its environment, which is modeled as a collection of harmonic oscillators of temperature T. This old clichéd system-environment model [3] has been re-considered many, many times by each next generation of physicists [4]. However, it is still difficult to find a transparent presentation of this fundamental issue of the quantum statistical physics focused on the kinetic energy. To achieve the aim, we try to use the simplest techniques and methods to make the paper consistent and self-contained, and we neglect many unnecessary and redundant aspects of the theory.
The remaining part of the paper is organized as follows. In Section 2, we describe the model and derive the integro-differential Generalized Langevin Equation (GLE); in Section 3, the fluctuation-dissipation relation is presented; the form of the dissipation integral kernel of GLE is described in Section 4; in Section 5, we convert the integro-differential Langevin equation into a set of differential equations; the equipartition theorem is discussed in Section 6; some selected physical regimes are analyzed in Section 7; we conclude the work with a brief résumé in Section 8; in Appendices A-C, we present some auxiliary calculations.

Hamiltonian Model and Generalized Langevin Equation
A paradigmatic model of a one-dimensional quantum Brownian motion consists of a particle of mass M subjected to the potential U(x) and interacting with a large number of independent oscillators, which form a thermostat (environment) E of temperature T being in an equilibrium canonical (Gibbs) state. The quantum-mechanical Hamiltonian of such a total system can be written in the form [3][4][5]: where the coordinate and momentum operators {x, p} refer to the Brownian particle and {q i , p i } are the coordinate and momentum operators of the i-th heat bath oscillator of mass m i and the eigenfrequency ω i . The parameter η i characterizes the interaction strength of the particle with the i-th oscillator. There is the counter-term, the last term proportional to x 2 , which is included to cancel a harmonic contribution to the particle potential. All coordinate and momentum operators obey canonical equal-time commutation relations. The next step is to write the Heisenberg equations of motion for all coordinate and momentum operators. For the Brownian particle, the Heisenberg equations are: where: For the environment operators, one gets: What we need in Equation (7) is the solution q i = q i (t), which can be obtained from Equations (9) and (10) with the result (see [6]): The following step is to integrate by parts the last term in Equation (11) and insert it into Equation (7) for p = p(t). Using Equation (6), after some algebra, one can obtain an effective equation of motion for the particle coordinate x(t). It is called a generalized Langevin equation and has the form: where U (x) denotes differentiation with respect to x, γ(t) is a dissipation function (damping or memory kernel) and F(t) denotes the random force, The dynamics of the quantum Brownian particle is therefore described by a stochastic integro-differential equation for the coordinate operator x(t). Probably Magalinskij [3] was the first, in 1959, to derive the integro-differential Equation (12) and formulated the problem in the above way. Next, from 1966, a series of papers were published on this topic, but a complete list of papers is too long to present here. We cite a part of them [7][8][9][10][11][12]. Generally, it is difficult to analyze this kind of integro-differential equation for operators. However, in the case of a free Brownian particle (when U(x) = 0) or for a harmonic oscillator (when U(x) ∝ x 2 ), Equation (12) can be solved exactly, at least in a formal way. In the classical case, equations like Equation (12) describe non-Markovian stochastic processes [5]. In the quantum case, there is no good definition of Markovian or non-Markovian processes, and therefore, a classification with respect to these notions is not constructive.

Fluctuation-Dissipation Theorem
In the standard approach, it is assumed that the initial state ρ(0) of the total system is uncorrelated, i.e., ρ(0) = ρ S ⊗ ρ T , where ρ S is an arbitrary state of the Brownian particle and ρ T is an equilibrium canonical state of the thermostat (environment) E of temperature T, namely, where: is the Hamiltonian of the thermostat. The operator-valued random force F(t) is a family of non-commuting operators whose commutators are c-numbers. Its mean value is zero, and the symmetrized correlation function: takes the form: The symmetrization is needed for C(t 1 , t 2 ) to be a real function with a correct limit in the classical case. We observe that it depends on the time difference C(t 1 , Statistical characteristics of the random force F(t) are similar to characteristics for a classical stationary Gaussian stochastic process, which models thermal equilibrium noise. Therefore, F(t) is called the Gaussian operator, which represents quantum thermal noise.
It is convenient to introduce the spectral function: Then, the damping kernel Equation (13) can be expressed as: and the correlation function Equation (19) reads: If we introduce the Fourier transforms of the dissipation and correlation functions, then we see that the following equality: holds. We obtain the relation between the dissipationγ(ω) and correlationĈ(ω) spectra. This is what is named the fluctuation-dissipation theorem [12][13][14]. In this relation, quantum effects are incorporated via the prefactor in r.h.s. of Equation (24).

Dissipation Function
The influence of the thermostat E on the Brownian particle is manifested through the correlation function Equation (19) of the thermostat. If it is a finite quantum system, then its energy spectrum is discrete, and all its correlation functions are almost periodic in time. On the other hand, if the thermostat is an infinitely-extended system, then its correlation functions decay. Let us observe that if E is finite, then the spectral function J(ω) is a completely singular distribution (the sum in Equation (20) is over discrete values of i), and if E is infinite, then J(ω) is, at least on some intervals, a continuous function of ω (in the thermodynamic limit for E, the summation in Equation (20) is replaced by an integral over a frequency). Below and up to the end of the paper, we assume that the environment E is an infinite system.
In order to start an analysis of GLE Equation (12), one has to specify at least one of three quantities: the dissipation kernel γ(t), or the correlation function C(t) of the random force F(t), or the spectral function J(ω). We assume γ(t) to have the form: with three non-negative parameters γ 0 , τ c and Ω. The parameter γ 0 is the system-environment coupling strength, and τ c defines decay or relaxation time and characterizes memory effects. Finally, Ω is the frequency in the relaxation process of the momentum (velocity). It is important to mention that τ c is, via the fluctuation-dissipation theorem, the correlation time of quantum thermal fluctuations and appears in the correlation function C(t 1 , t 2 ) defined by Equation (18). Why this form of the dissipation function? There are several reasons for this choice: -In the classical limit, it is a direct relation between the dissipation kernel γ(t) and the correlation function C(t) : C(t) = k B Tγ(t). Therefore, for Equation (25), the correlation function is sufficiently universal and has been considered for many systems. - When Ω = 0, it reduces to the exponential form and is known as a Drude model. Moreover, it has been considered in relation to colored noise problems. - When Ω = 0 and τ c → 0, then γ(t) tends to the Dirac delta function, and the integral term reduces to γ 0ẋ (t) (no memory effects), which is the well-known Stokes force with γ 0 interpreted as a friction (damping) coefficient. Here, the parameter γ 0 plays the role of coupling the strength of the Brownian particle to the thermostat. -This form allows converting the integro-differential Equation (12) into a set of differential equations, which can be handled by known methods. It is important from a technical point of view to have a method of any sort to analyze Equation (12).

Generalized Langevin Equation as a Set of Differential Equations
The integral part of Equation (12) is the convolution of γ(t) andẋ(t). It suggests applying integral transforms like the Laplace or Fourier ones to solve it. Here, we exploit another method, which is based on the observation that if γ(t) fulfills a linear ordinary differential equation with constant coefficients, then Equation (12) can be converted to a set of ordinary differential equations. Note that the function γ(t) in Form (25) fulfills a differential equation of second order, which is similar to the Newton equation for a damped harmonic oscillator. We introduce auxiliary variables (in fact, operators) u(t) and v(t) by the relations: Then, Equation (12) is converted to the following set of differential equations: To proceed further, we have to specify the form of the potential U(x). The simplest case is when U(x) = 0, i.e., for the free particle. The second exactly solvable case is the harmonic oscillator with U(x) = (1/2)Mω 2 0 x 2 . For other forms of U(x), the problem Equation (29) cannot be solved, and only mathematically uncontrolled approximations have been applied. The main elements of analysis of averaged kinetic energy are similar for both a free particle and a harmonic oscillator. However, a more pedagogical example is a free particle because the calculations are less tedious. It is the only reason for our choice U(x) = 0. In such a case, in order to calculate the averaged kinetic energy, it is sufficient to consider the reduced set of equations: It can be rewritten in the matrix form: where: and T denotes the transpose of a matrix, which switches the row into the column. The matrix A has the form: The solution of the the non-homogeneous linear differential Equation (31) reads [6]: where: The spectrum of the matrix A and its invariant subspaces determine the time dependence of Equation (35). Now, the only problem is to determine the exponential of the matrix At, i.e., the matrix R(t), which can be computed in many ways. The authors of the paper [15] say about 19 ways. As they write: "In practice, consideration of computational stability and efficiency indicates that some of the methods are preferable to others, but that none are completely satisfactory". The traditional way is to transform A into its Jordan canonical form. Here, we will use a less traditional method, namely the Putzer algorithm [16], in which the exponential of the matrix At can be computed knowing nothing more than the eigenvalues of the matrix A. Moreover, the algorithm does not require that the matrix A is diagonalizable. We think that this method is simple, elegant and suitable for presentation to students and younger researchers. It is described in Appendix A.

Average Kinetic Energy in Equilibrium
The operator of the kinetic energy H k (t) = p 2 (t)/2M is expressed by the momentum p(t), which is the first component of the vector X(t) determined by Equation (35). We calculate its average in the long time limit t → ∞ when a stationary state is approached. This stationary state is a thermal equilibrium state. The first component of X(t) is: where R 11 (t) is the first element of the matrix R(t). As is shown in Appendix A, elements of this matrix are exponentially decreasing functions of time. This means that the average value of the momentum p(t) → 0 as t → ∞. To evaluate the average kinetic energy, we consider the symmetrized momentum-momentum correlation function [p(t); p(u)] + . In the long time limit, the first two terms of Equation (37) do not contribute to it, and only the last term contributes, yielding: Now, we use the results of Section 3 and insert the expression for the correlation function Equation (23) of quantum thermal noise F(t) into Equation (38). The result is: In particular, for t = s, it is the second statistical moment of the momentum, We introduce new integration variables τ = t − t 1 and u = t − t 2 to convert Equation (40) into the form: Finally, in the long time limit t → ∞, average kinetic energy is obtained as: where: At this point, we use the fluctuation-dissipation relation Equation (24) to express the noise correlation spectrumĈ(ω) by the noise dissipation spectrumγ(ω), which, via Form (25) for the dissipation function γ(t) and its inverse Fourier transform, is given by: The function I(ω) in Equation (43) is calculated from the relation Equation (A11) in Appendix A. Its explicit form is given by Equation (A16) in Appendix B. The numerator of I(ω) cancels with the denominator in Equation (44), and finally, we obtain: where: We used the definition Equation (28) for µ = µ 0 ε. We rewrite it in this form because both µ 0 and ε have the unit of frequency (1/s). Moreover, µ 0 = γ 0 /M defines the rescaled coupling strength of the Brownian particle to the thermostat. The reciprocal 1/µ 0 = M/γ 0 has the unit of time, and in the case of a classical free Brownian particle, it is the relaxation time of the particle velocity v =ẋ (obtained from the Newton equation Mv = −γ 0 v) or, equivalently, it is the correlation time occurring in the velocity-velocity correlation function.
The function P(ω) has no poles on a real axis because the denominator in Equation (46) can be rewritten in the form x[(x + c 1 ) 2 + c 2 ] + c 3 > 0 for x = ω 2 , c 1 = ε 2 − Ω 2 − µ 0 ε, positive c 2 = 4ε 2 Ω 2 and positive c 3 = µ 2 0 ε 4 . As a consequence, the integral in Equation (45) exists for any fixed values of parameters. Various forms of the expression Equation (45) have previously been derived [17][18][19]. However, the influence of the dissipation function Equation (25) on system properties, in particular in the context of average kinetic energy, has not been discussed.
The expression Equation (45) with the integrand Equation (46) is a quantum version of the equipartition theorem. It differs from its classical counterpart, like in Equation (4). The dependence of E k on temperature is involved in the integrand of Equation (45) and cannot be presented in an explicit, simple way. We note that there are four characteristic times τ 1 = 1/µ 0 = M/γ 0 , τ 2 = τ c = 1/ε, τ 3 = 1/Ω and τ 4 =h/k B T (called the thermal time) and four corresponding frequencies µ 0 , ε, Ω and k B T/h (called the Matsubara frequency), which all influence the temperature dependence of E k . Now, we want to take a look at the relation Equation (45) from another point of view. To this aim, we use the expression Equation (2) for the average kinetic energy of a single harmonic oscillator and rewrite Equation (45) in the form: The function P(ω) is positive, P(ω) > 0, and normalizable (see Equation (A23) in Appendix C), Therefore, there exists a random variable ξ for which P is its probability distribution. From the mathematical point of view, Equation (47) is an average value of the function E k (ξ) of the random variable ξ (physicists frequently equate it with the integration variable). It allows presenting an interesting interpretation of the quantum equipartition theorem: the average kinetic energy E k of the Brownian particle is strongly related to the thermally-averaged kinetic energy E k (ω) per one degree of freedom of oscillators of the environment. Because frequencies ω of oscillators are random variables, E k (ω) has to be additionally averaged over all possible frequencies ω distributed according to the probability density P(ω) in which details of the particle-environment interaction are present via the parameters of the dissipation function γ(t).

Average Kinetic Energy in Terms of Series
From the relation Equation (45), it is difficult to draw conclusions on the dependence of the average kinetic energy on the system parameters. Therefore, we present another form of E k . To this aim, we exploit the series expansion [20]: x coth which allows calculating the integral in Equation (45). More details are presented in Appendix C. From Equation (A25) in this appendix, we get the expression: Here, the average kinetic energy is represented by an infinite series, and some information on E k can be inferred from this form. Since for n ≥ 1, all terms under the sum are non-negative, hence 2E k /k B T is a lower bound for the energy E k . Therefore, the kinetic energy in a quantum regime is always greater than the classical one. The term under the sum is a rational function of four characteristic energies k B T,hµ 0 ,hε,hΩ. The numerator and denominator are the products of energy to a power of three, like, e.g., (hµ 0 ) (hε) (k B T). It is easy to observe that each term under the sum is a non-increasing function with respect to Ω because it occurs only in the denominator. Moreover, it can be shown that partial derivatives of each term with respect to µ 0 and ε are non-negative, and it follows that all terms are non-decreasing with respect to µ 0 and ε, respectively. As a consequence, E k is a non-increasing function of Ω and a non-decreasing function of µ 0 and ε.
At a fixed temperature, the kinetic energy E k is finite for all finite values of parameters and behaves in the following way (see Figure 1): (i) E k increases monotonically to infinity when the coupling strength µ 0 increases to infinity. (ii) When the decay rate ε = 1/τ c grows from zero to infinity, E k grows from its classical value to infinity. In other words, for a very long decay time τ c , kinetic energy approaches its classical value Equation (4). When τ c → 0, the kinetic energy diverges. (iii) For Ω = 0 (the Drude model), E k starts from some maximal value for a given set of parameters, and next, it decreases when Ω increases.

High Temperature Regime
We focus now on the regime of high temperatures. In this regime, when T → ∞, we use the approximation: Then, Equation (45) reduces to the form: because the integral that occurs in Equation (45), according to Equation (48), is one. It is valid for any values of the parameters µ 0 , ε = 1/τ c and Ω, which characterize and are involved in the dissipation kernel γ(t) defined in Equation (25). In particular, Equation (52) holds for weak, as well as strong system-environment interactions. The same Form (52) is obtained from Equation (50). Indeed, each term under the sum behaves as T/T 3 ∼ 1/T 2 when T → ∞, and each term under the sum tends to zero. Only the first term survives, and Equation (50) is well approximated by Equation (52).
Let us now comment on the above approximations from the mathematical point of view. The approximation Equation (51) is valid whenhω/2k B T << 1. However, let us observe that for any large value of 2k B T, there is a larger value ofhω because the integration over ω in Equation (45) is to infinity and the inequalityhω/2k B T << 1 is not satisfied for all ω. Nevertheless, Equation (52) is a correct approximation of Equation (45), and it can be shown by rigorous mathematical means. (50) is uniformly convergent, we can take the limit T → ∞. Consequently, it implies Equation (52).

Low Temperature Regime
Let us note that the expression Equation (50) is not suitable for the calculation of the zero temperature limit because the indeterminate form 'zero times infinity' occurs. Equation (45) is more convenient. When temperature is low, T → 0, the cotangent function can be approximated as follows: We insert this expression into Equation (45) and obtain: where: is the average kinetic energy for temperature T = 0, i.e., when the thermostat is in a vacuum state and: is the first correction for small temperature T > 0. The kinetic energy E 0 at T = 0 is finite for all finite values of parameters. Its parameters dependence is visualized in Figure 1, and it follows that: (i) E 0 increases monotonically from zero to infinity when the coupling strength µ 0 increases from zero to infinity. (ii) When the decay rate ε = 1/τ c grows from zero to infinity, E 0 grows from zero to infinity. (iii) For Ω = 0 (the Drude model), E 0 starts from some maximal value for a given set of parameters, and next, it decreases when Ω increases.

Regime of Long Memory Time
The damping kernel γ(t) in the Langevin Equation (12) describes memory effects determined by the relaxation (decay) time τ c . For time scales shorter than τ c , memory effects can play an important role. For times longer than τ c , memory effects can be neglected, and the Ohmic model can be applied. Now, we consider the case of a long decay time τ c . More precisely, we assume that τ c is much longer than the thermal Matsubara timeh/k B T, namely, In other words,hε << 2πk B T, and then,hε + 2πnk B T ≈ 2πnk B T in Equation (50). In this regime, Equation (50) takes the form: We can use Formula (49) to rewrite the above equation in a more compact form as: For the Drude model, when Ω = 0, it reduces to the following equation: This is a surprising result because it looks like Equation (2) for the averaged kinetic energy of the oscillator with its redefined eigenfrequency ω 0 = √ εµ 0 . Remember that the relation Equation (57) should be satisfied, and this means that: e.g., for a temperature of 1 Kelvin, τ c >> 10 −12 s, while for 10 −4 Kelvin, τ c >> 10 −8 s. Therefore, for higher temperatures, it is easier to fulfil this condition.

Conclusions
In summary, in the framework of the Generalized Langevin Equation (GLE), we studied the kinetic energy of a quantum Brownian particle in an equilibrium state. We assumed a relatively general form of the integral kernel of GLE and analyzed kinetic energy in selected regimes like high or low temperature limits. In the high temperature limit, the standard result of the classical statistical mechanics of the equipartition energy E k = k B T/2 is valid independently of the strength of the system-environment interaction, decay of the dissipation kernel and its frequency parameter. In the zero temperature regime, when fluctuations of the environment vacuum influence the system, average kinetic energy of the free Brownian particle is non-zero, and its value depends on parameters of the dissipation function. In particular, it is an increasing function of the coupling strength quantified by the parameter γ 0 or the rescaled parameter µ 0 . It is also interesting that when the decay time τ c in the dissipation kernel tends to zero, average kinetic energy grows to infinity [19]. This means that the assumption of zero decay time is not physically correct, and memory effects have to be taken into account.
We propose a new interpretation of the relation Equation (45): the average kinetic equation of the Brownian particle equals the thermally-averaged kinetic energy per one degree of freedom of thermostat oscillators, additionally averaged over randomly-distributed oscillator frequencies.
Our conjecture is that this interpretation is not only for this particular system, but it may be universal, at least for some classes of systems.
We considered an exactly solved model of a quantum open system. In a general case, only approximate results can be derived, e.g., in the so-called quantum Smoluchowski regime, an effective evolution equation of the Fokker-Planck type has been derived [21,22] and applied to many problems such as quantum diffusion [23,24] or quantum Brownian motors [25]. Some extensions to include reservoirs consisting of non-linear oscillators have been proposed [26]. Finally, it is worth mentioning a novel and alternative approach to attack the problem of a quantum Brownian particle, which is based on an adjoint master equation for a generic operator of the system [27].
The integrand of the zero-th term I 0 is the same as Equation (46). The denominator of the integrand has six roots on the complex plane with three poles in the upper half-plane and three poles in the lower half-plane. We can choose a contour closed in the upper half-plane, and from the residue method, we obtain: For remaining terms, we obtain the expression for I n in the following form: Finally, 2E k k B T = 1 + ∞ ∑ n=1 2μ 0ε (ε + 2πn) µ 0ε (ε + 2πn) + 2πn (ε + 2πn) 2 + 2πΩ 2 n . (A25) Using this method, we obtain a form that is more convenient for discussing the general properties of kinetic energy with respect to system parameters. We presented this in Section 7.