(5.89)
Now let the response be still linear, but have an arbitrary frequency dispersion,
F
() q ,
(5.90)
where the function (), called the generalized susceptibility (in our case, of the environment) may be complex, i.e. have both the imaginary and real parts:
() ' ()
i " () .
(5.91)
Then Eq. (73) remains valid with the replacement ” ()/: 45
"
S
( )
T .
(5.92)
F
This fundamental relation46 may be used not only to calculate the fluctuation intensity from the known generalized responsibility (i.e. the deterministic response of the system to a small perturbation), but also in reverse – to calculate such linear response from the known fluctuations. The latter use is especially attractive at numerical simulations of complex systems, e.g., those based on molecular-dynamics approaches, because it circumvents the need in extracting a weak response to a small perturbation out of a noisy background.
Now let us discuss what generalization of Eq. (92) is necessary to make that fundamental result suitable for arbitrary temperatures, T ~ . The calculations we had performed were based on the apparently classical equation of motion, Eq. (63). However, quantum mechanics shows47 that a similar equation is valid for the corresponding Heisenberg-picture operators, so that repeating all the arguments leading to the Langevin equation (65), we may write its quantum-mechanical version
ˆ~
m ˆ q ˆ q ˆ q
Heisenberg-
Fˆ
.
(5.93)
det
F
Langevin
equation
This is the so-called Heisenberg-Langevin (or “quantum Langevin”) equation – in this particular case, for a harmonic oscillator.
The further operations, however, require certain caution, because the right-hand side of the equation is now an operator, and has some nontrivial properties. For example, the “values” of the Heisenberg operator, representing the same variable f( t) at different times, do not necessarily commute: ˆ~
ˆ~
ˆ~
ˆ~
f t
( ) f t'
( ) f t'
( ) f t
( ),
if t'
t .
(5.94)
45 Reviewing the calculations leading to Eq. (73), we may see that the possible real part ’() of the susceptibility just adds up to ( k – m2) in the denominator of Eq. (67), resulting in a change of the oscillator’s frequency 0.
This renormalization is insignificant if the oscillator-to-environment coupling is weak, i.e. if the susceptibility
() is small – as had been assumed at the derivation of Eq. (69) and hence Eq. (73).
46 It is sometimes called the Green-Kubo (or just the Kubo) formula. This is hardly fair, because, as the reader could see, Eq. (92) is just an elementary generalization of the Nyquist formula (81). Moreover, the corresponding works of M. Green and R. Kubo were published, respectively, in 1954 and 1957, i.e. after the 1951 paper by H.
Callen and T. Welton, where a more general result (98) had been derived. Much more adequately, the Green/Kubo names are associated with Eq. (102) below.
47 See, e.g., QM Sec. 4.6.
Chapter 5
Page 22 of 44
SM: Statistical Mechanics
As a result, the function defined by Eq. (46) may not be a symmetric function of the time delay t’ – t even for a stationary process, making it inadequate for the representation of the actual correlation function – which has to obey Eq. (50). This technical difficulty may be overcome by the introduction of the following symmetrized correlation function 48
1 ˆ~
ˆ~
ˆ~
ˆ~
1
ˆ~
ˆ~
K ( )
f ( t) f ( t ) f ( t ) f ( t)
f ( t), f ( t )
,
(5.95)
f
2
2
(where {…,…} denotes the anticommutator of the two operators), and, similarly, the symmetrical spectral density Sf(), defined by the following relation:
1
1
ˆ ˆ *
ˆ *
S () ( ')
f f
ˆ
ˆ
ˆ
,
,
(5.96)
f
ω'
fω' f
*
f
f
ω'
2
2
with Kf() and Sf() still related by the Fourier transform (59).
Now we may repeat all the analysis that was carried out for the classical case, and get Eq. (71) again, but now this expression has to be compared not with the equipartition theorem, but with its quantum-mechanical generalization (14), which, in our current notation, reads
q~2
0
coth
0 .
(5.97)
2
T
2
As a result, we get the following quantum-mechanical generalization of Eq. (92):
"
FDT
S ()
coth .
(5.98)
F
2
T
2
This is the much-celebrated fluctuation-dissipation theorem, usually referred to just as the FDT, first derived in 1951 by Herbert Bernard Callen and Theodore A. Welton – in a somewhat different way.
As natural as it seems, this generalization of the relation between fluctuations and dissipation poses a very interesting conceptual dilemma. Let, for the sake of clarity, temperature be relatively low, T
<< ; then Eq. (98) gives a temperature-independent result
Quantum
"
noise
S
( )
,
(5.99)
F
2
which describes what is frequently called quantum noise. According to the quantum Langevin equation (93), nothing but the random force exerted by the environment, with the spectral density (99) proportional to the imaginary part of susceptibility (i.e. damping), is the source of the ground-state
“fluctuations” of the coordinate and momentum of a quantum harmonic oscillator, with the r.m.s. values 1/ 2
1/ 2
1/ 2
1/ 2
m
~2
~
q q
,
2
0
p p
,
(5.100)
2
m
2
0
48 Here (and to the end of this section) the averaging … should be understood in the general quantum-statistical sense – see Eq. (2.12). As was discussed in Sec. 2.1, for the classical-mixture state of the system, this does not create any difference in either the mathematical treatment of the averages or their physical interpretation.
Chapter 5
Page 23 of 44
SM: Statistical Mechanics
and the total energy 0/2. On the other hand, the basic quantum mechanics tells us that exactly these formulas describe the ground state of a dissipation-free oscillator, not coupled to any environment, and are a direct corollary of the basic commutation relation
q ˆ, p ˆ i .
(5.101)
So, what is the genuine source of the uncertainty described by Eqs. (100)?
The best resolution of this paradox I can offer is that either interpretation of Eqs. (100) is legitimate, with their relative convenience depending on the particular application. One may say that since the right-hand side of the quantum Langevin equation (93) is a quantum-mechanical operator, rather than a classical force, it “carries the uncertainty relation within itself”. However, this (admittedly, opportunistic:-) resolution leaves the following question open: is the quantum noise (99) of the environment’s observable F directly, without any probe oscillator subjected to it? An experimental resolution of this dilemma is not quite simple, because usual scientific instruments have their own ground-state uncertainty, i.e. their own quantum fluctuations, which may be readily confused with those of the system under study. Fortunately, this difficulty may be overcome, for example, using unique frequency-mixing (“down-conversion”) properties of Josephson junctions. Special low-temperature experiments using such down-conversion49 have confirmed that the noise (99) is real and measurable.
Finally, let me mention an alternative derivation50 of the fluctuation-theorem (98) from the general quantum mechanics of open systems. This derivation is substantially longer than that presented above, but gives an interesting sub-product, the Green-Kubo formula
ˆ~
ˆ~
F ( t),F ( t ) iG( )
,
(5.102)
where G() is the temporal Green’s function of the environment, defined by the following relation:
t
F t
( ) G() q t( ) d G t( t' ) q t'
( ) dt' .
(5.103)
0
Plugging the Fourier transforms of all three functions of time participating in Eq. (103) into this relation, it is straightforward to check that this Green’s function is just the Fourier image of the complex susceptibility () defined by Eq. (90):
( )
ei
G
d ()
;
(5.104)
0
here 0 is used as the lower limit instead of (–) just to emphasize that due to the causality principle, Green’s function has to be equal zero for < 0.51
In order to reveal the real beauty of Eq. (102), we may use the Wiener-Khinchin theorem (59) to rewrite the fluctuation-dissipation theorem (98) in a similar time-domain form:
49 R. Koch et al., Phys. Rev. B 26, 74 (1982), and references therein.
50 See, e.g., QM Sec. 7.4.
51 See, e.g., CM Sec. 5.1.
Chapter 5
Page 24 of 44
SM: Statistical Mechanics
ˆ~
ˆ~
F ( t),F ( t ) 2 K ( )
F
,
(5.105)
where the symmetrized correlation function K () is most simply described by its Fourier transform, F
which is, according to Eq. (58), equal to S (), so that using the FDT, we get
F
"
K ( ) cos d
coth
.
(5.106)
F
2
T
2
0
The comparison of Eqs. (102) and (104), on one hand, and Eqs (105)-(106), on the other hand, shows that both the commutation and anticommutation properties of the Heisenberg-Langevin force operator at different moments of time are determined by the same generalized susceptibility () of the environment. However, the averaged anticommutator also depends on temperature, while the averaged commutator does not – at least explicitly, because the complex susceptibility of an environment may be temperature-dependent as well.
5.6. The Kramers problem and the Smoluchowski equation
Returning to the classical case, it is evident that Langevin equations of the type (65) provide means not only for the analysis of stationary fluctuations, but also for the description of arbitrary time evolution of (classical) systems coupled to their environments – which, again, may provide both dissipation and fluctuations. However, this approach to evolution analysis suffers from two major handicaps.
First, the Langevin equation does enable a straightforward calculation of the statistical average of the variable q, and its fluctuation variance – i.e., in the common mathematical terminology, the first and second moments of the probability density w( q, t) – as functions of time, but not of the probability distribution as such. Admittedly, this is rarely a big problem, because in most cases the distribution is Gaussian – see, e.g., Eq. (2.77).
The second, more painful drawback of the Langevin approach is that it is instrumental only for
“linear” systems – i.e., the systems whose dynamics may be described by linear differential equations, such as Eq. (65). However, as we know from classical dynamics, many important problems (for example, the Kepler problem of planetary motion52) are reduced to motion in substantially non-harmonic potentials U ef( q), leading to nonlinear equations of motion. If the energy of interaction between the system and its random environment is factorable – i.e. is a product of variables belonging to these subsystems (as it is very frequently the case), we may repeat all arguments of the last section to derive the following generalized version of the 1D Langevin equation: 53
U
( q, t)
~
q
m q
F ( t) ,
(5.107)
q
52 See, e.g., CM Secs. 3.4-3.6.
53 The generalization of Eq. (107) to higher spatial dimensionality is also straightforward, with the scalar variable q replaced with a multi-dimensional vector q, and the scalar derivative dU/ dq replaced with the vector U, where
is the del vector-operator in the q-space.
Chapter 5
Page 25 of 44
SM: Statistical Mechanics
valid for an arbitrary, possibly time-dependent potential U( q, t). Unfortunately, the solution of this equation may be very hard. Indeed, its Fourier analysis carried out in the last section was essentially based on the linear superposition principle, which is invalid for nonlinear equations.
If the fluctuation intensity is low, q << q, where q( t) is the deterministic solution of Eq.
(107) in the absence of fluctuations, this equation may be linearized54 with respect to small fluctuations q~ q q to get a linear equation,
2
~
q
m~ q~ t
( q~
) F t(),
with t
( )
U q t(), t.
(5.108)
q 2
This equation differs from Eq. (65) only by the time dependence of the effective spring constant ( t), and may be solved by the Fourier expansion of both the fluctuations and the function ( t). Such calculations may be more cumbersome than those performed above, but still be doable (especially if the unperturbed motion q( t) is periodic), and sometimes give useful analytical results.55
However, some important problems cannot be solved by linearization. Perhaps, the most apparent (and practically very important) example is the so-called Kramers problem 56 of finding the lifetime of a metastable state of a 1D classical system in a potential well separated from the region of unlimited motion with a potential barrier – see Fig. 10.
U ( q)
Iw
U 0
0
q
q
Fig. 5.10. The Kramers problem.
1 q'
2 q"
q
In the absence of fluctuations, the system, initially placed close to the well’s bottom (in Fig. 10, at q q 1), would stay there forever. Fluctuations result not only in a finite spread of the probability density w( q, t) around that point but also in a gradual decrease of the total probability W ( t) (
w q, t) dq
(5.109)
s
well'
bottom
to find the system in the well, because of a non-zero rate of its escape from it, over the potential barrier, due to thermal activation. What may be immediately expected of the situation is that if the barrier height,
U U ( q ) U ( q ) ,
(5.110)
0
2
1
54 See, e.g., CM Secs. 3.2, 5.2, and beyond.
55 See, e.g., QM Problem 7.8, and also Chapters 5 and 6 in the monograph by W. Coffey et al., cited above.
56 It was named after Hendrik Anthony (“Hans”) Kramers who, besides solving this conceptually important problem in 1940, has made several other seminal contributions to physics, including the famous Kramers-Kronig dispersion relations (see, e.g., EM Sec. 7.4) and the WKB (Wentzel-Kramers-Brillouin) approximation in quantum mechanics – see, e.g., QM Sec. 2.4.
Chapter 5
Page 26 of 44
SM: Statistical Mechanics
is much larger than temperature T,57 the Boltzmann distribution w exp{- U( q)/ T} should be still approximately valid in most of the well, so that the probability for the system to overcome the barrier should scale as exp{- U 0/ T}. From these handwaving arguments, one may reasonably expect that if the probability W( t) of the system’s still residing in the well by time t obeys the usual “decay law”
W
W
,
(5.111a)
then the lifetime has to obey the general Arrhenius law:
U 0
exp
.
(5.111b)
A
T
However, these relations need to be proved, and the pre-exponential coefficient A (usually called the
attempt time) needs to be calculated. This cannot be done by the linearization of Eq. (107), because this approximation is equivalent to a quadratic approximation of the potential U( q), which evidently cannot describe the potential well and the potential barrier simultaneously – see Fig. 10 again.
This and other essentially nonlinear problems may be addressed using an alternative approach to fluctuations, dealing directly with the time evolution of the probability density w( q, t). Due to the shortage of time/space, I will review this approach using mostly handwaving arguments, and refer the interested reader to special literature58 for strict mathematical proofs. Let us start from the diffusion of a free classical 1D particle with inertial effects negligible in comparison with damping. It is described by the Langevin equation (74) with Fdet = 0. Let us assume that at all times the probability distribution stays Gaussian:
1
q q 2
0
(
w q, t)
,
(5.112)
2
exp
1/ 2
2
q
( t)
2 q
( t)
where q 0 is the initial position of the particle, and q( t) is the time-dependent distribution width, whose growth in time is described, as we already know, by Eq. (77):
q
( t) 2 Dt1/2 .
(5.113)
Then it is straightforward to verify, by substitution, that this solution satisfies the following simple partial differential equation,59
2
w
w
D
,
(5.114)
2
t
q
with the delta-functional initial condition
(
w q )
0
, ( q q ) .
(5.115)
0
57 If U 0 is comparable with T, the system’s behavior also depends substantially on the initial probability distribution, i.e., does not follow the simple law (111).
58 See, e.g., either R. Stratonovich, Topics in the Theory of Random Noise, vol. 1., Gordon and Breach, 1963, or Chapter 1 in the monograph by W. Coffey et al., cited above.
59 By the way, the goal of the traditional definition (78) of the diffusion coefficient, leading to the front coefficient 2 in Eq. (77), is exactly to have the fundamental equations (114) and (116) free of numerical coefficients.
Chapter 5
Page 27 of 44
SM: Statistical Mechanics
The simple and important equation of diffusion (114) may be naturally generalized to the 3D motion:60
w
Equation
D 2
w .
(5.116) of diffusion
t
Now let us compare this equation with the probability conservation law,61
w
j 0 ,
(5.117a)
w
t
where the vector j w has the physical sense of the probability current density. (The validity of this relation is evident from its integral form,
d
3
2
wd r j d q 0
,
(5.117b)
w
dt V
S
which results from the integration of Eq. (117a) over an arbitrary time-independent volume V limited by surface S, and applying the divergence theorem62 to the second term.) The continuity relation (117a) coincides with Eq. (116), with D given by Eq. (78), only if we take
T
j D w w .
(5.118)
w
The first form of this relation allows a simple interpretation: the probability flow is proportional to the spatial gradient of the probability density (i.e., in application to N >> 1 similar and independent particles, just to the gradient of their concentration n = Nw), with the sign corresponding to the flow from the higher to lower concentrations. This flow is the very essence of the effect of diffusion. The second form of Eq. (118) is also not very surprising: the diffusion speed scales as temperature and is inversely proportional to the viscous drag.
The fundamental Eq. (117) has to be satisfied also in the case of a force-driven particle at negligible diffusion ( D 0); in this case
j wv ,
(5.119)
w
where v is the deterministic velocity of the particle. In the high-damping limit we are considering right now, v has to be just the drift velocity:
1
1
v F U (q) ,
(5.120)
det
where Fdet is the deterministic force described by the potential energy U(q).
Now that we have descriptions of j w due to both the drift and the diffusion separately, we may rationally assume that in the general case when both effects are present, the corresponding components (118) and (119) of the probability current just add up, so that
60 As will be discussed in Chapter 6, the equation of diffusion also describes several other physical phenomena –
in particular, the heat propagation in a uniform, isotropic solid, and in this context is called the heat conduction equation or (rather inappropriately) just the “heat equation”.
61 Both forms of Eq. (117) are similar to the mass conservation law in classical dynamics (see, e.g., CM Sec. 8.2), the electric charge conservation law in electrodynamics (see, e.g., EM Sec. 4.1), and the probability conservation law in quantum mechanics (see, e.g., QM Sec. 1.4).
62 See, e.g., MA Eq. (12.2),
Chapter 5
Page 28 of 44
SM: Statistical Mechanics
1
j
,
(5.121)
w
w
U T
w
and Eq. (117a) takes the form
w
Smoluchowski
w U T 2
w .
(5.122)
equation
t
This is the Smoluchowski equation,63 which is closely related to the drift-diffusion equation in multiparticle kinetics – to be discussed in the next chapter.
As a sanity check, let us see what the Smoluchowski equation gives in the stationary limit, w/ t
0 (which evidently may be eventually achieved only if the deterministic potential U is time-independent.) Then Eq. (117a) yields j w = const, where the constant describes the deterministic motion of the system as the whole. If such a motion is absent, j w = 0, then according to Eq. (121),
w
U
w U T w ,
0
i.e.
.
(5.123)
w
T
Since the left-hand side of the last relation is just (ln w), it may be easily integrated over q, giving U
U (q)
ln w
ln C,
w
i.e. (r) C exp
,
(5.124)
T
T
where C is a normalization constant. With both sides multiplied by the number N of similar, independent systems, with the spatial density n(q) = Nw(q), this equality becomes the Boltzmann distribution (3.26).
As a less trivial example of the Smoluchowski equation’s applications, let us use it to solve the 1D Kramers problem (Fig. 10) in the corresponding high-damping limit, m << A, where A (still to be calculated) is some time scale of the particle’s motion inside the well. It is straightforward to verify that the 1D version of Eq. (121),
1 U
w
I
,
(5.125a)
w
w
T
q
q
(where Iw is the probability current at a certain point q, rather than its density) is mathematically equivalent to
T
U ( q)
U ( q)
I exp
exp
,
(5.125b)
w
w
T q
T
so that we may write
U ( q)
T
U ( q)
I exp
exp
.
(5.126)
w
w
T
q
T
As was discussed above, the notion of metastable state’s lifetime is well defined only for sufficiently low temperatures
63 It is named after Marian Smoluchowski, who developed this formalism in 1906, apparently independently from the slightly earlier Einstein’s work, but in much more detail. This equation has important applications in many fields of science – including such surprising topics as statistics of spikes in neural networks. (Note, however, that in some non-physical fields, Eq. (122) is referred to as the Fokker-Planck equation, while actually, the latter equation is much more general – see the next section.)
Chapter 5
Page 29 of 44
SM: Statistical Mechanics
T U .
(5.127)
0
when the lifetime is relatively long: >> A. Since according to Eq. (111a), the first term of the continuity equation (117b) has to be of the order of W/, in this limit the term, and hence the gradient of Iw, are exponentially small, so the probability current virtually does not depend on q in the potential barrier region. Let us use this fact at the integration of both sides of Eq. (126) over that region: q"
U
( q)
T
U
( q q"
I
)
exp
dq
w exp
,
(5.128)
w
T
T
q'
q'
where the integration limits q’ and q” (see Fig. 10) are selected so that T U ( q' ) U ( q ), U ( q ) U ( q" ) U .
(5.129)
1
2
0
(Obviously, such selection is only possible if the condition (127) is satisfied.) In this limit, the contribution from the point q” to the right-hand side of Eq. (129) is negligible because the probability density behind the barrier is exponentially small. On the other hand, the probability at the point q’ has to be close to the value given by its quasi-stationary Boltzmann distribution (124), so that
U ( q' )
U ( q )
1
w( q' ) exp
w( q )exp
,
(5.130)
1
T
T
and Eq. (128) yields
q"
T
U ( q) U ( q )
1
I
w( q )
exp
.
(5.131)
w
1
dq
q'
T
Patience, my reader, we are almost done. The probability density w( q 1) at the well’s bottom may be expressed in terms of the total probability W of the particle being in the well by using the normalization condition
W
U ( q ) U ( q)
1
(
w q ) exp
dq ;
(5.132)
1
T
s
well'
bottom
the integration here may be limited to the region where the difference U( q) – U( q 1) is much smaller than U 0 – cf. Eq. (129). According to the Taylor expansion, the shape of virtually any smooth potential U( q) near the point q 1 of its minimum may be well approximated with a quadratic parabola: 2
d U
U ( q q ) U ( q )
1
( q q )2 ,
where
0.
(5.133)
1
1
2
1
1
2
q q 1
dq
With this approximation, Eq. (132) is reduced to the standard Gaussian integral:64
2
~
1/ 2
( q q )
2
q
2 T
1
1
1
~
W (
w q )
exp
dq (
w q ) exp
q
d (
w q )
.
(5.134)
1
2
1
T
2
1
T
s
well'
1
bottom
To complete the calculation, we may use a similar approximation for the barrier top:
64 If necessary, see MA Eq. (6.9b) again.
Chapter 5
Page 30 of 44
SM: Statistical Mechanics
U ( q q ) U ( q ) U ( q )
2
( q q )2 U ( q )
2
U
( q q )2 ,
2
1
2
2
2
1
0
2
2
(5.135)
2
d U
where
,
0
2
2
q q 2
dq
and work out the remaining integral in Eq. (131), because in the limit (129) it is dominated by the contribution from a region very close to the barrier top, where the approximation (135) is asymptotically exact. As a result, we get
q"
1/ 2
U ( q) U ( q )
U
2 T
1
0
exp
dq exp
.
(5.136)
q'
T
T
2
Plugging Eq. (136), and the w( q 1) expressed from Eq. (134), into Eq. (131), we finally get 1/ 2
( )
1 2
U 0
I W
exp
.
(5.137)
w
2
T
This expression should be compared with the 1D version of Eq. (117b) for the segment [–, q’].
Since this interval covers the region near q 1 where most of the probability density resides, and Iq(-) =
0, this equation is merely
dW
I ( q' ) 0 .
(5.138)
dt
w
In our approximation, Iw( q’) does not depend on the exact position of the point q’, and is given by Eq.
(137), so that plugging it into Eq. (138), we recover the exponential decay law (111a), with the lifetime
obeying the Arrhenius law (111b), and the following attempt time:
Kramers
2
formula
2
,
where
.
(5.139)
A
for high
1/ 2
1 2 1/ 2
damping
1 2
,
1 2
,
1 2
Thus the metastable state lifetime is indeed described by the Arrhenius law, with the attempt time scaling as the geometric mean of the system’s “relaxation times” near the potential well bottom (1) and the potential barrier top (2).65 Let me leave for the reader’s exercise to prove that if the potential profile near well’s bottom and/or top is sharp, the expression for the attempt time should be modified, but the Arrhenius decay law (111) is not affected.
5.7. The Fokker-Planck equation
Formula (139) is just a particular, high-damping limit of a more general result obtained by Kramers. In order to get all of it (and much more), we need to generalize the Smoluchowski equation to arbitrary values of damping . In this case, the probability density w is a function of not only the particle’s position q (and time t) but also of its momentum p – see Eq. (2.11). Thus the continuity equation (117) needs to be generalized to the 6D phase space {q, p}. Such generalization is natural: 65 Actually, 2 describes the characteristic time of the exponential growth of small deviations from the unstable fixed point q 2 at the barrier top, rather than their decay, as near the stable point q 1.
Chapter 5
Page 31 of 44
SM: Statistical Mechanics
w
j j 0,
(5.140)
q
q
p
p
t
where j q (which was called j w in the last section) is the probability current density in the coordinate space, and q (which was denoted as in that section) is the usual vector operator in the space, while j p
is the current density in the momentum space, and p is the similar vector operator in that space: 3
3
n
n
.
(5.141)
q
j
p
j
j1
q
p
j
j1
j
At negligible fluctuations ( T 0), j p may be composed using the natural analogy with j q – see Eq. (119). In our new notation, that relation reads,
p
j wq w ,
(5.142)
q
m
so it is natural to take
j wp w
,
(5.143a)
p
F
where the (statistical-ensemble) averaged force F includes not only the contribution due to the potential’s gradient, but also the drag force –v provided by the environment – see Eq. (64) and its discussion:
p
j w( U v) w( U ) .
(5.143b)
p
q
q
m
As a sanity check, it is straightforward to verify that the diffusion-free equation resulting from the combination of Eqs. (140), (142) and (143),
w
p
p
,
(5.144)
drift
q
w
p
w U
q
t
m
m
allows the following particular solution:
w q
( , p, t) q q t p p t, (5.145)
where the statistical-averaged coordinate and momentum satisfy the deterministic equations of motion, p
p
q
,
p U
,
(5.146)
m
q
m
describing the particle’s drift, with the usual deterministic initial conditions.
In order to understand how the diffusion should be accounted for, let us consider a statistical ensemble of free ( qU = 0, 0) particles that are uniformly distributed in the direct space q (so that
qw = 0), but possibly localized in the momentum space. For this case, the right-hand side of Eq. (144) vanishes, i.e. the time evolution of the probability density w may be only due to diffusion. In the corresponding limit F 0, the Langevin equation (107) for each Cartesian coordinate is reduced to
~
~
q
m F ( t),
i.e. p F ( t) .
(5.147)
j
j
j
j
Chapter 5
Page 32 of 44
SM: Statistical Mechanics
The last equation is identical to the high-damping 1D equation (74) (with Fdet = 0), with the replacement q pj/, and hence the corresponding contribution to w/ t may be described by the last term of Eq.
(122), with that replacement:
w
2
T
D
w
2
2
w T 2
w
t diffusion
p /
.
(5.148)
p
p
Now the reasonable assumption that in the arbitrary case the drift and diffusion contributions to w/ t just add up immediately leads us to the full Fokker-Planck equation:66
Fokker-
w
p
Planck
w w U p T 2
w .
(5.149)
equation
t
q
m
p
q
m
p
As a sanity check, let us use this equation to calculate the stationary probability distribution of the momentum of particles with an arbitrary damping but otherwise free, in the momentum space, assuming (just for simplicity) their uniform distribution in the direct space, q = 0. In this case, Eq.
(149) is reduced to
p
p
2
.
(5.150)
p
w
T w ,
0
i.e.
p
p
w T
w
p
0
m
m
The first integration over the momentum space yields
p
p 2
w T w j ,
i.e. w
T w j ,
(5.151)
m
p
w
p 2 m
p
w
where j w is a vector constant describing a possible general probability flow in the system. In the absence of such flow, j w = 0, we get
p 2
w
2
2
p
p
p
T
T ln w ,
0
giving w const exp
,
(5.152)
p
p
2 m
w
2 m
2 mT
i.e. the Maxwell distribution (3.5). However, the result (152) is more general than that obtained in Sec.
3.1, because it shows that the distribution stays the same even at non-zero damping. It is easy to verify that in the more general case of an arbitrary stationary potential U(q), Eq. (149) is satisfied with the stationary solution (3.24), also giving j w = 0.
It is also straightforward to show that if the damping is large (in the sense assumed in the last section), the solution of the Fokker-Planck equation tends to the following product
2
p
(
w q,p, t) const exp
w(q ,t) ,
(5.153)
2 mT
where the direct-space distribution w(q, t) obeys the Smoluchowski equation (122).
Another important particular case is that of a quasi-periodic motion of a particle, with low damping, in a soft potential well. In this case, the Fokker-Planck equation describes both diffusion of the effective phase of such (generally nonlinear, “anharmonic”) oscillator, and slow relaxation of its 66 It was first derived by Adriaan Fokker in 1913 in his PhD thesis, and further elaborated by Max Planck in 1917.
(Curiously, A. Fokker is more famous for his work on music theory, and the invention and construction of several new keyboard instruments, than for this and several other important contributions to theoretical physics.) Chapter 5
Page 33 of 44
SM: Statistical Mechanics
energy. If we are only interested in the latter process, Eq. (149) may be reduced to the so-called energy diffusion equation,67 which is easier to solve.
However, in most practically interesting cases, solutions of Eq. (149) are rather complicated.
(Indeed, the reader should remember that these solutions embody, in the particular case T = 0, all classical dynamics of a particle.) Because of this, I will present (rather than derive) only one more of them: the solution of the Kramers problem (Fig. 10). Acting almost exactly as in Sec. 6, one can show68
that at virtually arbitrary damping (but still in the limit T << U 0), the metastable state’s lifetime is again given by the Arrhenius formula (111b), with the attempt time again expressed by the first of Eqs. (139), but with the reciprocal time constants 1/1,2 replaced with
2
1/2
η
m
'
2
,
for
,
(5.154)
,
1 2
,
1 2
,12
,
1 2
2 m
2 m
1/ , for m ,
η
,
1 2
,
1 2
where 1,2 (1,2/ m)1/2, and 1,2 are the effective spring constants defined by Eqs. (133) and (135). Thus, in the important particular limit of low damping, Eqs. (111b) and (154) give the famous formula Kramers
2
U 0
formula
exp
.
(5.155)
for low
1
2 1/ 2
T
damping
This Kramers’ result for the classical thermal activation of the dissipation-free system over a potential barrier may be compared with that for its quantum-mechanical tunneling through the barrier.69
The WKB approximation for the latter effect gives the expression
2
2
( q)
exp
,
(5.156)
Q
A
2
( q) dq
,
with
U ( q) E
2 m
2
( )0
q
showing that generally, the classical and quantum lifetimes of a metastable state have different dependences on the barrier shape. For example, for a nearly-rectangular potential barrier, the exponent that determines the classical lifetime (155) depends (linearly) only on the barrier height U 0, while that defining the quantum lifetime (156) is proportional to the barrier width and to the square root of U 0.
However, in the important case of “soft” potential profiles, which are typical for the case of barely emerging (or nearly disappearing) quantum wells (Fig. 11), the classical and quantum results are closely related.
U ( q)
q
1
U 0
q
q
0
2
Fig. 5.11. Cubic-parabolic potential
profile and its parameters.
67 An example of such an equation, for the particular case of a harmonic oscillator, is given by QM Eq. (7.214).
The Fokker-Planck equation, of course, can give only its classical limit, with n, n e >> 1.
68 A detailed description of this calculation (first performed by H. Kramers in 1940) may be found, for example, in Sec. III.7 of the review paper by S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943).
69 See, e.g., QM Secs. 2.4-2.6.
Chapter 5
Page 34 of 44
SM: Statistical Mechanics
Indeed, such potential profile U( q) may be well approximated by four leading terms of its Taylor expansion, with the highest term proportional to ( q – q 0)3, near any point q 0 in the vicinity of the well. In this approximation, the second derivative d 2 U/ dq 2 vanishes at the inflection point q 0 = ( q 1 + q 2)/2, exactly between the well’s bottom and the barrier’s top (in Fig. 11, q 1 and q 2). Selecting the origin at this point, as this is done in Fig. 11, we may reduce the approximation to just two terms:70
b 3
U ( q) aq q .
(5.157)
3
(For the particle’s escape into the positive direction of the q-axis, we should have a, b > 0.) An easy calculation gives all essential parameters of this cubic parabola: the positions of its minimum and maximum:
q q a b
(5.158)
2
1
/ 1/2,
the barrier height over the well’s bottom:
1/ 2
3
4 a
U U ( q ) U ( q )
,
(5.159)
0
2
1
3
b
and the effective spring constants at these points:
2
d U
2 ab
.
(5.160)
1
2
1/2
2
dq q ,12
The last expression shows that for this potential profile, the frequencies 1,2 participating in Eq.
(155) are equal to each other, so that this result may be rewritten as
Soft well:
2
U
(
2
)
0
2
ab 1/ 2
thermal
exp
,
with
.
(5.161)
lifetime
T
0
m
0
On the other hand, for the same profile, the WKB approximation (156) (which is accurate when the height of the metastable state energy over the well’s bottom, E – U( q 1) 0/2, is much lower than the barrier height U 0) yields71
Soft well:
2
1/2
36 U
quantum
0
0
exp
.
(5.162)
Q
lifetime
864 U
5
0
0
0
The comparison of the dominating, exponential factors in these two results shows that the thermal activation yields a lower lifetime (i.e., dominates the metastable state decay) if the temperature is above the crossover value
36
T
2
.
7
.
(5.163)
c
0
0
5
70 As a reminder, a similar approximation arises for the P( V) function, at the analysis of the van der Waals model near the critical temperature – see Problem 4.6.
71 The main, exponential factor in this result may be obtained simply by ignoring the difference between E and U( q 1), but the correct calculation of the pre-exponential factor requires taking this difference, 0/2, into account
– see, e.g., the model solution of QM Problem 2.43.
Chapter 5
Page 35 of 44
SM: Statistical Mechanics
This expression for the cubic-parabolic barrier may be compared with the similar crossover for a quadratic-parabolic barrier,72 for which T c = 2 0 6.28 0. We see that the numerical factors for the quantum-to-classical crossover temperature for these two different soft potential profiles are close to each other – and much larger than 1, which could result from a naïve estimate.
5.8. Back to the correlation function
Unfortunately, I will not have time/space to either derive or even review solutions of other problems using the Smoluchowski and Fokker-Planck equations, but have to mention one conceptual issue. Since it is intuitively clear that the solution w(q, p, t) of the Fokker-Planck equation for a system provides the complete statistical information about it, one may wonder how it may be used to find its temporal characteristics that were discussed in Secs. 4-5, using the Langevin formalism. For any statistical average of a function taken at the same time instant, the answer is clear – cf. Eq. (2.11): f q t
( ),p t
( ) f (q,p) w(q,p, t d 3
) qd 3 p ,
(5.164)
but what if the function depends on variables taken at different times, for example as in the correlation function Kf() defined by Eq. (48)?
To answer this question, let us start from the discrete-variable case when Eq. (164) takes the form (2.7), which, for our current purposes, may be rewritten as
f t f W t() .
(5.165)
m
m
m
In plain English, this is a sum of all possible values of the function, each multiplied by its probability as a function of time. But this implies that the average f( t) f( t’) may be calculated as the sum of all possible products fmfm’, multiplied by the joint probability to measure outcome m at moment t, and outcome m’ at moment t’. The joint probability may be represented as a product of Wm( t) by the conditional probability W( m’, t’ m, t). Since the correlation function is well defined only for stationary systems, in the last expression we may take t = 0, i.e. look for the conditional probability as the solution, Wm’(), of the equation describing the system’s probability evolution, at time = t’ – t (rather than t’), with the special initial condition
W (0)
.
(5.166)
m'
m' , m
On the other hand, since the average f( t) f( t +) of a stationary process should not depend on t, instead of Wm(0) we may take the stationary probability distribution Wm(), independent of the initial conditions, which may be found as the same special solution, but at time . As a result, we get Correlation
function:
f t f t f W f W
.
(5.167) discrete
m
m
m' m'
m, m'
system
This expression looks simple, but note that this recipe requires solving the time evolution equations for each Wm’() for all possible initial conditions (166). To see how this recipe works in practice, let us revisit the simplest two-level system (see, e.g., Fig. 4.13, which is reproduced in Fig. 12
72 See, e.g., QM Sec. 2.4.
Chapter 5
Page 36 of 44
SM: Statistical Mechanics
below in a notation more convenient for our current purposes), and calculate the correlation function of its energy fluctuations.
W ( t)
E
1
1
Fig. 5.12. Dynamics of a two-level system.
W 0 t
E 0
0
The stationary probabilities of the system’s states (i.e. their probabilities for ) have been calculated in problems of Chapter 2, and then again in Sec. 4.4 – see Eq. (4.68). In our current notation (Fig. 12),
W
W
0
1
1
1
,
,
/ T
/
T
1 e
e
1
(5.168)
that
so
E W W
0
0
1
.
/ T
e
1
To calculate the conditional probabilities Wm’( ) with the initial conditions (167) (according to Eq.
(168), we need all four of them, for { m, m’} = {0, 1}), we may use the master equations (4.100), in our current notation reading
dW
dW
1
0
W W .
(5.169)
0
1
d
d
Since Eq. (169) conserves the total probability, W 0 + W 1 = 1, only one probability (say, W 1) is an independent variable, and for it, Eq. (169) gives a simple, linear differential equation dW
1 W ,
where
,
(5.170)
1
d
which may be readily integrated for an arbitrary initial condition:
W
( ) W ( )
0 e
W ()
,
(5.171)
1
1
1
1 e
where W 1() is given by the second of Eqs. (168). (It is straightforward to verify that the solution for W 0() may be represented in a similar form, with the corresponding change of the state index.) Now everything is ready to calculate the average E( t) E( t +) using Eq. (167), with fm,m’ = E 0,1.
Thanks to our (smart :-) choice of the energy reference, of the four terms in the double sum (167), all three terms that include at least one factor E 0 = 0 vanish, and we have only one term left to calculate: 2
E( t) E( t ) E W () E W ( )
E W () W 0
e
W () 1
e
1
1
1
1
W (0) 1
1
1
1
1
W (0) 1
1
1
2
1
T
(5.172)
T
e
2
1
e
T
e
e
2 1
/
/
e
1
/
e
1
/
e T
.
1
Chapter 5
Page 37 of 44
SM: Statistical Mechanics
From here and the last of Eqs. (168), the correlation function of energy fluctuations is73
~
~
K ( ) E( t) E( t )
E
E( t) E( t) E( t ) E( t)
/
e T
(5.173)
E( t) E( t ) E 2
2
e
e T
,
/
2
1
so that its variance, equal to KE(0), does not depend on the transition rates and . However, since the rates have to obey the detailed balance relation (4.103), / = exp{/ T}, for this variance we may formally write
KE 0
/ T
e
/
,
(5.174)
2
/ T
e
2
1
/
2
1
2
2
so that Eq. (173) may be represented in a simpler form:
Energy
2
K
( )
e
.
(5.175) fluctuations:
E
2
two-level
system
We see that the correlation function of energy fluctuations decays exponentially with time, with the net rate . Now using the Wiener-Khinchin theorem (58) to calculate its spectral density, we get
2
1
Γ Γ
Γ
Δ
Γ Γ
2
S (
)
Δ
e
cos
.
(5.176)
E
d
2
2
2
Γ
Γ Γ
0
Such Lorentzian dependence on frequency is very typical for discrete-state systems described by master equations. It is interesting that the most widely accepted explanation of the 1/ f noise (also called the “flicker” or “excess” noise), which was mentioned in Sec. 5, is that it is a result of thermally-activated jumps between states of two-level systems with an exponentially-broad statistical distribution of the transition rates . Such a broad distribution follows from the Kramers formula (155), which is approximately valid for the lifetimes of both states of systems with double-well potential profiles (Fig.
13), for a statistical ensemble with a smooth statistical distribution of the energy barrier heights U 0.
Such profiles are typical, in particular, for electrons in disordered (amorphous) solid-state materials, which indeed feature high 1/ f noise.
U
U
0
Fig. 5.13. Typical double-
well potential profile.
0
q
Returning to the Fokker-Planck equation, we may use the following evident generalization of Eq. (167) to the continuous-variable case:
73 The step from the first line of Eq. (173) to its second line utilizes the fact that our system is stationary, so that
E( t + ) = E( t) = E() = const.
Chapter 5
Page 38 of 44
SM: Statistical Mechanics
Correlation
function:
f ( t) f ( t )
3
3
3
3
d qd p d q'd p' f (q , p)
w q , p, f (q ', p ' )
w q ', p ' ,
continuous
,
(5.177)
system were both probability densities are particular values of the equation’s solution with the delta-functional initial condition
(
w q ' , p ' )
0
, (q ' - q) (p ' - p) .
(5.178)
For the Smoluchowski equation valid in the high-damping limit, the expressions are similar, albeit with a lower dimensionality:
f ( t) f ( t )
3
3
d q d q' f (q)
w q, f (q ' )
w q ' ,
,
(5.179)
(
w q ' )
0
, (q ' - q) .
(5.180)
To see this formalism in action, let us use it to calculate the correlation function Kq() of a linear relaxator, i.e. an overdamped 1D harmonic oscillator with m0 << . In this limit, as Eq. (65) shows, the oscillator’s coordinate, averaged over the ensemble of environments, obeys a linear equation,
q q 0 ,
(5.181)
which describes its exponential relaxation from the initial position q 0 to the equilibrium position q = 0, with the reciprocal time constant = /:
t
q t
( ) q e
.
(5.182)
0
The deterministic equation (181) corresponds to the quadratic potential energy U( q) = q 2/2, so that the 1D version of the corresponding Smoluchowski equation (122) takes the form
2
w
w
wq T
.
(5.183)
2
t
q
q
It is straightforward to check, by substitution, that this equation, rewritten for the function w( q’,), with the 1D version of the delta-functional initial condition (180), w( q’,0) = ( q’ – q), is satisfied with a Gaussian function:
2
1
q' q ( )
(
w q' , )
,
(5.184)
2
exp
1/ 2
2
q( )
2 q ( )
with its center q() moving in accordance with Eq. (182), and a time-dependent variance 2
2
2
2
2
T
q
( ) q ()
1 e
,
where q () q
.
(5.185)
(As a sanity check, the last equality coincides with the equipartition theorem’s result.) Finally, the first probability under the integral in Eq. (179) may be found from Eq. (184) in the limit (in which
q() 0), by replacing q’ with q:
1
2
q
(
w q,)
(5.186)
2
exp
.
1/ 2 q
()
2 2
q
()
Now all ingredients of the recipe (179) are ready, and we can spell it out, for f ( q) = q, as Chapter 5
Page 39 of 44
SM: Statistical Mechanics
2
1
2
q
q'
qe
q( t) q( t )
dq dq' q exp
q'
.
(5.187)
2
exp
2
2 q( ) q()
2 q ( )
2 q ( )
The integral over q’ may be worked our first, by replacing this integration variable with ( q” + qe-) and hence dq’ with dq” :
1
2
2
q
q"
q( t) q( t )
q exp
dq q" qe
dq" (5.188)
2
exp
.
2 q( ) q
()
2 q
()
2 2
q
( )
The internal integral of the first term in the parentheses equals zero (as that of an odd function in symmetric integration limits), while that with the second term is the standard Gaussian integral, so that 1
q 2
2
T
2
q t
( ) q t
(
)
e
q exp