Thermodynamic theory of equilibrium fluctuations

arXiv:1507.05662v1 [cond-mat.stat-mech] 20 Jul 2015

Y. Mishin
Department of Physics and Astronomy, George Mason University, MSN 3F3,
Fairfax, VA 22030, USA;
Email: ymishin@gmu.edu

July 22, 2015

Abstract
The postulational basis of classical thermodynamics has been expanded to in-
corporate equilibrium fluctuations. The main additional elements of the proposed
thermodynamic theory are the concept of quasi-equilibrium states, a definition of non-
equilibrium entropy, a fundamental equation of state in the entropy representation, and
a fluctuation postulate describing the probability distribution of macroscopic parame-
ters of an isolated system. Although these elements introduce a statistical component
that does not exist in classical thermodynamics, the logical structure of the theory
is different from that of statistical mechanics and represents an expanded version of
thermodynamics. Based on this theory, we present a regular procedure for calculations
of equilibrium fluctuations of extensive parameters, intensive parameters and densi-
ties in systems with any number of fluctuating parameters. The proposed fluctuation
formalism is demonstrated by four applications: (1) derivation of the complete set of
fluctuation relations for a simple fluid in three different ensembles; (2) fluctuations in
finite-reservoir systems interpolating between the canonical and micro-canonical en-
sembles; (3) derivation of fluctuation relations for excess properties of grain boundaries
in binary solid solutions, and (4) derivation of the grain boundary width distribution
for pre-melted grain boundaries in alloys. The last two applications offer an efficient
fluctuation-based approach to calculations of interface excess properties and extraction
of the disjoining potential in pre-melted grain boundaries. Possible future extensions
of the theory are outlined.

Keywords: Thermodynamics, fluctuation, entropy, grain boundary, pre-melting

1 Introduction
1.1 Historical background and goal of the work
Fluctuations of thermodynamic properties play a crucial role in many physical phe-
nomena and diverse applications. Fluctuations are especially important in nanometer-scale

1
systems where they can lead to a large variability of mechanical and functional properties
and create noises affecting performance of devices. Fluctuations are unavoidable in molecu-
lar dynamics and Monte Carlo simulations of materials, where the system dimensions rarely
exceed a few nanometers. In fact, in many atomistic calculations, equilibrium properties of
interest are extracted by analyzing statistical fluctuations of other properties [1–3]. Exam-
ples include the calculations of elastic coefficients of solid materials from strain and/or stress
fluctuations in molecular dynamics [4–7] or Monte Carlo [8, 9] simulations; calculation of
partial molar properties of solutions from concentration fluctuations [10, 11]; and calcula-
tion of the interface free energy of solid-liquid [12, 13] and solid-sold [14–23] interfaces from
capillary fluctuations or fluctuations of the interface width.
Presently, fluctuations are primarily discussed in the statistical-mechanical literature and
typically in the context of specific (usually, very simple) models. At the same time, the ther-
modynamics community traditionally relies on classical thermodynamics [24–28] which, by
the macroscopic and equilibrium nature of this discipline, disregards fluctuations and oper-
ates solely in terms of static properties. In fact, the very term “fluctuation” has a temporal
connotation incompatible with the time-independent character of classical thermodynamics.
There is, however, a third direction in thermal physics that pursues a generalized form of
thermodynamics that incorporates fluctuations of thermodynamic parameters around equi-
librium. By contrast to the statistical-mechanical approach, such theories seek to introduce
fluctuations directly into the thermodynamic framework via additional assumptions, postu-
lates or similar elements of the logical structure. The goal of such theories is to expand the
scope of thermodynamics by introducing statistical elements while preserving the traditional
axiomatic approach that distinguishes this discipline. It is this direction in thermal physics
that constitutes the subject of the present paper.1
The logical foundations of classical thermodynamics have been the subject of research
over the past hundred or more years, starting with the works of Carathéodory [29] and
Ehrenfest [30] and continuing in modern times [26–28, 31–36]. In spite of the fundamental
importance of the traditional three laws of thermodynamics, they are essentially a reflection
of the historical development of the discipline and do not constitute an autonomous and
logically complete structure. The four-postulate structure proposed by Callen [26, 27] is de-
tached from the historical context, deeply thought-through, but still far from complete. The
most rigorous and complete postulational basis of classical thermodynamics has been for-
mulated by Tisza [28]. His theory, called the Macroscopic Thermodynamics of Equilibrium
(MTE), presents an elegant logical structure comprising a set of interconnected definitions,
postulates and corollaries. This rigor comes at a price: Tisza’s MTE is restricted to a certain
class of rather simple thermodynamic systems. Nevertheless, it demonstrates an approach
that can be applied for the construction of similar postulational structures for other classes
of systems.
Unfortunately, attempts to create a similarly rigorous thermodynamic formalism that
1
Although we do not wish to enter into terminological discussions, we point out that an expanded ther-
modynamic theory that incorporates statistical elements such as fluctuations could be called statistical
thermodynamics. This term would distinguish it from both statistical mechanics and classical thermody-
namics. Unfortunately, the term statistical thermodynamics is already used with several different meanings,
most notably as a collective reference to statistical mechanics and thermodynamics (which we refer to as
thermal physics).

2
would include fluctuations have not been very successful. The thermodynamic fluctuation
theory by Greene and Callen [37, 38] and Callen’s Postulate II′ [26]2 (which generalizes his
entropy Postulate II to include fluctuations) turned out to be insufficient. They required
additional assumptions, such as the approximation of average values of thermodynamic
properties by their most probable values, and relied on the assumption that there is no dis-
tinction between the “canonical thermodynamics” and “microcanonical thermodynamics”.
The latter assumption was later criticized by Tisza and Quay [39] from the standpoint of
statistical mechanics. Tisza and Quay’s [39] own treatment of fluctuations was an amalga-
mation of statistical mechanics and thermodynamics, and in this sense, did not achieve the
goal either. In fact, the authors [39] admit that their theory, called the Statistical Thermo-
dynamic of Equilibrium, is only a stepping stone toward a future more general theory whose
development encounters major difficulties.
There are two major challenges in the development of a thermodynamic theory of fluctu-
ations: (1) appropriate definition of non-equilibrium entropy that would seamlessly connect
with the system of other definitions and postulates of thermodynamics, and (2) formulation
of a general fluctuation law of thermodynamic parameters. Historically, the first work ad-
dressing these issues was Einstein’s 1910 paper on critical opalescence [40]. In his theory,
the point of departure was Boltzmann’s equation

S = k ln W + const, (1)

where k is Boltzmann’s constant and W is the probability associated with the entropy S (our
notations are slightly different from Einstein’s). Einstein considered a completely isolated
thermodynamic system with a fixed energy E. He assumed that all equilibrium and non-
equilibrium (fluctuated) thermodynamic states of the system can be described by a set of
macroscopic parameters λ1 , ..., λn . For a discrete set of such parameters, he proposed to
interpreted W as the fractions of time spent by the fluctuating system in each of the states
compatible with the given energy E. Based on this interpretation, Einstein inverted Eq.(1)
to obtain the probability of a macroscopic state and thus a given set of parameters λ1 , ..., λn :
 
S
W = const · exp . (2)
k

For continuos state variables, the probability that the parameters λ1 , ..., λn be found between
λ1 and λ1 + dλ1 ,..., λn and λn + dλn is
 
S − S0
dW = const · exp dλ1 ...dλn , (3)
k

where S0 is the maximum entropy. Einstein’s theory does not make any reference to micro-
states of the system. Instead, it relates the entropy directly to the probabilities of macro-
scopic states. As such, this theory belongs more to the realm of thermodynamics than
statistical mechanics.3
2
Callen’s Postulate II′ [26] was abandoned in a later edition of his book [27].
3
The first statistical-mechanical theory of fluctuations was developed by Gibbs [41], whose work was
apparently unknown to Einstein at the time. Gibbs’ treatment of fluctuations was based explicitly on classical

3
 Equations (2) and (3) essentially constitute a fluctuation law that can be built into
the formalism of thermodynamics. The fluctuations described by these equations occur
in an isolated system, thus the respective entropy has a micro-canonical meaning. For
canonical systems, the distribution function of macroscopic parameters can be derived by
considering a relatively small (but macroscopic) subsystem of an isolated system and treating
the remaining part (complementary system) as a reservoir.
The situation is complicated by the fact that the most advanced thermodynamic theories
of fluctuations [26, 39] postulate the fluctuation law directly in a canonical form and derive
properties of isolated systems as a liming case of canonical. This approach is different from
Einstein’s [40] and is rooted in statistical-mechanical arguments. The treatment of a canon-
ical system as an asymptotic case when the complementary system tends to infinity has
certain conceptual disadvantages in rigorous formulations of statistical mechanics. Alterna-
tive theories [39, 43, 44] postulate the canonical distribution for an infinite reservoir from
the outset and derive properties of isolated and finite systems only later. Thus, postulating
thermodynamic fluctuations in a canonical form is more consistent with the aforementioned
statistical-mechanical theories.
In the present paper, we take the view that, if a thermodynamic fluctuation theory is
to be relatively autonomous, then the simplicity of its internal logical structure is more
important than consistency with all aspects the formal structure of statistical mechanics.
We thus propose to build the theory on a fluctuation law postulated for an isolated system.
The law which will be adopted in this work will look similar to Eq.(3), except for certain
refinements in justification and interpretation.
We emphasize that this paper is focused on equilibrium fluctuations and will not discuss
driven or other strongly non-equilibrium systems. Furthermore, given that our approach
is purely thermodynamic, we are only interested in statistical distributions of fluctuating
parameters without attempting to describe the actual dynamics of the fluctuations [45, 46].

1.2 Organization of the paper

The outline of the paper is a follows. In Secs. 2-7 we formulate the basic assumptions
of the theory describing equilibrium fluctuations in an isolated system. We then consider
the usual quadratic expansion of the entropy around equilibrium and arrive at the Gaussian
law of fluctuations and a set of relations for mean-square fluctuations and covariances of
thermodynamic properties (Sec. 8). In Sec. 9 we derive the fluctuation law for generalized
canonical systems, treating the system of interest and the reservoir as a combined isolated
system. We derive the Gaussian law of canonical fluctuations in two forms, called the entropy
scheme and the energy scheme, and present the fluctuation formalism in both schemes. To
illustrate the calculation techniques, we compute fluctuations of all properties of a simple
fluid in three different ensembles. In Sec. 10 we analyze systems connected to a finite-
size reservoir. The finite-reservoir ensemble interpolates between the micro-canonical and
canonical ensembles and has been implemented in recent computer simulations [21, 47–55].
dynamics of particles, with micro-states defined in the momentum-coordinate phase space. Einstein’s theory
[40] is obviously more general, even though Eqs.(2) and (3) per se do not suggest any means of calculating
the probabilities or the entropy for any particular system. For a more detailed comparison of the Gibbs and
Einstein approaches see Tisza and Quay [39] and a more recent paper by Rudoi and Sukhanov [42].

4
This section demonstrates the ease with which the theory can handle finite-size systems as
opposed to alternate theories [26, 39] built directly on the canonical formalism.
Sections 11 and 12 are devoted to applications of the theory to fluctuations in a single-
phase interface, called a grain boundary (GB), in a binary solid solution. This is a non-trivial
problem, which is addressed by conceptually partitioning the system into an imaginary
perfect crystal (grain) and an excess system representing the GB, each described by their
own fundamental equation. The goal is to describe fluctuations in the excess system. In
Sec. 11 we develop a set of equations linking fluctuations of some interface properties to
equilibrium values of other properties, such as the excess heat capacity, excess elastic moduli
and others. An interesting result here is the prediction of fluctuations of the interface free
energy, an important excess quantity that is usually treated as a static property. Next,
we address the problem of width fluctuations of a pre-melted GB represented as a thin
liquid layer subject to a disjoining pressure. We derive the distribution function of the GB
width expressed through the disjoining potential describing interactions between the two
solid-liquid interfaces bounding the liquid layer. While the functional form of the width
distribution is known for single-component systems [16–18, 20], it has not been previously
derived for binary mixtures. Finally, in Sec. 13 we summarize the paper and outline future
work.
Although the paper is focused on thermodynamics, in Secs. 7 and 9.2 and a few other
parts of the text we do discuss statistical distributions of micro-states and other aspects of
statistical mechanics. It should be emphasized that these discussions are not part of the
proposed fluctuation theory. They are only included here to provide a statistical-mechanical
interpretation of certain thermodynamic statements or demonstrate their consistency with
principles of statistical mechanics. In addition, we slightly deviate from the traditional
macroscopic terminology of classical thermodynamics and talk about the number of particles
(instead of moles) and ensembles (instead of types of imposed constraints). Obviously, this
is only a matter of terminology that does not compromise the thermodynamic nature of the
theory.
Finally, we generally prefer a smooth development of the logic and do not follow the
mathematical style of presentation [28, 29, 39] in which all statements are broken into
definitions, postulates, corollaries and other deductive steps.

2 Definition of equilibrium fluctuations

Consider an isolated thermodynamic system, i.e., a system incapable of exchanging heat
or matter with its environment or performing work on the environment. Some physical
properties of an isolated system are strictly fixed by conservation laws (e.g., the total energy,
total momentum, total amounts of chemical elements4 ), whereas other properties can vary.
Thermodynamics postulates that all non-conserved properties averaged over a certain time
scale tT D , which we call the thermodynamic time scale, eventually stop varying and remain
constant for as long as the system remains isolated. The quiescent state of an isolated
system in which all properties remain constant on the tT D time scale is called the state of
thermodynamic equilibrium.
4
In a system without nuclear reactions.

5
 Even after thermodynamic equilibrium has been reached, non-conserved properties con-
tinue to vary on a shorter time scale, tf ≪ tT D , called the fluctuation time scale. Such
continual variations of properties of an equilibrium system are called equilibrium fluctua-
tions.
For any physical property X, let X be its value averaged over the thermodynamic time
scale tT D . During the equilibrium fluctuations, instantaneous values of X randomly devi-
ate from X. The goal of the fluctuation theory is to predict the probability distribution
of any property X around its average X and compute the moments and other statistical
characteristics of this distribution.

3 The fundamental equation

Consider an isolated thermodynamic system whose equilibrium properties averaged over
the tT D time scale are fully defined by a set of n conserved extensive parameters X1 , ..., Xn .
(Tisza [28] refers to such conserved extensive properties as “additive invariants”.) For ex-
ample, all equilibrium properties of a simple fluid enclosed in an isolated rigid box are fully
defined by its energy E, volume V , and the number of particles N.5 As long as the pa-
rameters X1 , ..., Xn remain fixed, the system remains in a state of equilibrium. One can
temporarily break the isolation, change all or some of the parameters X1 , ..., Xn , and then
isolate and re-equilibrate the system to a new state. By repeating such perturbation/re-
equilibration steps, one can vary the parameters X1 , ..., Xn and measure all equilibrium
thermodynamic properties of the system as functions of X1 , ..., Xn .
In particular, the described procedure can be applied to measure the system entropy S as
a function of the conserved parameters X1 , ..., Xn . Here, the entropy is treated as a property
defined in statistical mechanics. Namely, for an equilibrium isolated system, S is identified
with the micro-canonical entropy S = k ln Ωmax , where Ωmax is the maximum number of
micro-states compatible with the given set of parameters X1 , ..., Xn (see Sec. 7 for a more
detailed statistical-mechanical interpretation of entropy). The entropy S so defined will be
referred to as the “equilibrium entropy” in order to distinguish it from the “non-equilibrium
entropy” defined later.
The function
S = S (X1 , ..., Xn ) (4)
is called the fundamental equation of the system in the entropy representation [24, 26–28].
More accurately, Eq.(4) is the fundamental equation of a particular phase of the system.
Different phases are specified by different fundamental equations (4). For a simple fluid, the
fundamental equation has the form S = S(E, V, N).
It is important to recognize that, at this point of the development, the entropy S has
been defined only for an equilibrium isolated system. Furthermore, the arguments X1 , ..., Xn
of the fundamental equation (4) are understood as conserved properties that remain fixed
(no fluctuations!) in every given state of the isolated system.
5
This can be considered a definition of the simple fluid. Physically, it is a homogeneous single-component
fluid without electric, magnetic or other contributions to its thermodynamic properties.

6
4 The non-equilibrium entropy
When a non-equilibrium isolated system approaches equilibrium, it is often possible to
define a property called the non-equilibrium entropy. Consider two examples.
Example 1. When the system is close enough to equilibrium, it can be mentally parti-
tioned into small but macroscopic subsystems that can be considered as isolated and equi-
librium for a period of time on the order of tq . Such subsystems are called quasi-equilibrium
and the entire isolated system is said to be in a quasi-equilibrium state.6 On the quasi-
equilibrium time scale tq , each subsystem α can be described by a fundamental equation
Sα = Sα (X1α , ..., Xnα ), where X1α , ..., Xnα are extensive properties assumed to be fixed. The
non-equilibrium entropy Ŝ is defined as the sum of the entropies of all quasi-equilibrium
subsystems: X
Ŝ (λ1 , ..., λm , X1 , ..., Xn ) ≡ Sα (X1α , ..., Xnα ) . (5)
α

Here, λ1 , ..., λm are so-called internal parameters describing the distribution of the conserved
extensive properties X1 , ..., Xn of the entire isolated system over its quasi-equilibrium sub-
systems.7 As the system evolves towards equilibrium, these internal parameters vary on
the time scale ≫ tq , as do the local parameters X1α , ..., Xnα of the subsystems and thus the
non-equilibrium entropy Ŝ. Note that the quasi-equilibrium system as a whole does not
have a fundamental equation; its properties depend not only on the total amounts of the
extensive quantities X1 , ..., Xn but also on their distribution among the subsystems.
As a specific example, consider an isolated cylinder with rigid side walls, rigid top and
bottom, and a movable piston dividing the cylinder in two compartments 1 and 2 (Fig. 1).
The compartments are filled with different amounts of the same simple fluid obeying a
fundamental equation S = S(E, V, N). The system is initially not in equilibrium. Heat and
particles can diffuse through the piston, but so slowly that each compartment maintains its
own thermodynamic equilibrium at all times. Thus, as the system drifts toward equilibrium,
it always remains in quasi-equilibrium. During this quasi-equilibrium process, the total
energy E = E1 + E2 , total volume V = V1 + V2 and the total number of particles N =
N1 + N2 remain fixed. These properties are the additive invariants denoted as Xi . The
non-equilibrium entropy of the system equals

Ŝ(E1 , V1 , N1 , E, V, N) = S (E1 , V1 , N1 ) + S (E − E1 , V − V1 , N − N1 ) , (6)

| {z }
λ

where E1 , V1 and N1 play the role of the internal λ-parameters describing the distribution
of the fixed quantities E, V and N between the compartments.
Example 2. Many systems reach equilibrium with respect to some thermodynamic prop-
erties while other properties are still in the process of equilibration. In such cases, we can talk
about quasi-equilibrium as equilibrium with respect to some properties without equilibrium
6
On the time scale tq , the quasi-equilibrium system can be treated as if it was an equilibrium system
with subsystems separated by isolating walls.
7
To avoid notational confusion, we note that our parameters λ1 , ..., λm are independent variables de-
scribing the macroscopic states subject to isolation constraints. In Einstein’s work [40], λ1 , ..., λn are
not independent variables. His Eq.(3) is only valid on the (n − 1)-dimensional constant-energy surface
E(λ1 , ..., λn ) = const in the n-dimensional parameter space of λ1 , ..., λn .

7
 E 1 V1 N 1 E 2 V2 N 2

Figure 1: A model of an isolated system: a rigid cylinder divided in two compartments

containing different amounts of the same substance separated by a movable piston.

with respect to other properties. For example, chemical reactions usually occur much slower
in comparison with thermal and mechanical equilibration. There is a certain time scale tq on
which chemical reactions can be considered as “frozen” and the system can be treated as if it
were in full equilibrium. Accordingly, the system can be assigned a non-equilibrium entropy
Ŝ computed by ignoring the reactions and treating the amounts of all chemical components
as fixed parameters. Of course, on the time scale much longer than tq , the entropy slowly
drifts as the system evolves towards chemical equilibrium. The non-equilibrium entropy so
defined is a function of not only the conserved parameters X1 , ..., Xn but also some progress
variables λ1 , ..., λm describing the extents of the chemical reactions. The latter thus play
the role of the internal parameters specifying the quasi-equilibrium states.
As a simple illustration, consider an isolated box containing a mixture of three chemically
different gases A, B and C. In the absence of chemical reactions, this mixture is described
by a fundamental equation S = S(E, V, NA , NB , NC ). Let the system initially contain the
(i) (i) (i)
amounts of the components NA , NB and NC . Now suppose that the components can
transform to each other via a chemical reaction A + B ⇋ C. Assuming that the reaction is
slow, the system always maintains mechanical and thermal equilibrium, while its chemical
composition gradually changes. Let ∆NA be the change in the number of particles A relative
to the initial state. The non-equilibrium entropy is defined by
(i) (i) (i) (i) (i) (i)
Ŝ(∆NA , E, V, NA , NB , NC ) = S(E, V, NA − ∆NA , NB − ∆NA , NC + ∆NA ). (7)
| {z }
λ

Here, ∆NA characterizes the progress of the chemical reaction and plays the role of an
(i) (i) (i)
internal parameter λ, whereas E, V, NA , NB , NC is a set of fixed extensive parameters
X1 , ..., Xn .

5 The second law of thermodynamics

In the extended thermodynamics including fluctuations, the second law (entropy postu-
late) can be formulated as two statements:

1. The non-equilibrium entropy Ŝ of an isolated system averaged over the thermodynamic

time scale tT D , which we denote S̃, increases with time and reaches a maximum value
S when the system arrives at equilibrium.

8
 2. The individual values of Ŝ fluctuate but can never exceed the equilibrium value S.

Note that this formulation of the second law is stronger than the frequently used weak
formulation. The latter only states that entropy should increase but does not guarantee
that an equilibrium state will be eventually reached.8 Clearly, to talk about equilibrium
fluctuations we must first ensure that the system actually arrives at equilibrium.
Once equilibrium has been reached, the average X of any thermodynamic property
X remains constant while instantaneous values of X can randomly deviate from X on
the fluctuation time scale tf . This includes fluctuations of the non-equilibrium entropy Ŝ.
However, the latter is a special case because it can only fluctuate down from its maximum
value S. As a result, the average entropy S is always slightly smaller than S.
The time evolution of an isolated system approaching equilibrium is illustrated schemat-
ically in Fig. 2. The plot assumes the existence of three different time scales9 such that

tT D ≫ tf ≫ tq . (8)

It should be remembered that Ŝ is only defined on the quasi-equilibrium time scale tq . The
inequality (8) reflects the important assumption of the fluctuation theory that all states
arising during the fluctuations are quasi-equilibrium states. As such, they can always be
assigned a non-equilibrium entropy Ŝ. The assumption that the approach to equilibrium
(relaxation) and equilibrium fluctuations realize the same, quasi-equilibrium, states of the
system bears a certain analogy with the fluctuation-dissipation theorem in statistical me-
chanics [56–58].
The non-equilibrium entropy Ŝ reaches its maximum value S for some set of internal
parameters λ01 , ..., λ0m called the equilibrium internal parameters. The latter must satisfy
the necessary conditions of extremum,
!
∂ Ŝ
= 0, i = 1, ..., m, (9)
∂λi
λj6=i ,X1 ,...,Xn

and the relation

Ŝ λ01 , ..., λ0m , X1 , ..., Xn = S (X1 , ..., Xn ) , (10)
where the right-hand side is the fundamental equation of the substance. This is illustrated
schematically in Fig. 3(a) for a single internal parameter λ.
As an example, consider again the isolated cylinder with two compartments discussed
in Sec. 4 (Fig. 1). We have three internal parameters: λ1 = E1 , λ2 = V1 and λ3 = N1 .
Applying Eq.(9) to Ŝ given by Eq.(6), we have

∂ Ŝ ∂S(E1 , V1 , N1 ) ∂S(E2 , V2 , N2 ) 1 1
= − = − = 0, (11)
∂E1 ∂E1 ∂E2 T1 T2
8
The weak formulation of the second law of thermodynamics has a “negative” (prohibitive) character
stating what cannot happen - the entropy cannot decrease; whereas the strong formulation makes a “positive”
statement predicting what will happen - the system will reach equilibrium and the entropy a maximum.
9
Generally, the time scales may depend on the observable. For the present discussion, it will suffice to
represent each time scale by one “characteristic” order of magnitude.

9
 Entropy
S
S
~ tq
S ^ tf
S t TD
Quasi-equilibrium
Equilibrium
Non-equilibrium
time
Figure 2: Schematic temperature dependence of entropy of an isolated system. The system
is initially in a highly non-equilibrium state in which the non-equilibrium entropy Ŝ is
undefined. As the system evolves towards equilibrium, it eventually reaches the quasi-
equilibrium stage in which the non-equilibrium entropy Ŝ can be defined on a time scale
tq . The subsequent time dependence of Ŝ is shown by a broken line to emphasize that
individual values of Ŝ can only be measured within time intervals on the order of tq . The
values of Ŝ can never exceed the maximum (micro-canonical) entropy S. By the second
law of thermodynamics, Ŝ averaged over the thermodynamic time scale tT D (S̃, shown by
the dashed blue line) monotonically increases with time and eventually levels out when the
system reaches equilibrium. In the equilibrium state, Ŝ can randomly deviate down from
its maximum value S on the fluctuation time scale tf such that tq ≪ tf ≪ tT D , always
returning to the maximum value when the fluctuation is over.

^
S W
S Wm

λ0 λ λ0 λ
(a) (b)

Figure 3: Schematic diagram of equilibrium fluctuations in an isolated system with one

internal parameter λ. (a) Non-equilibrium entropy Ŝ reaches its maximum S at the equilib-
rium value of the internal parameter λ0 . The arrows symbolize equilibrium fluctuations of λ
on either side of λ0 accompanied by downward deviations of Ŝ from S. (b) The probability
density of λ has a sharp maximum at λ = λ0 of a height Wm .

10
 ∂ Ŝ ∂S(E1 , V1 , N1 ) ∂S(E2 , V2 , N2 ) p1 p2
= − = − = 0, (12)
∂V1 ∂V1 ∂V2 T1 T2
∂ Ŝ ∂S(E1 , V1 , N1 ) ∂S(E2 , V2 , N2 ) µ1 µ2
= − =− + = 0, (13)
∂N1 ∂N1 ∂N2 T1 T2
where we used the standard thermodynamic relations [25–27, 58] ∂S/∂E = 1/T , ∂S/∂V =
p/T and ∂S/∂N = −µ/T (T being temperature, p pressure and µ chemical potential). We
thus recover the well-known thermodynamic equilibrium conditions T1 = T2 , p1 = p2 and
µ1 = µ2 . These conditions must be satisfied for the equilibrium parameter set (E10 , V10 , N10 ).
As another example, the non-equilibrium entropy of a gas mixture with the chemical
reaction A + B ⇋ C is defined by Eq.(7) with a single internal parameter λ = ∆NA . The
necessary condition of equilibrium reads
(i) (i) (i)
∂S(E, V, NA − ∆NA , NB − ∆NA , NC + ∆NA )
∂∆NA
∂S(E, V, NA , NB , NC ) ∂S(E, V, NA , NB , NC ) ∂S(E, V, NA , NB , NC )
= − − +
∂NA ∂NB ∂NC
= 0. (14)

We obtained the standard chemical equilibrium condition µA + µB = µC , which must be

satisfied at some λ0 = ∆NA0 . The equilibrium amounts of the chemical components are then
(i) (i) (i)
NA0 = NA − ∆NA0 , NB0 = NB − ∆NA0 and NC0 = NC + ∆NA0 .

6 The general law of fluctuations

During equilibrium fluctuations, the internal parameters randomly deviate from the equi-
librium values λ01 , ..., λ0m on the time scale of tf . These fluctuations of the internal parameters
are accompanied by deviations of the non-equilibrium entropy Ŝ down from its maximum
value S [Fig. 3(a)].
The Fundamental Postulate of the thermodynamic theory of fluctuations states that the
probability density W (λ1 , ..., λm ) of internal parameters of an equilibrium isolated system is
" #
Ŝ(λ1 , ..., λm )
W (λ1, ..., λm ) = A exp , (15)
k

where A is a constant. To simplify the notations, we do not display the fixed parameters
X1 , ..., Xn as arguments of W and Ŝ.
By definition, W (λ1 , ..., λm )dλ1 ...dλm is the probability of finding the internal parameters
in the infinitesimal region dλ1 ...dλm of the parameter space. The constant A can be found
from the normalization condition
ˆ
W (λ1 , ..., λm )dλ1 ...dλm = 1, (16)

11
where the integration extends over the entire parameter space. Since Ŝ reaches its maximum
value
S = Ŝ(λ01 , ..., λ0m ) (17)
at the equilibrium parameter set λ01 , ..., λ0m [cf. Eq.(10)], W has a peak of the height
 
S
Wm = A exp (18)
k

at λ01 , ..., λ0m [Fig. 3(b)]. In other words, λ01 , ..., λ0m is the most probable parameter set. Using
this fact, Eq.(15) can be recast as
" #
Ŝ(λ1 , ..., λm ) − S
W (λ1 , ..., λm ) = Wm exp . (19)
k

Equation (19) is the central equation of the thermodynamic fluctuation theory. It is

the departure point for thermodynamic analysis of all statistical properties of equilibrium
fluctuations.
The average (expectation) value of any internal parameter λi is obtained by
" #
Ŝ(λ1 , ..., λm ) − S
ˆ
λi = Wm λi exp dλ1 ...dλm . (20)
k

Note that this average does not generally coincide with the most probable value λ0i . It
is common to use the symbol ∆ to denote deviations of fluctuating properties from their
average values. For example,
∆λi ≡ λi − λi . (21)
The covariance (second correlation moment) of any pair of internal parameters can be com-
puted by " #
Ŝ(λ1 , ..., λm ) − S
ˆ
∆λi ∆λj = Wm ∆λi ∆λj exp dλ1 ...dλm , (22)
k
while the mean-square fluctuation of λi is
" #
Ŝ(λ1 , ..., λm ) − S
ˆ
(∆λi )2 = Wm (∆λi )2 exp dλ1 ...dλm . (23)
k

7 Statistical-mechanical interpretation of fluctuations

In this section we will discuss the general law of fluctuations formulated above from the
viewpoint of statistical-mechanics.
As before, we consider an isolated system whose equilibrium state is defined by a set of
fixed parameters X1 , ..., Xn . If multiple observations are made, the system can be found in
different macro-states, each specified by a set of internal parameters λ1 , ..., λm . Let us first
consider a system with discrete λ-parameters. Each parameter set λ1 , ..., λm is implemented

12
by a certain number Ω of micro-states. The latter number is sometimes called the degeneracy
of the macro-state. The macro-states are analogs of the equilibrium and quasi-equilibrium
states introduced in the thermodynamic theory (Sec. 4). The degeneracy of a given macro-
state is a function of the internal parameters λ1 , ..., λm and the fixed parameters X1 , ..., Xn :

Ω = Ω (λ1 , ..., λm , X1 , ..., Xn ) . (24)

It is postulated that for an isolated equilibrium system, there is one unique macro-state
for which the degeneracy reaches a maximum value Ωmax . Denoting the internal parameters
corresponding to this macro-state λ01 , ..., λ0m , the maximum degeneracy depends solely on
the fixed parameters X1 , ..., Xn :


Ωmax (X1 , ..., Xn ) ≡ Ω λ01 , ..., λ0m , X1 , ..., Xn . (25)

For a given set of fixed parameters X1 , ..., Xn , the degeneracy Ω of an equilibrium isolated
system fluctuates with time from its maximum value Ωmax down to smaller values. These
equilibrium fluctuations follow the Fundamental Postulate of Statistical Mechanics stating
that the probability P of a single micro-state is the same for all micro-states and depends
only on the fixed parameters X1 , ..., Xn :

P = P (X1 , ..., Xn ) . (26)

This statement is often referred to as the micro-canonical distribution. It immediately

follows that the probability of any macro-state is simply P Ω. Thus, macro-states with a
larger degeneracy have a proportionately higher probability to occur. The macro-state with
the largest degeneracy Ωmax is the most probable one.
Analysis of fluctuations relies solely on the fact that P is the same for all micro-states
and does not require a specific knowledge of P . We note, however, that P can be found
from the normalization condition stating that the sum of probabilities of all possible macro-
states of the system be unity. In the approximation where fluctuations are neglected, only
the most probable macro-state with the highest degeneracy Ωmax is considered. Then the
normalization condition reads P Ωmax = 1, giving
1
P = . (27)
Ωmax
On the other hand, when fluctuations are significant, all possible degeneracies from Ω = 1
to Ωmax must be included. From the normalization condition
Ω
X max

PΩ = 1 (28)
Ω=1

we obtain
2
P = , (29)
Ωmax (Ωmax + 1)
which for large systems can be approximated by 2/Ω2max .

13
 Regardless of the specific value of P , the product P Ω is the probability of a given macro-
state and thus the probability
Wλ1 ,...,λm = P Ω (λ1 , ..., λm ) (30)
of the corresponding set of internal parameters. Here, as in Sec. 6, we simplify the notations
by suppressing the fixed parameters X1 , ..., Xn as arguments of P and Ω. The probabil-
ity Wλ1 ,...,λm reaches a maximum value, denoted Wm , at the most probable parameter set
λ01 , ..., λ0m , for which

Wm = P Ω λ01 , ..., λ0m = P Ωmax . (31)
Combining Eqs.(30) and (31),
Ω (λ1 , ..., λm )
Wλ1 ,...,λm = Wm . (32)
Ωmax
We now define the non-equilibrium entropy of an isolated system by
Ŝ (λ1 , ..., λm ) ≡ k ln Ω (λ1 , ..., λm ) (33)
and the equilibrium entropy by
S ≡ k ln Ωmax . (34)
It immediately follows that the non-equilibrium entropy can never exceed its equilibrium
value S. It also follows that in equilibrium, Ŝ can occasionally fluctuate down from S but
eventually returns to S as the most probable value.
The probability of any given set of internal parameters given by Eq.(32) can now be
rewritten in the form
" #
Ŝ(λ1 , ..., λm ) − S
Wλ1 ,...,λm = Wm exp . (35)
k
This equation treats λ1 , ..., λm as discrete parameters. For continuous parameters, Eq.(35)
remains valid but gives the probability density of the internal parameters. The obtained
Eq.(35) is identical to Eq.(19) of the thermodynamic theory of fluctuations (Sec. 6).

8 The Gaussian law of fluctuations

In macroscopic systems, the peak of W at the equilibrium parameter set λ01 , ..., λ0m is
extremely sharp. In most cases, one can safely replace the exponent in Eq.(19) by its Taylor
expansion at λ01 , ..., λ0m truncated at quadratic terms:
m
! m
!
2
X ∂ Ŝ 1 X ∂ Ŝ
Ŝ(λ1 , ..., λm ) − S = (λi − λ0i ) + (λi − λ0i )(λj − λ0j ), (36)
i=1
∂λ i 2 i,j=1
∂λ i ∂λ j
0 0

where the subscript 0 indicates that the derivatives are taken at λ01 , ..., λ0m . The linear terms
vanish by the extremum condition (9). As a result, Eq.(19) reduces to the multi-variable
Gaussian distribution
" m
#
1X
W (λ1, ..., λm ) = Wm exp − Λij (λi − λ0i )(λj − λ0j ) (37)
2 i,j=1

14
with the symmetrical matrix !
2
1 ∂ Ŝ
Λij ≡ − (38)
k ∂λi ∂λj
0
which we call the stability matrix of the system. For normally stable systems this ma-
trix must be positive-definite.10 Applying the normalization condition (16), the Gaussian
distribution can be rewritten as
s " #
m
|Λ| 1X
W (λ1 , ..., λm ) = exp − Λij (λi − λ0i )(λj − λ0j ) , (39)
(2π)m 2 i,j=1

where Λ is the determinant of Λij .11

We can now make use of the well-known mathematical properties of the normal distri-
bution. The average value of any parameter λi coincides with its most probable value,

λi = λ0i , i = 1, ..., m. (40)

For covariances of internal parameters we have

∆λi ∆λj = Λ−1

ij , i, j = 1, ..., m, (41)

where Λ−1ij is the inverse of the stability matrix Λij . In particular, the mean-square fluctua-
tion of a parameter λi equals

(∆λi )2 = Λ−1
ii , i = 1, ..., m. (42)

We can also introduce the following linear combinations of internal parameters:

m
X
ϕi ≡ Λij λj , i = 1, ..., m. (43)
j=1

It can be shown that

∆ϕi = 0, i = 1, ..., m, (44)
∆ϕi ∆λj = δij , i, j = 1, ..., m, (45)
∆ϕi ∆ϕj = Λij , i, j = 1, ..., m, (46)
and therefore
(∆ϕi )2 = Λii , i = 1, ..., m. (47)
Thus Λij is the covariance matrix of the parameters ϕi .
An important property of Eqs.(41), (42), (46) and (47) is that their left-hand sides
contain fluctuations, whereas the right-hand sides are computed in the state of equilibrium
10
We follow the classification [28] dividing all systems into normally stable, critically stable, unstable and
metastable. In normally stable systems, all fluctuations are finite.
11
The Gaussian distribution (39) can be interpreted as a consequence of the central limit theorem for
a system composed of many statistically independent subsystems, see for example Chapter V of Khinchin
[44].

15
and can be obtained by the standard formalism of equilibrium thermodynamics. In other
words, these equations relate fluctuations to equilibrium properties of the system.
As a simple illustration, let us revisit the isolated system divided in two compartments
(Fig. 1). Suppose the system is in equilibrium. Consider a particular type of fluctuations
in which the compartments can only exchange energy at fixed volumes and numbers of
particles. Thus, the only internal parameter is λ1 = E1 . The stability matrix (38) reduces
to the scalar !
1 ∂ 2 Ŝ
Λ=− . (48)
k ∂E12
0
From Eq.(11),
∂ Ŝ 1 1
= − (49)
∂E1 T1 T2
and thus    
1 1 ∂T1 1 ∂T2 1 1 1
Λ= + = + , (50)
k T12 ∂E1 T22 ∂E2 0 kT02 cv N1 N2
where T0 is the equilibrium temperature of both compartments and
 
1 ∂E
cv = (51)
N ∂T V,N
is the specific heat of the substance per particle at the equilibrium temperature T0 . Applying
Eq.(42) we obtain the energy fluctuation of compartment 1:
  1 N1 N2
(∆E1 )2 = = kT02 cv . (52)
V1 ,N1 Λ N1 + N2
When the second compartment is much larger than the first, N2 ≫ N1 , this equation
reduces to  
(∆E1 )2 = N1 kT02 cv . (53)
V1 ,N1

In this limit, the second compartment serves as a heat bath (reservoir) for the first. The
relative fluctuation (variance) of E1 becomes
  1/2
(∆E1 )2
V1 ,N1 T0 (cv k)1/2
= √ , (54)
E1 ε N1
where ε = E1 /N1 is the equilibrium energy per particle. This
√ relation shows that the relative
fluctuation of energy decreases with the system size as 1/ N .

9 Fluctuations in canonical systems

9.1 General relations
We now return to the general fluctuation law formulated in Secs. 6 and will apply it to a
relatively small (but still macroscopic) subsystem of an equilibrium isolated system. We will

16
 ISOLATED SYSTEM

Reservoir
r r
X ,..., X
1 n

X1,..., Xn

X1,..., Xn
Figure 4: A subsystem chosen inside an isolated system characterized by a set extensive
properties X̃1 , ..., X̃n . These properties are partitioned between the subsystem and the com-
plementary system (reservoir). The partitioning is specified by the parameter sets X1 , ..., Xn
and X1r , ..., Xnr , respectively.

adopt the terminology in which we refer to the selected subsystem as the “canonical system”,
or simply “system”, and to the remaining part of the isolated system as the “complementary
system” or “reservoir”.
An example of the situation was already given by the system with two compartments
(Fig. 1) when the second compartment is much larger than the first. We will now con-
sider a more general case in which the entire isolated system is characterized by a set of
conserved extensive properties (additive invariants) X̃1 , ..., X̃n . These properties are parti-
tioned between our system, X1 , ..., Xn , and the reservoir, X1r , ..., Xnr (Fig. 4). For some of
these properties, the partitioning is fixed once and for all, while m remaining properties can
“flow” between the system and the reservoir as a result of equilibrium fluctuations. Such
fluctuations are subject to the conservation constraints
Xi + Xir = X̃i = const, i = 1, ..., m. (55)
We assume that the exchange of the quantities Xi occurs slowly enough that we can
consider both the system and the reservoir as quasi-equilibrium systems. Accordingly, they
obey the fundamental equations
S = S (X1 , ..., Xn ) (56)
for the system and
Sr = Sr (X1r , ..., Xnr ) (57)
for the reservoir. Note that these two equations are generally different, allowing the system
and the reservoir to be composed of different substances.
In the subsequent analysis, a special role will be played by the entropy derivatives
∂S
Fi ≡ , (58)
∂Xi
called the thermodynamic “forces” conjugate to the extensive variables Xi . Similar deriva-
tives ∂Sr /∂Xir define the forces Fir in the reservoir. For example, for a simple fluid with the

17
fundamental equation S(E, V, N), the thermodynamic forces are F1 = 1/T , F2 = p/T and
F3 = −µ/T . Note that the thermodynamic forces are intensive variables.
We define the non-equilibrium entropy of the entire isolated system as
 
Ŝ(X1 , ..., Xm , Xm+1 , ..., Xn , X̃1 , ..., X̃n ) ≡ S (X1 , ..., Xn )+Sr X̃1 − X1 , ..., X̃n − Xn , (59)
| {z }
λ
where X1 , ..., Xm play the role of the fluctuating internal parameters λi , whereas the remain-
ing parameters Xm+1 , ..., Xn , X̃1 , ..., X̃n remain fixed. The equilibrium value of this entropy
is


S̃ = S X10 , ..., Xm 0
, Xm+1 , ..., Xn
 
+ Sr X̃1 − X10 , ..., X̃m − Xm 0
, X̃m+1 − Xm+1 , ..., X̃n − Xn , (60)
where X10 , ..., Xm 0
are the equilibrium values of the fluctuating parameters. To simplify
the notations, from now on we will suppress the fixed parameters in all equations, writing
Eqs.(59) and (60) in the form
 
Ŝ (X1 , ..., Xm ) ≡ S (X1 , ..., Xm ) + Sr X̃1 − X1 , ..., X̃m − Xm , (61)

 
S̃ = S X10 , ..., Xm
0
+ Sr X̃1 − X10 , ..., X̃m − Xm0
. (62)
Before analyzing fluctuations, we will apply the equilibrium conditions (9) to the internal
parameters λi = Xi (i = 1, ..., m):
!    
∂ Ŝ ∂S ∂Sr
= − r
= Fi0 − (Fir )0 = 0, i = 1, ..., m, (63)
∂λi ∂Xi 0 ∂Xi 0
0
where, as usual, the symbol 0 indicates that the derivatives are computed for the equilibrium
parameters X10 , ..., Xm0
. We thus conclude that in equilibrium, the system and the reservoir
have equal values of the thermodynamic forces conjugate to the fluctuating parameters.
We now turn to fluctuations. The probability density of the internal parameters
X1 , ..., Xm is given by Eq.(19), which takes the form
" #
Ŝ(X1 , ..., Xm ) − S̃
W (X1 , ..., Xm ) = Wm exp
k
 
S(X1 , ..., Xm ) − S(X10 , ..., Xm
0
)
= Wm exp
k
" #
Sr (X̃1 − X1 , ..., X̃m − Xm ) − Sr (X̃1 − X10 , ..., X̃m − Xm
0
)
× exp .(64)
k
In the second exponential factor, we will expand Sr around the equilibrium parameters
X10 , ..., Xm0
:
m
X
0 0
Sr (X̃1 − X1 , ..., X̃m − Xm ) − Sr (X̃1 − X1 , ..., X̃m − Xm ) = − Fi0 (Xi − Xi0 )
i=1
m  
1X ∂ 2 Sr
+ (Xi − Xi0 )(Xj − Xj0 ). (65)
2 i,j=1 ∂Xir ∂Xjr 0

18
 As will be shown later, the second derivatives appearing in the quadratic terms scale as
1/N r , where N r is the number of particles in the reservoir. Since the reservoir is assumed
to be much larger than our system, the quadratic terms in Eq.(65) can be neglected in
comparison with the linear terms. Neglecting the quadratic terms constitutes what can be
called the “reservoir approximation”, in which
m
X
Sr (X̃1 − X1 , ..., X̃m − Xm ) − Sr (X̃1 − X10 , ..., X̃m − Xm
0
)=− Fi0 (Xi − Xi0 ). (66)
i=1

Inserting this equation in Eq.(64), we obtain

 m 
X
0
 S(X1 , ..., Xm ) − S(X1 , ..., Xm ) −
0
Fi0 (Xi − Xi0 ) 
 i=1

W (X1 , ..., Xm ) = Wm exp  . (67)

 k 


Note that this equation contains only properties of the canonical system and does not depend
on properties of the reservoir. The role of the reservoir is only to impose the thermodynamic
forces Fi0 .
Equation (67) is an important relation that will serve as the starting point for deriving
all fluctuation properties of canonical systems.
In some canonical systems, there can be additional internal parameters, say Y1 , Y2 , ..., Yl ,
that can fluctuate without affecting the reservoir parameters X1r , ..., Xnr . In such cases, the
set of internal parameters is X1 , ..., Xm , Y1 , Y2 , ..., Yl and Eq.(61) is replaced by

 
Ŝ (X1 , ..., Xm , Y1 , Y2, ..., Yl ) ≡ S (X1 , ..., Xm , Y1, Y2 , ..., Yl ) + Sr X̃1 − X1 , ..., X̃m − Xm .
(68)
The equilibrium conditions with respect to the parameters decoupled from the reservoir is
!  
∂ Ŝ ∂S 0
= ≡ Fi+m = 0, (69)
∂Yi ∂Yi 0
0

where i = 1, 2, ..., l. It is easy to see that the additional parameters can be incorporated into
the parameter list X1 , ..., Xm by simply assigning them zero thermodynamic forces. Accord-
ingly, Eq.(67) remains valid, except that the sum with Fi0 excludes the terms corresponding
to the decoupled parameters Y1 , Y2 , ..., Yl .

9.2 The generalized canonical distribution

Before proceeding with applications of Eq.(67), we return for a moment to the statistical-
mechanical interpretation of fluctuations discussed in Sec. 7. For discrete internal parameters
X1 , ..., Xm , Eq.(67) gives the probability WX1 ,...,Xm of a given parameter set X1 , ..., Xm .
Recall that on the quasi-equilibrium time scale tq , our system can be treated as if it were
equilibrium and isolated with fixed values of X1 , ..., Xm . Thus, its fundamental equation

19
S(X1 , ..., Xm ) gives the equilibrium entropy S = k ln Ωmax . The probability P of a single
micro-state can be obtained from the relation WX1 ,...,Xm = P Ωmax , which gives
 
WX1 ,...,Xm S(X1 , ..., Xm )
P = = WX1 ,...,Xm exp − . (70)
Ωmax k

Combining this relation with Eq.(67) we obtain

 m 
X
0
 S(X1 , ..., Xm ) +
0
Fi0 (Xi − Xi0 ) 
 i=1

−
P = Wm exp  , (71)
 k 


which can be rewritten as !

m
1X 0
P = A exp − F Xi , (72)
k i=1 i
where  m 
X
 S(X10 , ..., Xm
0
) − Fi0 Xi0 
 i=1

−
A = Wm exp   (73)
 k 


is a constant.
Equation (72) constitutes the generalized canonical distribution describing several sta-
tistical ensembles.12
Assuming that one of the extensive parameters, say X1 , is energy, Eq.(72) can be rewrit-
ten in the form  m 
X
 E+ T0 Fi0 Xi 
 i=2

P = A exp  − . (74)

 kT0



As an illustration, suppose our system is a simple fluid. Either one or two of its extensive
parameters E, V and N can fluctuate with the third parameter fixed (otherwise the system
is undefined). For example, suppose only energy can fluctuate while V and N remain fixed.
From Eq.(74), the probability of a micro-state is
 
E
P = A exp − , fixed V, N, (75)
kT0
12
The way we arrived at Eq.(72) is not intended to be a rigorous statistical-mechanical derivation of
the canonical distribution. The goal was only to show that the thermodynamic fluctuation relation (67) is
consistent with the distribution of micro-states known from statistical mechanics.

20
a relation which is known as the canonical, or NV T , distribution. If both energy and volume
are allowed to fluctuate, we have X2 = V with F2 = p/T , giving the NpT distribution
 
E + p0 V
P = A exp − , fixed N. (76)
kT0

If the energy and the number of particles can fluctuate at a fixed volume (open system),
we have X2 = N and F2 = −µ/T . Equation (74) gives the grand-canonical, or µV T ,
distribution  
E − µ0 N
P = A exp − , fixed V. (77)
kT0
Of course, the pre-exponential factor A is different in all three distributions and is related
to the respective partition functions through the probability normalization condition.
Note that the distributions (75), (76) and (77) contain the temperature T0 , pressure
p0 and the chemical potential µ0 corresponding to the exact equilibrium with the reservoir.
These properties are fixed parameters imposed by the reservoir, by contrast to the fluctuating
temperature, pressure and chemical potential inside the system.
We re-emphasize the important difference between Eq.(72) and the previously derived
Eq.(67): Eq.(72) gives the probability of a single micro-state of the system, whereas Eq.(67)
is the probability distribution of the parameters X1 , ..., Xm characterizing different macro-
states.

9.3 The Gaussian law of canonical fluctuations

We now return to the probability distribution of the fluctuating parameters given by
Eq.(67). Let us expand the entropy of the system around its equilibrium value, keeping only
linear and quadratic terms:
m m  
X 1X ∂2S
0
S(X1 , ..., Xm ) − S(X10 , ..., Xm ) = Fi0 (Xi − Xi0 ) + (Xi − Xi0 )(Xj − Xj0 ).
i=1
2 i,j=1 ∂Xi ∂Xj 0
(78)
Inserting this expansion in Eq.(67), we arrive at the multi-variable Gaussian distribution
" m
#
1X
W (X1 , ..., Xm ) = Wm exp − Λij (Xi − Xi0 )(Xj − Xj0 ) , (79)
2 i,j=1

with the stability matrix

   
1 ∂2S 1 ∂Fi
Λij ≡ − =− . (80)
k ∂Xi ∂Xj 0 k ∂Xj 0

In this approximation, the average values X i are numerically equal to the most probable
values Xi0 . Using the standard properties of the Gaussian distribution (Sec. 8), we obtain
the fluctuation relations
∆Xi ∆Xj = Λ−1 ij , i, j = 1, ..., m, (81)

21
and thus
(∆Xi )2 = Λ−1
ii , i = 1, ..., m. (82)
We can also introduce the parameters
m
X
ϕi ≡ Λij Xj , i = 1, ..., m, (83)
j=1

for which m  
1X ∂Fi 1
ϕi − ϕi = − (Xj − Xj0 ) ≡ − ∆Fi . (84)
k j=1 ∂Xj 0 k
Here, ∆Fi is the deviation of Fi from its equilibrium value evaluated in the quadratic ap-
proximation (78). Using the properties (46) and (47) of the Gaussian distribution, we obtain
the fluctuation relations for thermodynamic forces:
∆Xi ∆Fj = −kδij , i, j = 1, ..., m, (85)
∆Fi ∆Fj = k 2 Λij , i, j = 1, ..., m, (86)
and therefore
(∆Fi )2 = k 2 Λii , i = 1, ..., m. (87)
Equations (82) and (87) show that for normally stable states, the derivatives (∂Xi /∂Fi )0
and (∂Fi /∂Xi )0 must be positive. At critical points one of these derivatives tends to infinity
and the respective fluctuation diverges (unless the reservoir has a finite size, see Sec. 10).
Another form of the foregoing equations is obtained by introducing the thermodynamic
potential
m
X
Θ(F1 , ..., Fm ) ≡ S − Fi Xi (88)
i=1

as the Legendre transform of the entropy (Massieu function) [25–27]. Then,

∂Θ
Xi = − (89)
∂Fi
and  
∂2Θ
Λ−1
ij =k . (90)
∂Fi ∂Fj 0
For normally stable states (∂ 2 Θ/∂Fi2 )0 > 0 for all i = 1, ..., m.
It is interesting to evaluate canonical fluctuations of entropy relative to its equilibrium
value S(X10, ..., Xm
0
). Using the quadratic approximation (78), we have
m m
X 1 X
∆S = Fi0 ∆Xi − k Λij ∆Xi ∆Xj . (91)
i=1
2 i,j=1

Averaging this equation over the distribution (79) and taking into account Eq.(81), we
obtain m m
X
0 1 X 1
∆S = Fi ∆Xi − k Λij ∆Xi ∆Xj = − k. (92)
i=1
2 i,j=1
2

22
Thus, in the Gaussian approximation the canonical entropy equals
1
S̄ = S(X10 , ..., Xm
0
) − k. (93)
2
Considering that entropy scales in proportion to the total of number of particles N, the
relative fluctuation ∆S/S̄ ∝ N −1 . Using the same Gaussian approximation, the mean-
square fluctuation of entropy scales in proportion to N,
m
X m
X
(∆S)2 = Fi0 Fj0 ∆Xi ∆Xj = Fi0 Fj0 Λ−1
ij ∝ N. (94)
i,j=1 i,j=1

where we neglected higher-order moments of the distribution. Thus, the relative fluctuation
 1/2
(∆S)2 /S̄ ∝ N −1/2 . Both estimates show that, with increasing size of the system,
the canonical entropy S̄ tends to the micro-canonical entropy S satisfying the fundamental
equation (4). While the two entropies are conceptually different, they converge to the same
numerical value in the thermodynamic limit.

9.4 Canonical fluctuations in the energy scheme

The treatment of fluctuations based on the fundamental equation (4) is called the entropy
scheme. In practical situations, energy is always one of the extensive parameters fluctuating
between the canonical system and the reservoir.13 Suppose X1 = E and thus F1 = 1/T .
Equation (4) can be inverted to obtain

E = E (Z1 , ..., Zn ) , (95)

where Z1 = S and Zi = Xi , i = 2, ..., n. This equation is called the fundamental equation in

the energy representation, or the energy scheme [28]. In the energy scheme, the thermody-
namic forces Pi = ∂E/∂Zi conjugate to the extensive arguments Z1 , ..., Zn are P1 = T and
Pi = −T Fi , i = 2, ..., n. For example, for a simple fluid E = E(S, V, N) and thus P1 = T ,
P2 = −p and P3 = µ.
In many applications, we are more interested in fluctuations of the P -forces rather than
F -forces. The formalism presented below expresses canonical fluctuations in terms of the
variable sets Z1 , ..., Zn and P1 , ..., Pn corresponding to the energy scheme.
As before, we consider a system coupled to a reservoir and described by n extensive
properties X1 , ..., Xn . Out of them, m properties X1 , ..., Xm can flow between the system
and the reservoir while the remaining properties Xm+1 , ..., Xn are fixed. Assuming X1 = E,
we transform the variable set X1 , ..., Xm to Z1 , ..., Zm and perform a Legendre transformation
to obtain the thermodynamic potential appropriate for the the energy scheme:
m
X
Φ(P1 , ..., Pm ) ≡ E − Pi Z i . (96)
i=1

13
While “adiabatic ensembles” in which the system is adiabatically isolated but can still exchange particles
with the reservoir are helpful theoretical constructions, their practical implementation in experiments or
atomistic simulations is unfeasible.

23
This potential has the property
∂Φ
= −Zi . (97)
∂Pi
Returning to the parameter distribution (67), we can now recast it in the form
 m 
X
 E − E0 − Pi0 (Zi − Zi0 )   
 i=1
 R
W = Wm exp −
 ≡ Wm exp − kT0 ,
 (98)
 kT0 

where we denote m
X
R ≡ E − E0 − Pi0 (Zi − Zi0 ). (99)
i=1

It can be shown that R has the meaning of the reversible work required for the creation
of the fluctuated state of the canonical system at a fixed value of the total entropy of the
system and the reservoir [24, 58].
Equation (98) is the canonical distribution function of the extensive parameters
Z1 , ..., Zm . For small fluctuations, we can expand E − E0 in these parameters and neglect
all terms beyond quadratic to obtain the multi-variable Gaussian distribution
" m
#
1X
W (Z1 , ..., Zm ) = Wm exp − Kij (Zi − Zi0 )(Zj − Zj0 ) , (100)
2 i,j=1

where the matrix    

1 ∂2E 1 ∂Pi
Kij ≡ = (101)
kT0 ∂Zi ∂Zj 0 kT0 ∂Zj 0

can be called [28] the generalized stiffness matrix. Its inverse

 
−1 ∂Zi ∂2Φ
Kij = kT0 = −kT0 (102)
∂Pj 0 ∂Pi ∂Pj

can be then called [28] the generalized compliance matrix.

Using the standard properties of the Gaussian distribution (Sec. 8), we obtain at once
the covariances and mean-square fluctuations of all extensive and intensive properties of the
canonical system in the energy scheme:

∆Zi ∆Zj = Kij−1 , i, j = 1, ..., m, (103)

(∆Zi )2 = Kii−1 , i = 1, ..., m, (104)

∆Zi ∆Pj = kT0 δij , i, j = 1, ..., m, (105)
∆Pi ∆Pj = (kT0 )2 Kij , i, j = 1, ..., m, (106)
(∆Pi )2 = (kT0 )2 Kii , i = 1, ..., m. (107)

24
 Having these fluctuation relations, it is straightforward to find the covariance of any
two thermodynamic properties y = y(Z1, ..., Zm ) and z = z(Z1 , ..., Zm ). Indeed, for small
fluctuations,
m  
X ∂y
∆y = ∆Zi , (108)
i=1
∂Z i 0
m  
X ∂z
∆z = ∆Zi , (109)
i=1
∂Zi 0
from which
m     m    
X ∂y ∂z X ∂y ∂z
∆y∆z = ∆Zi ∆Zj = Kij−1 . (110)
i,j=1
∂Z i 0 ∂Z j 0 i,j=1
∂Z i 0 ∂Z j 0

Similarly, for any two functions y = y(P1, ..., Pm ) and z = z(P1 , ..., Pm ),
m     m    
X ∂y ∂z 2
X ∂y ∂z
∆y∆z = ∆Pi ∆Pj = (kT0 ) Kij . (111)
i,j=1
∂P i 0 ∂P j 0 i,j=1
∂P i 0 ∂P j 0

The same calculation method can be applied when y is a function Z1 , ..., Zm while z is a
function of P1 , ..., Pm , or when both are functions of mixed sets of Z’s and P ’s.

9.5 Canonical fluctuations in a simple fluid

To demonstrate the foregoing formalism for computing canonical fluctuations, we will
apply it to a simple fluid coupled to a reservoir. We will consider three different ensembles,
deriving equations for the mean-square fluctuations and covariances of all thermodynamic
properties in each case. The calculations will be conducted in the energy scheme, although
the entropy scheme could be applied just as well. We are working with the extensive and
intensive parameter sets (Z1 , Z2 , Z3) = (S, V, N) and (P1 , P2 , P3 ) = (T, −p, µ), respectively.
To simplify the notations, we will suppress the index 0 but it is implied that all derivatives
are taken at the equilibrium state.
While a number of fluctuation relations for fluids can be found in the literature, in this
paper we present a complete set of such relations for each ensemble. For the NV T ensem-
ble, there are 5 fluctuating parameters (E, S, T, p, µ). Out of 25 covariances that can be
formed among them, only 15 are distinct due to the symmetry of the covariance. Similarly,
in the NpT and µV T ensembles, there are 6 fluctuating parameters (E, S, V, T, p, µ) and
(E, S, N, T, p, µ), respectively.14 Taking into account the symmetry of covariances, 21 dis-
tinct fluctuation relations exist for each of these ensembles. All these fluctuation relations
will be derived below in a systematic manner enabled by the proposed formalism.
The equilibrium properties appearing in right-hand sides of the fluctuation relations will
be presented in simplest possible form. They can be further transformed to many other
14
Recall that all intensive properties of a canonical system are allowed to fluctuate, including those which
are imposed by the reservoir. For example, in the µV T ensemble, the reservoir imposes a certain temperature
T0 and chemical potential µ0 ; however, the actual temperature and chemical potential inside the system
fluctuate around the imposed values.

25
forms expressing them through specific heats, moduli or other experimentally accessible
properties. We do not pursue such transformations as they lie outside the fluctuation topic
and depend on the problem at hand.

9.5.1 NV T ensemble
Suppose the volume and number of particles in the system are fixed while the energy and
entropy can fluctuate (canonical ensemble). The stiffness and compliance matrices reduce
to the scalars  
1 ∂T 1
K11 = = (112)
kT0 ∂S V,N Nkcv
and
−1
K11 = Nkcv . (113)
By Eq.(104),
(∆S)2 = K11
−1
= Nkcv . (114)
Applying Eq.(110), we find the energy fluctuation
 2
∂E
(∆E)2 = −1
K11 = NkT02 cv (115)
∂S V,N

and the energy-entropy covariance

 
∂E −1
∆E∆S = K11 = NkT0 cv , (116)
∂S V,N

where we used the relation (∂E/∂S)V,N = T0 . For temperature fluctuations, we follow

Eq.(107) to obtain
kT02
(∆T )2 = (kT0 )2 K11 = . (117)
Ncv
We next derive fluctuations of the chemical potential and pressure. Both are dependent
variables that can be expressed as functions of temperature and the number density of
particles, ρ ≡ N/V . The latter is fixed, leaving only T as an independent argument.
Applying Eq.(111) we obtain
 2  2
2 ∂µ kT02 ∂µ
(∆µ)2 = (kT0 ) K11 = , (118)
∂T ρ Ncv ∂T ρ

2  2
2 ∂p kT02 ∂p
(∆p)2 = (kT0 ) K11 = , (119)
∂T ρ Ncv ∂T ρ
       
2 ∂µ ∂p kT02 ∂µ ∂p
∆µ∆p = (kT0 ) K11 = , (120)
∂T ρ ∂T ρ Ncv ∂T ρ ∂T ρ
   
2 ∂µ kT02 ∂µ
∆µ∆T = (kT0 ) K11 = , (121)
∂T ρ Ncv ∂T ρ

26
    
2 ∂p kT02 ∂p
∆p∆T = (kT0 ) K11 = . (122)
∂T ρ Ncv ∂T ρ

For “cross-fluctuations” between extensive and intensive properties, we first apply

Eq.(105) to obtain
∆S∆T = kT0 . (123)
The remaining “cross-fluctuations” are easily computed by the method outlined by Eqs.(108)-
(111):  
∂E
∆E∆T = ∆S∆T = kT02 , (124)
∂S V,N
   
∂p ∂p
∆S∆p = ∆S∆T = kT0 , (125)
∂T ρ ∂T ρ
     
∂E ∂p 2 ∂p
∆E∆p = ∆S∆T = kT0 , (126)
∂S V,N ∂T ρ ∂T ρ
   
∂µ ∂µ
∆S∆µ = ∆S∆T = kT0 , (127)
∂T ρ ∂T ρ
     
∂E ∂µ 2 ∂µ
∆E∆µ = ∆S∆T = kT0 . (128)
∂S V,N ∂T ρ ∂T ρ
Table 1 summarizes all fluctuation relation derived for the NV T ensemble.

9.5.2 NpT ensemble

Now suppose that the volume of the system can also fluctuate while N =const. The
stiffness matrix is now 2 × 2 and its K11 component is still given by Eq.(112). For the
remaining components we have
   
∂S ∂p
 
1 ∂T 1 ∂V ∂T ρ
K12 = =−  T,N = − , (129)
kT0 ∂V S,N kT0 ∂S Nkcv
∂T V,N

where we used the Maxwell relation

   
∂S ∂p
= , (130)
∂V T,N ∂T V,N

and    2  
1 ∂p 1 ∂p 1 ∂p
K22 =− = − , (131)
kT0 ∂V S,N Nkcv ∂T ρ kT0 ∂V T,N

where we used the identity

     2
∂p ∂p T ∂p
= − . (132)
∂V S,N ∂V T,N Ncv ∂T V,N

27
    
kT02 ∂µ ∂p
(∆E)2 = NkT02 cv ∆µ∆p =
Ncv ∂T ρ ∂T ρ

(∆S)2 = Nkcv ∆S∆T = kT0

(∆E)(∆S) = NkT0 cv ∆E∆T = kT02

 
2
kT02 ∂p
(∆T ) = ∆S∆p = kT0
Ncv ∂T ρ

 2  
kT02 ∂µ ∂p
(∆µ)2 = ∆E∆p = kT02
Ncv ∂T ρ ∂T ρ

 2  
kT02 ∂p ∂µ
(∆p)2 = ∆S∆µ = kT0
Ncv ∂T ρ ∂T ρ

   
kT02 ∂p ∂µ
∆T ∆p = ∆E∆µ = kT02
Ncv ∂T ρ ∂T ρ

 
kT02 ∂µ
∆µ∆T =
Ncv ∂T ρ

Table 1: Complete set of fluctuation relations for a simple fluid in the NV T ensemble.

28
Note that the derivatives (∂p/∂V )T,N and (∂p/∂V )S,N are related to the isothermal and
adiabatic moduli, respectively.
Thus, the stiffness matrix of the fluid takes the form
   
∂p
1 −

 ∂T ρ 

1  
K=      (133)
Nkcv   2  
 ∂p ∂p Ncv ∂p 
− −
∂T ρ ∂T ρ T0 ∂V T,N

with the compliance matrix

 2     
∂p Ncv ∂p ∂p
 − 

∂V
  ∂T ρ T0 ∂V T,N ∂T ρ
 
K −1 = −kT0  . (134)
∂p T,N    
 ∂p 
1
∂T ρ

We are ready to calculate the fluctuations. To find the entropy and volume fluctuations,
we apply Eq.(103):
   2
−1 ∂V ∂p
(∆S)2 = K11 = Nkcv − kT0 ≡ Nkcp , (135)
∂p T,N ∂T ρ
 
∂V−1
(∆V = )2
= −kT0 K22 , (136)
∂p T,N
     
−1 ∂V ∂p ∂V
∆S∆V = K12 = −kT0 = kT0 , (137)
∂p T,N ∂T ρ ∂T p,N
where cp is the specific heat at constant pressure. Fluctuation relations involving energy are
calculated by Eq.(110):
 2  2   
∂E −1 ∂E ∂E ∂E−1 −1
(∆E)2 = +K11 +2 K22 K12
∂S
V,N S,N ∂V ∂S V,N ∂V S,N
 "   #2
∂V ∂p p0
= NkT02 cv − kT03 − . (138)
∂p T,N ∂T ρ T0

   
∂E −1 ∂E −1
∆E∆S = K11 + K12
∂S V,N ∂V S,N
    "  #
∂V ∂p ∂p p0
= NkT0 cv − kT02 − . (139)
∂p T,N ∂T ρ ∂T ρ T0

29
    
∂E ∂E
−1 −1
∆E∆V = + K12 K22
∂S
V,N ∂V S,N
  "  #
∂V ∂p p 0
= −kT02 − . (140)
∂p T,N ∂T ρ T0

Turning to fluctuations of intensive parameters, we use Eqs.(106) and (107) to obtain

kT02
(∆T )2 = (kT0 )2 K11 = , (141)
Ncv
"   2 #
∂p T0 ∂p
(∆p)2 = (kT0 )2 K22 = −kT0 − , (142)
∂V T,N Ncv ∂T ρ
 
kT02 ∂p
∆T ∆p = −(kT0 )2 K12 = . (143)
Ncv ∂T ρ
To derive fluctuation relations involving the chemical potential, we treat it as a function of
T and p. Using Eq.(111),
"   2     #
2
∂µ ∂µ ∂µ ∂µ
(∆µ)2 = (kT0 )2 K11 + K22 − 2 K12
∂T p ∂p T ∂T p ∂p T
 
= (kT0 )2 s20 K11 + v02 K22 + 2s0 v0 K12
"   #2  
kT02 ∂p 2 ∂p
= v0 − s0 − kT0 v0 , (144)
Ncv ∂T ρ ∂V T,N

where  
∂µ
s0 = − (145)
∂T p

is the equilibrium entropy per particle and

 
∂µ
v0 = (146)
∂p p

is the equilibrium volume per particle (v0 = 1/ρ0 ). Similarly,

"    # "   #
2
∂µ ∂µ kT 0 ∂p
∆µ∆T = (kT0 )2 K11 − K12 = v0 − s0 , (147)
∂T p ∂p T Ncv ∂T ρ

"     #
∂µ ∂µ
∆µ∆p = (kT0 )2 − K12 + K22 (148)
∂T p ∂p T
    "   #
∂p kT02 ∂p ∂p
= −kT0 v0 + v0 − s0 . (149)
∂V T,N Ncv ∂T ρ ∂T ρ

30
 Equation (105) gives the “cross-fluctuations”

∆S∆T = −∆V ∆p = kT0 , (150)

∆S∆p = 0, (151)
∆V ∆T = 0. (152)
The remaining “cross-fluctuations” are computed in a manner similar to Eqs.(108)-(111):
 
∂µ
∆S∆µ = ∆S∆T = −kT0 s0 , (153)
∂T p
 
∂E
∆E∆T = ∆S∆T = kT02 , (154)
∂S V,N
 
∂E
∆E∆p = ∆V ∆p = kT0 p0 , (155)
∂V S,N
       
∂E ∂µ ∂E ∂µ
∆E∆µ = ∆S∆T + ∆V ∆p = kT0 (p0 v0 − T0 s0 ) ,
∂S V,N ∂T p ∂V S,N ∂p T
  (156)
∂µ
∆V ∆µ = ∆V ∆p = −kT0 v0 . (157)
∂p T
The fluctuation relations for the NpT ensemble are summarized in Table 2.

9.5.3 µV T ensemble
We now turn to the grand-canonical ensemble, in which the volume of the system is fixed
while the energy and the number of particles can fluctuate. The equilibrium temperature
T0 and the chemical potential µ0 are fixed by the reservoir. The independent fluctuating
variables are S and N, with the conjugate thermodynamic forces P1 = T , and P2 = µ.
Energy is a dependent extensive parameter. The chemical potential and pressure are also
dependent parameters treated as functions of T and ρ.
All fluctuation relations can be obtained from those for the NpT ensemble by simply
swapping the variables V → N, p → −µ. In particular, the derivative (∂p/∂T )ρ is replaced
by −(∂µ/∂T )ρ and (∂p/∂V )T,N by −V −1 (∂µ/∂ρ)T . Furthermore, it easy to verify that s0
is replaced by s0 /v0 and v0 by 1/v0 . The stiffness matrix becomes
   
∂µ
 1 ∂T ρ 
1  


K=      (158)
N0 kcv   2 
 ∂µ ∂µ ρ0 cv ∂µ 
+
∂T ρ ∂T ρ T0 ∂ρ T

31
    2  
∂V ∂p ∂p
(∆S)2 = Nkcv − kT0 ∆µ∆p = −kT0 v0
∂p T,N ∂T ρ ∂V T,N
2
  "   #
kT0 ∂p ∂p
+ v0 − s0
Ncv ∂T ρ ∂T ρ
 
∂V
∆S∆V = kT0 ∆S∆T = kT0
∂T p,N

 
∂V
(∆V )2 = −kT0 ∆V ∆p = −kT0
∂p T,N

(∆E)2 = NkT02 cv ∆S∆p = 0

  "  #2
∂V ∂p p 0
−kT03 −
∂p T,N ∂T ρ T0

∆E∆S = NkT0 cv ∆V ∆T = 0
    "  #
∂V ∂p ∂p p 0
−kT02 −
∂p T,N ∂T ρ ∂T ρ T0

  "  #
∂V ∂p p0
∆E∆V = −kT02 − ∆S∆µ = −kT0 s0
∂p T,N ∂T ρ T0

kT02
(∆T )2 = ∆E∆T = kT02
Ncv
"   2 #
∂p T0 ∂p
(∆p)2 = −kT0 − ∆E∆p = kT0 p0
∂V T,N Ncv ∂T ρ

 
kT02 ∂p
∆T ∆p = ∆E∆µ = kT0 (p0 v0 − T0 s0 )
Ncv ∂T ρ

"   #2  
2
kT0 ∂p 2 ∂p
(∆µ)2 = v0 − s0 − kT0 v0 ∆V ∆µ = −kT0 v0
Ncv ∂T ρ ∂V T,N
"   #
kT02 ∂p
∆µ∆T = v0 − s0
Ncv ∂T ρ

Table 2: Complete set of fluctuation relations for a simple fluid in the NpT ensemble.

32
  −1  2  
∂µ ∂µ kT0 ∂µ
(∆S)2 = ρ0 kcv V + kT0 V ∆µ∆p =
∂ρ T ∂T ρ ρ0 v0 V ∂ρ
2
  " T   #
kT0 ∂µ ∂µ
+ s0 +
ρ0 cv v0 V ∂T ρ ∂T ρ
 −1  
∂µ ∂µ
∆S∆N = −kT0 V ∆S∆T = kT0
∂ρ T ∂T ρ

 −1
∂µ
(∆N)2 = kT0 V ∆N∆µ = kT0
∂ρ T

(∆E)2 = ρ0 kT02 cv V ∆S∆µ = 0

 −1 "  #2
∂µ ∂µ µ 0
+kT03 − V
∂ρ T ∂T ρ T0

∆E∆S = ρ0 kT0 cv V ∆N∆T = 0

 −1   "  #
∂µ ∂µ ∂µ µ 0
−kT02 − V
∂ρ T ∂T ρ ∂T ρ T0

 −1 "  #
∂µ ∂µ µ0 s0
∆E∆N = −kT02 − V ∆S∆p = kT0
∂ρ T ∂T ρ T0 v0

kT02
(∆T )2 = ∆E∆T = kT02
ρ0 cv V
"    2 #
kT 0 ∂µ T0 ∂µ
(∆µ)2 = ρ0 + ∆E∆µ = kT0 µ0
ρ0 V ∂ρ T cv ∂T ρ
 
kT02 ∂µ kT0
∆T ∆µ = ∆E∆p = (µ0 + T0 s0 )
ρ0 cv V ∂T ρ v0

  "   #2
kT0 ∂µ kT02 ∂µ kT0
(∆p)2 = 2 + s0 + ∆N∆p =
v0 V ∂ρ T ρ0 cv v02 V ∂T ρ v0
"  #
kT02 ∂µ
∆p∆T = s0 +
ρ0 cv v0 V ∂T ρ

Table 3: Complete set of fluctuation relations for a simple fluid in the µV T ensemble.

33
with the inverse
 2    
∂µ ρ0 cv ∂µ ∂µ
 −1  + − 
∂µ  ∂T ρ T0 ∂ρ T ∂T ρ
 
K −1 = kT0 V  . (159)
∂ρ T    
 ∂µ 
− 1
∂T ρ

The results of the calculations are summarized in Table 3.

The following scaling properties of fluctuations are evident from Tables 1, 2 and 3:

• Covariances of extensive properties scale in proportion to the system size (N or V ).

• Covariances of intensive properties scale as inverse of the system size (N −1 or V −1 ).

• “Cross-fluctuations” between extensive and intensive properties are independent of the

system size.

The first of these properties explains why we could neglect the quadratic terms in the
expansion for the reservoir entropy, justifying the reservoir approximation (66) discussed in
Sec. 9.1.

10 Finite-reservoir ensembles
10.1 General considerations
We now revisit the system coupled to a reservoir discussed in Sec. 9.1 (Fig. 4). The
system and the reservoir are described by extensive parameters X1 , ..., Xn and X1r , ..., Xnr ,
respectively, with fixed total amounts X̃i = Xi + Xir . Suppose m of these parameters
are allowed to flow back and forth between the system and the reservoir, whereas the re-
maining (n − m) parameters are fixed. Both the system and the reservoir are assumed
to be single-phase substances with different fundamental equations S = S(X1 , ..., Xn ) and
Sr = Sr (X1r , ..., Xnr ), respectively. The probability distribution of the fluctuating parameters
of the system is given by Eq.(64), which is repeated here for convenience:
 
S(X1 , ..., Xm ) − S(X10 , ..., Xm
0
)
W (X1 , ..., Xm ) = Wm exp
k
" #
Sr (X̃1 − X1 , ..., X̃m − Xm ) − Sr (X̃1 − X10 , ..., X̃m − Xm
0
)
× exp (160)
.
k

In Sec. 9.1, the entropy of the reservoir was approximated by the linear expansion (66)
at the equilibrium parameter set X10 , ..., Xm
0
. This linear approximation was justified by the
large size of the reservoir in comparison with the system. We will now lift this assumption

34
and adopt a more general expansion, Eq.(65), which includes quadratic terms:
m
X
Sr (X̃1 − X1 , ..., X̃m − Xm ) − Sr (X̃1 − X10 , ..., X̃m − 0
Xm ) = − Fi0 (Xi − Xi0 )
i=1
m
1 X r
− k Λ (Xi − Xi0 )(Xj − Xj0 ), (161)
2 i,j=1 ij

where  
1 ∂ 2 Sr
Λrij ≡− (162)
k ∂Xir ∂Xjr 0
is the stability matrix of the reservoir. As usual, the derivatives are taken at the equilibrium
parameter set X10 , ..., Xm 0
.
Before analyzing fluctuations, let us examine the micro-state probability distribu-
tion P (X1 , ..., Xm ). Recall that the latter is obtained from W (X1 , ..., Xm ) by dropping
S(X1 , ..., Xm ) in the first line of Eq.(160) [see Eq.(70) in Sec. 9.2]. Thus,
" #
−S(X10 , ..., Xm
0
) + Sr (X̃1 − X1 , ..., X̃m − Xm ) − Sr (X̃1 − X10 , ..., X̃m − Xm0
)
P = Wm exp
k
" m m
#
1X 0 1X r
= A exp − F Xi − Λ (Xi − Xi0 )(Xj − Xj0 ) , (163)
k i=1 i 2 i,j=1 ij

where A is a constant that can be found from the probability normalization condition.
Thus, in addition to the usual linear terms Fi0 Xi appearing in the generalized canonical
distribution (72), we now have a quadratic form representing the effect of the finite reservoir
on the micro-state probability of the system.
The matrix elements Λrij scale as 1/Nr . When the size of the reservoir increases to
infinity, the quadratic terms in Eq.(163) vanish and we return to the generalized canonical
distribution (72). In the opposite limit when the reservoir becomes extremely small in
comparison with our system, the probability distribution is dominated by the quadratic
terms. The distribution becomes a very sharp Gaussian peak. Accordingly, all micro-
states whose Xi are even slightly deviated from Xi0 become extremely improbable. The
system behaves as if isolated and Eq.(163) approaches the micro-canonical distribution. In
other words, Eq.(163) interpolates between the generalized canonical and micro-canonical
distributions; we can smoothly transition from one distribution to the other by varying the
reservoir size and thus scaling the reservoir stability matrix (162).
We now return to the parameter distribution function W (X1 , ..., Xm ) and expand the
entropy of our system in a manner similar to Eq.(161):
m m
X 1 X
S(X1 , ..., Xm ) − S(X10 , ..., Xm
0
) = Fi0 (Xi − Xi0 ) − k Λij (Xi − Xi0 )(Xj − Xj0 ), (164)
i=1
2 i,j=1

with the usual stability matrix

 
1 ∂2S
Λij ≡ − . (165)
k ∂Xi ∂Xj 0

35
Inserting the expansions (161) and (164) in Eq.(160), we obtain the multi-variable Gaussian
distribution
" m
#
1 X
W (X1 , ..., Xm ) = Wm exp − Λ̃ij (Xi − Xi0 )(Xj − Xj0 ) , (166)
2 i,j=1

with the stability matrix

Λ̃ij = Λij + Λrij . (167)
We thus arrive at the remarkable result that the finite size of the reservoir does not
change the Gaussian form of the parameter distribution of the system. However, the reser-
voir does contribute to the matrix elements of the Gaussian distribution in an additive
manner. We can, therefore, easily compute the mean-square fluctuations and covariances of
the parameters X1 , ..., Xm using the known Eqs.(81) and (82) with Λij replaced by the full
stability matrix Λ̃ij :
∆Xi ∆Xj = Λ̃−1
ij , i, j = 1, ..., m, (168)
(∆Xi )2 = Λ̃−1
ii , i = 1, ..., m. (169)
A number of other interesting relations can be derived. It follows from Eqs.(164) and
(161) that
m
∂S 0
X
Fi = = Fi − k Λij (Xj − Xj0 ), (170)
∂Xi j=1
m
∂Sr X
Fir = = Fi
0
+ k Λrij (Xj − Xj0 ). (171)
∂Xir j=1

Averaging these equations over the distribution and using the fact that for a Gaussian
distribution Xj = Xj0 , we obtain
Fir = Fi = Fi0 . (172)
In other words, not only the most probable but also average values of the thermodynamic
forces are the same in the system and in the reservoir. Furthermore, subtracting Eq.(170)
from Eq.(171) we have
m
X
r
Fi − Fi = k Λ̃il (Xl − Xl0 ). (173)
l=1

Starting from this equation, we can generate a set of fluctuation relations involving the
difference (Fir − Fi ). For example, multiplying Eq.(173) by ∆Xj and averaging over the
parameter distribution, we have

(Fir − Fi )∆Xj = kδij , i, j = 1, ..., m, (174)

where we used Eq.(168). Likewise, multiplying Eq.(173) by (Fjr − Fj ) and averaging over
the distribution, we obtain

(Fir − Fi )(Fjr − Fj ) = k 2 Λ̃ij , i, j = 1, ..., m, (175)

36
where we used Eq.(174). In particular,

(Fir − Fi )2 = k 2 Λ̃ii , i = 1, ..., m. (176)

These equations describe fluctuations of the differences between the values of the intensive
parameters inside our system and in the reservoir.
We can also compute fluctuations inside our system. To this end, we rewrite Eq.(170)
as m
X
∆Fi = −k Λip ∆Xp . (177)
p=1

Multiplying this equation by ∆Xj and averaging over the distribution, we obtain the co-
variances m
X
∆Fi ∆Xj = −k Λip Λ̃−1
pj , i, j = 1, ..., m. (178)
p=1

Similarly, multiplying Eq.(177) by ∆Fj and averaging over the distribution, we have
m
X
∆Fi ∆Fj = k 2 Λip Λpq Λ̃−1
qj , i, j = 1, ..., m, (179)
p,q=1

m
X
2
(∆Fi )2 =k Λip Λpq Λ̃−1
qi , i = 1, ..., m. (180)
p,q=1

As expected, in the limit of an infinitely large reservoir when Λ̃ij becomes identical to Λij ,
Eqs.(178), (179) and (180) reduce to Eqs.(85) (86) and (87) for the generalized canonical
ensemble.

10.2 Application to a simple fluid

To demonstrate the finite-ensemble formalism, consider a fixed volume V of a simple
fluid embedded in another (generally, different) simple fluid serving as a reservoir. This
situation was already discussed in Sec. 9.5.3 in the context of the grand-canonical ensemble.
This time, however, the reservoir will be treated as an infinite source of heat (thermostat)
but a finite source of particles. We thus have two fluctuating extensive properties, X1 = E
and X2 = N, with the conjugate thermodynamic forces F1 = 1/T and F2 = −µ/T .
Because the reservoir has an infinite heat capacity, we have Λr11 = 0. Assuming for
simplicity that the chemical potential of the reservoir µr is temperature-independent, we
additionally have Λr12 = Λr21 = 0. Thus, the reservoir stability matrix is
 
r 0 0
Λij = (181)
0 Λr22
with a single nonzero parameter Λr22 .
The micro-state probability distribution for our system can be found from Eq.(163):
 
E − µ0 N 1 r 2
P = A exp − − Λ22 (N − N0 ) . (182)
kT0 2

37
As a side note, a similar equation can be obtained for a closed system coupled to a thermostat
with a finite heat capacity. In this case, the role of the number of particles is played by the
energy, so that the additional term contributed by the thermostat is quadratic in (E − E0 ):
 
E 1 r 2
P = A exp − − Λ (E − E0 ) . (183)
kT0 2 11

This situation is known in the literature as the Gaussian ensemble [47, 50]. By analogy,
formula (182) can also be referred to as the Gaussian ensemble distribution.
Equation (182) can be rewritten the form
 
E − µ̂N
P = A exp − , (184)
kT0

where
µ̂ ≡ a − bN (185)
with
a = µ0 + Λr22 kT0 N0 , (186)
1
b = Λr22 kT0 . (187)
2
Equation (184) looks similar to the grand-canonical distribution, except that the fixed chem-
ical potential of the reservoir µ0 is replaced an “effective” chemical potential µ̂ linearly depen-
dent on the current number of particles N. In the limit of Λr22 → 0 (infinitely large reservoir),
this dependence disappears and we return to the standard grand-canonical distribution with
the chemical potential µ0 .
Equations (184)-(187) can be implemented in Monte Carlo computer simulations us-
ing an algorithm similar to that of the grand-canonical ensemble [1] but with a chemical
potential adjustable “on the fly”. For binary mixtures, this approach was implemented in
the semi-grand canonical mode as “variance-constrained” Monte Carlo [54, 55]. A simi-
lar semi-grand canonical Monte Carlo method, called feedback Monte Carlo, was recently
employed for interface free energy calculations by the capillary fluctuation approach [21].
Related Monte Carlo algorithms were earlier developed for closed systems coupled to a finite-
capacity thermostat. The initially proposed “dynamical ensemble” method [59] developed
for such systems later evolved into more sophisticated and computationally more efficient
iterative Monte Carlo schemes based on Eq.(183) [52, 53].
We now return to the analysis of fluctuations. To simplify the calculations, we will neglect
the temperature dependence of the chemical potential, treating it solely as a function of the
particle density ρ. Following the calculation methods presented in Sec. 9.5, it is easy to
obtain the stability matrix of the fluid in the entropy scheme,
 
1 −µ0
1  
Λij =     (188)
ρ0 cv kT02 V  ∂µ 2

−µ0 ρ0 cv T0 + µ0
∂ρ T

38
and its inverse
   
 −1 ρ0 cv T0 ∂µ + µ20 µ0 
∂µ  ∂ρ
Λ−1
ij = kT0 V  T . (189)
∂ρ T  
µ0 1
Accordingly,
 
1 −µ0
1  
Λ̃ij =     (190)
ρ0 cv kT02 V  ∂µ 
−µ0 ρ0 cv T0 + µ20 + ρ0 cv kT02 Λr22 V
∂ρ T

and
   
∂µ
ρ0 cv T0 ∂ρ + µ20 + ρ0 cv kT02 Λr22 V µ0 
kT0 V
Λ̃−1
ij =    T . (191)
∂µ  
+ kT0 V Λr22 µ0 1
∂ρ T

We can now compute the mean-square fluctuations of the energy and the number of
particles:
 
2 ∂µ
ρ0 cv kT0 V + kT0 V (µ20 + ρ0 cv kT02 Λr22 V )
∂ρ T
(∆E)2 = Λ̃−1
11 =   , (192)
∂µ r
+ kT0 V Λ22
∂ρ T
kT0 V
(∆N)2 = Λ̃−1
22 =   . (193)
∂µ
+ kT0 V Λr22
∂ρ T

Similarly, for the covariance ∆E∆N we have

kT µ V
∆E∆N = Λ̃−1
12 =   0 0 . (194)
∂µ r
+ kT0 V Λ22
∂ρ T
Fluctuations of the chemical potential could be found directly from Eq.(180). Instead,
we will take a shortcut by taking advantage of the fact that µ = µ(ρ) and therefore
 
∂µ ∆N
∆µ = . (195)
∂ρ T V
Taking a square of this equation and averaging over the distribution, we obtain
 2
kT0 ∂µ
 2
2
∂µ (∆N)2 V ∂ρ T
(∆µ) = =   , (196)
∂ρ T V 2 ∂µ r
+ kT0 V Λ22
∂ρ T

39
where we used Eq.(193). Applying the same method,
 
∂µ
  kT0
(∆N)2
∂µ T ∂ρ
∆µ∆N = =  , (197)
T ∂ρ
V ∂µ r
+ kT0 V Λ22
∂ρ T
 
∂µ
  µ0 kT0
∂µ ∆E∆N ∂ρ T
∆µ∆E = =  . (198)
∂ρ T V ∂µ r
+ kT0 V Λ22
∂ρ T
Continuing along this line, all other covariances can be readily calculated.
Consider limiting cases of the above relations. When Λr22 → 0 (the reservoir is an infinite
source of particles), we exactly recover all results for the grand-canonical ensemble listed
in Table 3. When Λr22 → ∞ (the reservoir cannot supply or absorb particles), the system
behaves as virtually closed (NV T ). Accordingly, the energy fluctuation (192) reduces to

(∆E)2 = ρ0 kT02 cv V (199)

in agreement with the NV T result (Table 1). All other mean-square fluctuations and co-
variances appearing in Eqs.(193)-(198) tend to zero as 1/Λr22 . As evident from Table 1, in
the NV T ensemble they are proportional to (∂µ/∂T )ρ and thus generally not zero. In the
present calculations, they vanished because we neglected the temperature dependence of µ
for the sake of simplicity. More accurate but rather tedious calculations taking into account
the temperature dependence of µ exactly recover the fluctuation relations from Table 1 when
Λr22 → ∞.
It is interesting to note that at the limit of normal thermodynamic stability when
(∂µ/∂ρ)T → 0, the canonical fluctuation (∆N)2 diverges to infinity (Table 3). However,
when the system is coupled to a finite-size reservoir with a positive coefficient Λr22 > 0, the
latter stabilizes the fluid and its compositional fluctuations remain finite and equal to 1/Λr22 ,
see Eq.(193). The finite size of the reservoir smooths the singularity.

11 Fluctuations of grain boundary properties

GBs are interfaces separating regions of the same crystalline solid phase (called grains)
with different crystallographic orientations. GBs can strongly impact physical and me-
chanical properties of crystalline materials, especially in nano-structured systems. While
thermodynamics and kinetics of GBs has been studied extensively for several decades [60],
it is only recently that fluctuations of GB properties have become the subject of dedicated
research, primarily by atomistic computer simulations [15–22].
In preparation for the discussion of GB fluctuations, we will give a briefly account of GB
thermodynamics focusing on a plane GB in a binary substitutional solid solution.

40
 GB

Figure 5: Two grains with different crystallographic orientations separated by a plane grain
boundary (GB). The shaded areas bounded by dashed lines mark a region containing the
GB and a comparison region inside one of the grains containing the same total number of
atoms. These two regions are selected for calculations of GB excess quantities.

11.1 Grain boundary thermodynamics

Consider a bicrystal (i.e., a system of two grains) enclosed in a rectangular box with
rigid impermeable walls (Fig. 5). Since the system is closed, the total number of particles
N is fixed. We will assume that the grains possess appropriate crystallographic symmetries
such that the mechanical stresses existing in the bicrystal do not create a difference in strain
energy densities in the grans that would cause GB migration.
Before discussing GB thermodynamics, we will first specify thermodynamic properties
of the homogeneous single-crystalline solid solution forming the grains. We postulate the
following fundamental equation of this solution in the energy representation:15
Eg = Eg (Sg , Vg , Nsg , N2g ). (200)
Here, the suffix g indicates that the properties refer to a region in one of the grains, such
as the shaded layer shown in Fig. 5, N2g is the number of atoms of species 2 and Nsg is the
number of lattice sites in the region. Vacancies are neglected, so that the number of atoms
of species 1 is N1g = Nsg − N2g . The differential of Eq.(200) is
dEg = T dSg − pdVg + ϕdNsg + MdN2g , (201)
where T = ∂Eg /∂Sg is temperature, p = −∂Eg /∂Vg is pressure (negative of the stress
component normal to the GB plane), M = ∂Eg /∂N2g is the diffusion potential of species 2
relative to species 1 [63–65], and ϕ = ∂Eg /∂Nsg is a thermodynamic potential conjugate to
the number of sites Nsg . On the other hand, the energy defined by Eq.(200) is a homogeneous
function of first degree in all four arguments. Using the Euler theorem for homogeneous
functions [66] we have
Eg = T Eg − pVg + ϕNsg + MN2g . (202)
15
Generally, the fundamental equation of a stressed solid depends on the deformation gradient relative
to a chosen reference state [24, 61, 62]. In the particular case considered here, the cross-section of the
solid is fixed and only deformation in the normal direction is permitted. Under such conditions, only the
normal component of the deformation gradient is relevant and varies as a linear function of the volume. It
is convenient to represent this component by the volume itself.

41
Combining Eqs.(201) and (202) we obtain the Gibbs-Duhem equation

− Sg dT + Vg dp − Nsg dϕ − N2g dM = 0. (203)

To develop GB thermodynamics, we define the GB excess of any extensive property X

by
X̃ ≡ X − Xg , (204)
where X is the value of the property for a bicrystalline layer containing the GB and Xg is the
value of the same property for a homogeneous grain layer containing the same total number
of atoms. Examples of such layers are shown by shaded regions in Fig. 5. The choice of these
layers is totally arbitrary and does not affect X̃ as long as their boundaries are unaffected
by the GB or the walls. According to this definition of excess, a bicrystalline layer can be
thought of as composed of two subsystems: (1) the reference grain layer containing the same
total number of atoms and (2) the excess system representing the GB and described by the
excess properties S̃, Ẽ, Ṽ , etc.
Just as the grain thermodynamics has been formulated above starting from a postulated
fundamental equation (200), we will postulate a fundamental equation of the GB in the
form
Ẽ = Ẽ(S̃, Ṽ , Ñ2 , A). (205)
Note that Ñ1 is not an independent variable because

Ñ1 + Ñ2 = 0 (206)

by our definition of excesses. The excess energy (205) is a homogeneous function of first
degree in all four arguments. Applying the Euler theorem [66] we have

∂ Ẽ ∂ Ẽ ∂ Ẽ ∂ Ẽ
Ẽ = S̃ + Ṽ + Ñ2 + A. (207)
∂ S̃ ∂ Ṽ ∂ Ñ2 ∂A

On the other hand, the differential of Eq.(205) is

∂ Ẽ ∂ Ẽ ∂ Ẽ ∂ Ẽ
dẼ = dS̃ + dṼ + dÑ2 + dA. (208)
∂ S̃ ∂ Ṽ ∂ Ñ2 ∂A

Combining Eqs.(207) and (208) we obtain the GB version of the Gibbs-Duhem equation:

∂ Ẽ ∂ Ẽ ∂ Ẽ ∂ Ẽ
S̃d + Ṽ d + Ñ2 d + Ad = 0. (209)
∂ S̃ ∂ Ṽ ∂ Ñ2 ∂A

Next, we find the conditions of thermodynamic equilibrium between the GB and the
grains. To this end, we consider a variation of the total energy (Ẽ + Eg ) of a bicrystalline
layer at ax fixed composition, volume and entropy. In equilibrium, this variation (dẼ + dEg )
must be zero [24]. Here, dEg is given by Eq.(201) and dẼ by Eq.(208). The variation is
subject to the following constraints:

dÑ1 + dN1g = 0, (210)

42
 dÑ2 + dN2g = 0, (211)
dṼ + dVg = 0, (212)
dS̃ + dSg = 0, (213)
as well as dA = 0. Furthermore, by adding Eqs.(210) and (211) and taking into account
Eq.(206), it follows that the number of sites in the reference grain region is conserved:
dNsg = dN1g + dN2g = 0. Inserting these constraints into the total energy variation, we arrive
at the equilibrium condition
! ! !
∂ Ẽ ∂ Ẽ ∂ Ẽ
− T dS̃ + + p dṼ + − M dÑ2 = 0. (214)
∂ S̃ ∂ Ṽ ∂ Ñ2

Because the variations dS̃, dṼ and dÑ2 are independent and can be positive or negative,
the coefficients multiplying them must be zero. We thus obtain the equilibrium values of
the derivatives of Ẽ: !
∂ Ẽ
= T0 , (215)
∂ S̃ 0
!
∂ Ẽ
= M0 , (216)
∂ Ñ2 0
!
∂ Ẽ
= −p0 . (217)
∂ Ṽ 0
As usual, the subscript 0 labels the equilibrium state.
The remaining derivative !
∂ Ẽ
γ0 ≡ (218)
∂A
0
is called the equilibrium GB free energy. The precise thermodynamic meaning of this prop-
erty will be discussed later.

11.2 Quasi-equilibrium grain boundary states

So far, we have only specified the derivatives of Ẽ with respect to S̃, Ñ2 , Ṽ and A
for equilibrium states of the system. Namely, such derivatives are equal to the equilibrium
temperature, diffusion potential, negative of pressure in the grains and the equilibrium
GB free energy, respectively. For the analysis of GB fluctuations, we must consider these
derivatives without assuming equilibrium between the GB and the grains. Such derivatives
form a new set of variables T , p, M and γ defined as follows:

∂ Ẽ
≡ T, (219)
∂ S̃
∂ Ẽ
≡ M, (220)
∂ Ñ2

43
 ∂ Ẽ
≡ −p, (221)
∂ Ṽ
∂ Ẽ
≡ γ. (222)
∂A
These new variables can be interpreted as, respectively, the local temperature, negative
of local pressure, local diffusion potential and the GB free energy when the GB is not in
equilibrium with the grains. Indeed, by postulating the fundamental equation (205) we
implicitly assumed that the GB follows this equation even in the absence of equilibrium
with the grains. One can think of such GB states as being disconnected from the actual
grains and equilibrated with imaginary grains with a different set of intensive parameters
T , M and p. In terms of the variables (219)-(222), the differential form of the fundamental
equation (205) becomes
dẼ = T dS̃ − pdṼ + MdÑ2 + γdA. (223)
The existence of quasi-equilibrium states in which the GB obeys a fundamental equation
without being in equilibrium with the grains is an important supposition of the interface
fluctuation theory. As discussed in Secs. 4 and 5, the entire thermodynamic theory of
equilibrium fluctuations is built on the assumption that all fluctuated states are states of
quasi-equilibrium. The theory discussed here is an application of this principle to GBs.
Suppose the bicrystal has been initially equilibrated and then the GB deviates from this
initial equilibrium to a nearby quasi-equilibrium state where it is no longer in equilibrium
with the grains. Denoting the equilibrium excess entropy, volume and number of atoms S̃0 ,
Ṽ0 and Ñ20 , respectively, the Gibss-Duhem equation (209) for this variations becomes

∂ Ẽ ∂ Ẽ ∂ Ẽ ∂ Ẽ
S̃0 d + Ṽ0 d + Ñ20 d + Ad
∂ S̃ ∂ Ṽ ∂ Ñ2 ∂A
0
= S̃0 dT − Ṽ0 dp + Ñ2 dM + Adγ = 0. (224)

This equation can be rewritten in the form

S̃0 Ṽ0 Ñ 0
dγ = − dT + dp − 2 dM, (225)
A A A
which will be used later.
In the particular case when the states of the grains also vary in such a manner that the
entire system undergoes an equilibrium process, Eq.(225) becomes

S̃0 Ṽ0 Ñ 0
dγ0 = − dT0 + dp0 − 2 dM0 . (226)
A A A
This is one of the forms of the Gibbs adsorption equation for a fixed GB area. This equi-
librium equation should not be confused with Eq.(225) describing fluctuations in which the
GB deviates from equilibrium with the grains.
We can now understand the meaning of the parameter γ. Let us rewrite Eq.(207) in the
form
Ẽ = T S̃ − pṼ + M Ñ2 + γA, (227)

44
from which
Ẽ − T S̃ + pṼ − M Ñ2
γ= . (228)
A
As was shown by Gibbs [24], an expression of the form E −T S +pV −Σi µi Ni is the reversible
work required for the creation of a new system with extensive properties E, S, V and Ni
by drawing the energy and matter from an infinite reservoir with a temperature T , pressure
p and chemical potentials µi . Thus, the numerator in Eq.(228) can be interpreted as the
reversible work of GB formation from an imaginary infinite reservoir with the intensive
parameters T , p and M. In the present case, the term µ1 Ñ1 + µ2 Ñ2 reduces to M Ñ2 due
to the constraint imposed by Eq.(206). In short, γ is the reversible work of GB formation
per unit area. Accordingly, the equilibrium GB free energy γ0 is the reversible work of GB
formation per unit area from a reservoir with the same intensive parameters T0 , p0 and M0
as in the actual grains:
Ẽ − T0 S̃ + p0 Ṽ − M0 Ñ2
γ0 = . (229)
A
In other words, this is the reversible work of GB formation between two infinitely large
grains.

11.3 Grain boundary fluctuations

We are now in a position to analyze GB fluctuations. We will compute generalized
canonical fluctuations of GB properties treating the grains as an infinite reservoir cou-
pled to the GB. Although the foregoing discussion utilized the energy scheme, the latter
is impractical for GB fluctuations.16 Instead, we will now switch to the entropy scheme in
which the independent additive invariants fluctuating between the GB and the grains are
X1 = Ẽ, X2 = Ṽ and X3 = Ñ2 . The conjugate thermodynamic forces are readily found
from Eq.(223): F1 = 1/T , F2 = p/T and F3 = −M/T . We can also rewrite Eq.(225) in
terms of the new variables:
X3
dγ = Bi dFi , (230)
i=1

with the coefficients

T0 (Ẽ0 − γ0 A)
B1 ≡ , (231)
A
T0 Ṽ0
B2 ≡ , (232)
A
T0 Ñ20
B3 ≡ − . (233)
A
16
The energy scheme involves fluctuations of the excess entropy S̃, which is not accessible by experiments
or atomistic simulations. In the entropy scheme, S̃ is replaced by the more accessible excess energy Ẽ.

45
 In the Gaussian approximation to fluctuations, the distribution of the excess parameters
Xi is given by Eq.(79) with the stability matrix
 ! ! !
∂ 2 S̃ ∂ 2 S̃ ∂ 2 S̃
 
 ∂ Ẽ 2 0 ∂ Ẽ∂ Ṽ 0 ∂ Ẽ∂ Ñ2 0 
 
 
 ! ! !
1 ∂ 2 S̃ ∂ 2 S̃ ∂ 2 S̃ 

Λij = −  . (234)
k  ∂ Ẽ∂ Ṽ ∂ Ṽ 2 ∂ Ṽ ∂ Ñ2 
 0 0 0
 
 ! ! ! 
∂ 2 S̃ ∂ 2 S̃ ∂ 2 S̃
 
 
2
∂ Ẽ∂ Ñ2 0 ∂ Ṽ ∂ Ñ2 0 ∂ Ñ2 0

The elements of this matrix represent equilibrium excess properties of the GB. For example,
!   
∂ 2 S̃ ∂ 1 1
= =− 2 , (235)
∂ Ẽ 2 0 ∂ Ẽ T 0 T0 Ac̃v

where !
1 ∂ Ẽ
c̃v ≡ (236)
A ∂T
0

is the GB heat capacity (more accurately, the excess heat capacity per unit GB area com-
puted at fixed Ṽ /A and Ñ2 /A). Similarly,
!  
∂ 2 S̃ ∂ p 1 p0
= =− − 2 , (237)
∂ Ṽ 20
∂ Ṽ T 0 T0 Ṽ0 κ̃T T0 Ṽ0 α̃

where !
1 ∂ Ṽ
κ̃T ≡ − (238)
Ṽ0 ∂p
0
is the isothermal compressibility of the GB and
!
1 ∂ Ṽ
α̃ ≡ (239)
Ṽ0 ∂T
0

is the GB thermal expansion coefficient.

Equilibrium GB properties can be computed directly by atomistic simulations. For
example, κ̃T can be found by computing the excess GB volume as a function of pressure
at a fixed temperature followed by numerical differentiation [67]. However, analysis of
GB fluctuations offers a more efficient approach in which the entire set of equilibrium GB
properties can be extracted from a single simulation run using fluctuation relations. Indeed,
the excess properties Ẽ, Ṽ and Ñ2 are readily accessible by atomistic simulations. By
computing their covariances, the inverse stability matrix can be found from the fluctuation
relation

46
 ∆Xi ∆Xj = Λ−1
ij , i, j = 1, 2, 3. (240)
The latter can be then inverted to obtain the stability matrix Λij and thus the equilibrium
GB properties such as c̃v , κ̃T , α̃, etc. Furthermore, knowing the stability matrix, fluctuations
of thermodynamic forces can be predicted from the fluctuation relations

∆Fi ∆Fj = k 2 Λij , i, j = 1, ..., 3. (241)

Knowing these covariances, straightforward algebraic manipulations can be applied to find

(∆M)2 , (∆T )2 and (∆p)2 , as well as the covariances ∆M∆T , ∆M∆p and ∆p∆T . As a
general rule,
3    
X
2 ∂Pi ∂Pj
∆Pi ∆Pj = k Λnm , i, j = 1, ..., 3, (242)
m,n=1
∂Fm 0 ∂Fn 0

where the transformation matrix ∂Pi /∂Fm is

 2 
−T0 p0 T0 −M0 T0
   
∂Pi  
= 0 −T0 0 . (243)
∂Fm 0  


0 0 −T0

To our knowledge, the proposed calculation scheme has not been implemented.
As an interesting application, one can predict the equilibrium mean-square fluctuation
(∆γ)2 of the GB free energy, where ∆γ ≡ γ − γ0 . The fluctuation form of Eq.(230) is
3
X
∆γ = Bi ∆Fi . (244)
i=1

Taking a square of this equation and averaging over the distribution, we obtain
3
X
(∆γ)2 = Bi Bj ∆Fi ∆Fj . (245)
i,j=1

Applying Eq.(241) for ∆Fi ∆Fj , this equation finally becomes

3
X
2
(∆γ)2 =k Bi Bj Λij . (246)
i,j=1

This calculation requires the knowledge of Λij and the equilibrium values of the excess GB
energy, volume and segregation as ingredients for the parameters Bi . The GB free energy
is usually treated as a static property. An estimate of its fluctuation (∆γ)2 would provide
a critical assessment of this approximation.

47
 w

Liquid

GB
Figure 6: Sharp interface model of GB pre-melting. The pre-melted GB is modeled by a
liquid layer of a width w bounded by two sharp interfaces endowed with thermodynamic
properties of actual solid-liquid interfaces between bulk phases.

12 Fluctuations in pre-melted grain boundaries

In this section we address another fluctuation topic associated with GBs. At high temper-
atures approaching the melting point of the solid, many GBs become atomically disordered
and can develop a liquid-like structure reminiscent of a thin liquid film [60]. Our goal is
to describe equilibrium thermodynamic properties of this film and derive the distribution
function of its width. The analysis is based on the sharp interface model of pre-melting
(Fig. 6). The grains are formed by the same binary solid solution phase as discussed in
Sec. 11.

12.1 Thermodynamic background

12.1.1 Equilibrium state of pre-melted grain boundary
We will consider the same bicrystal in a rigid box as discussed in Sec. 11, but this
time we will adopt a particular model of GB structure, namely, a liquid layer of a certain
width w bounded by two solid-liquid interfaces (Fig. 6). We start by describing equilibrium
thermodynamic properties of this boundary.
As in Sec. 11, thermodynamic properties of the solid solution forming the grains are
fully described by the fundamental equation (200) in the energy representation. The binary
liquid solution (index L) is described by its own fundamental equation

EL = EL (SL , VL , N1L , N2L ) (247)

with the differential form

dEL = TL dSL − pL dVL + µ1 dN1L + µ2 dN2L . (248)

Here, TL and pL are the temperature and pressure of the liquid and µ1 and µ2 are the
chemical potentials of the components.
To describe thermodynamics of the solid-liquid interfaces bounding the liquid layer, we
will first take a detour and consider a single plane interface between the bulk solid and liquid

48
phases. We can choose an arbitrary geometric dividing surface partitioning our system in
two parts with volumes Vg and VL . We then compute the excess of any extensive property
X by
X̃ = X − Xg − XL , (249)
where X is the value of the property for the solid-liquid system with a volume (Vg + VL )
and Xg and XL are properties of the phases computed as if the volumes Vg and VL remained
homogeneous all the way to the dividing surface [24]. By this definition, Ṽ = 0. Because
the densities of all properties are generally different on either side of the interface, the
excess defined by Eq.(249) depends on the exact position of the dividing surface within the
interface region. We will adopt the rule by which we always place the dividing surface in
such a manner that the excess of the total number of atoms N vanishes: Ñ = 0. The
fundamental equation of the solid-liquid interface is formulated in the form
Ẽ = Ẽ(S̃, Ñ2 , A), (250)
with the differential form
∂ Ẽ ∂ Ẽ ∂ Ẽ
dẼ = dS̃ + dÑ2 + dA. (251)
∂ S̃ ∂ Ñ2 ∂A
Returning to the pre-melted GB modeled by a thin layer of the liquid phase (Fig. 6), we
will assume that the solid-liquid interfaces bounding this layer follow the same fundamental
equation (250) as if they were interfaces between two bulk phases. This is, of course, an
approximations. In reality, the liquid layer width is on the order of a nanometer and its
thermodynamic properties in this narrow confinement are different from properties of bulk
liquid. In addition, the solid-liquid interfaces separated by such a short distance can interact
with each other due to their finite thickness comparable with the width of the layer. The
difference between the actual GB structure and this highly idealized liquid-layer structure
is accounted for by introducing an interaction between the two solid-liquid interfaces called
the “disjoining” interaction. The disjoining interaction (index d) is included in the model by
adding an extra energy term Ed and an associated entropy Sd , and postulating the following
fundamental equation of the disjoining interaction:
Ed = Ed (Sd , VL ). (252)
The dependence on VL has been added to make the disjoining interaction a function of the
distance between the solid-liquid interfaces, which for a fixed cross-section considered here,
is proportional to VL . Both Ed and Sd are proportional to the GB area and have the same
dimensions as the excess quantities Ẽ and S̃, respectively. Both Ed and Sd tend to zero if
the width of the liquid layer increases to infinity (bulk solid-liquid system).
To summarize the model, the GB is represented by a homogeneous liquid layer separated
from the grains by two geometrically sharp solid-liquid interfaces to which we assign ther-
modynamic properties of actual solid-liquid interfaces between bulk phases (Fig. 6). These
sharp interfaces coincide with the dividing surfaces separated by a distance w and interact
with each other by disjoining forces that follow the fundamental equation (252).
We will next derive the conditions of thermodynamic equilibrium by considering a vari-
ation of state of a bicrystalline layer indicated in Fig. 6 by a shaded region with dashed

49
boundaries. During this variation, the layer remains isolated from the rest of the system
and keeps a constant value of its entropy. Under these constraints, the variation of the total
energy
Ê = Eg + EL + Ẽ + Ed (253)
must be zero [24]. The terms of this equation are represented by the fundamental equations
(200), (247), (250) and (252), respectively. Accordingly,

dÊ = T dSg − pg dVg + ϕdNsg + MdN2g

+ TL dSL − pL dVL + µ1 dN1L + µ2 dN2L
∂ Ẽ ∂ Ẽ ∂ Ẽ
+ dS̃ + dÑ2 + dA
∂ S̃ ∂ Ñ2 ∂A
∂Ed ∂Ed
+ dSd + dVL = 0 (254)
∂Sd ∂VL
subject to the following constraints:

dSg + dSL + dS̃ + dSd = 0, (255)

dVg + dVL = 0, (256)

dN1g + dN1L + dÑ1 = 0, (257)
dN2g + dN2L + dÑ2 = 0, (258)
dÑ1 + dÑ2 = 0, (259)
dN1g + dN2g − dNsg = 0, (260)
as well as dA = 0. Here, −pg is the component of the mechanical stress in the grains normal
to the GB.
Out of 14 differentials appearing in the right-hand sides of the Eqs.(254)-(260), 7 are
eliminated by the constraints, leaving 7 independent variations:
!  
∂ Ẽ ∂Ed
dÊ = (TL − T ) dSL + − T dS̃ + − T dSd
∂ S̃ ∂Sd
 
∂Ed
+ − pL + pg dVL + (ϕ − µ1 ) dN1g
∂VL
!
∂ Ẽ
+ (ϕ − µ2 + M) dN2g + + µ1 − µ2 dÑ2 = 0. (261)
∂ Ñ2

Equating the differential coefficients to zero, we obtain the following equilibrium conditions:
!  
∂ Ẽ ∂Ed
T = TL = = ≡ T0 , (262)
∂ S̃ ∂Sd 0
0
 
∂Ed
p0L = p0g + , (263)
∂VL 0

50
 !
∂ Ẽ
= M0 = µ02 − µ01 , (264)
∂ Ñ2 0

ϕ0 = µ01 , (265)
where T0 is the uniform temperature throughout the equilibrium system. As usual, the
superscript 0 marks the equilibrium state.
Equation (263) is the mechanical equilibrium condition, which can be written in the form

p0L = p0g + p0d , (266)

where p0d is the equilibrium value of the disjoining pressure defined by

∂Ed
pd ≡ . (267)
∂VL
Thus, the liquid phase within the GB layer is subject to an additional pressure pd on top of
the pressure pg existing in the grains.

12.1.2 Calculation of the disjoining pressure

Calculation of the equilibrium disjoining pressure p0d can be simplified by reformulation
of the equilibrium conditions in the density form. The fundamental equation (200) of the
solid solution can be rewritten in the density form

εg = εg (sg , vg , cg ), (268)

where εg , sg and vg are, respectively, the energy, entropy and volume per lattice site (which
is same as per atom since we neglect vacancies), and cg is the site (atomic) fraction of species
2. Accordingly, T = ∂εg /∂sg , pg = −∂εg /∂vg and M = ∂εg /∂cg . From Eq.(202) we have

ϕ = εg − T sg + pg vg − Mcg , (269)

while the Gibbs-Duhem equation (203) becomes

sg dT − vg dpg + dϕ + cg dM = 0. (270)

It is more convenient to express all properties of the grain as functions of the variable
set (T, pg , cg ). The thermodynamic potential appropriate for these variables is

ω(T, pg , cg ) ≡ ϕ + Mcg = εg − Tg sg + pg vg . (271)

Its differential form is

dω = −sg dT + vg dpg + Mdcg , (272)
which is easily derivable from Eq.(270). Similarly, for the liquid solution we will use the
variable set (TL , pL , cL ) and, accordingly, the Gibbs free energy

g(TL , pL , cL ) ≡ εL − TL sL + pL vL = (1 − cL )µ1 + cL µ2 (273)

51
with the differential
dg = −sL dTL + vL dpL + (µ2 − µ1 )dcL , (274)
where all densities are counted per atom.
We now reformulate the equilibrium conditions (264) and (265) in the form

M(T0 , p0g , c0g ) = µ2 (T0 , p0L , c0L ) − µ1 (T0 , p0L , c0L ), (275)

ω(T0 , p0g , c0g ) = g(T0 , p0L , c0L ) + M(T0 , p0g , c0g )(c0g − c0L ), (276)
respectively. This is a system of two non-linear equations with respect to five equilibrium
parameters (T0 , p0g , p0L , c0g , c0L ). It follows that our solid-liquid system has three degrees of
freedom. Let us choose the independent variable set (T0 , p0g , c0g ). For any given state of the
grains defined by these parameters, we can solve Eqs.(275) and (276) to find the composition
c0L and pressure p0L in the liquid layer, and thus the disjoining pressure p0d = p0L − p0g .
For practical applications, we will replace this exact calculation scheme by a simple
approximation. At given temperature T0 and pressure p0g , let c∗g and c∗L be the solid and
liquid compositions corresponding to the bulk solid-liquid coexistence (p0d = 0). These two
compositions must satisfy the equations17

M(T0 , p0g , c∗g ) = µ2 (T0 , p0g , c∗L ) − µ1 (T0 , p0g , c∗L ), (277)

ω(T0 , p0g , c∗g ) = g(T0 , p0g , c∗L ) + M(T0 , p0g , c∗g )(c∗g − c∗L ). (278)
Note that the pressure continuity condition p0L = p0g valid for a two-phase bulk system with
a plane interface eliminates one degree of freedom, leaving our system with two degrees of
freedom in compliance with the Gibbs phase rule. For a fixed pressure p0g , Eqs.(277) and
(278) define the solidus and liquidus lines, c∗g (T0 ) and c∗L (T0 ), on the phase diagram of the
system.
We seek an approximate form of Eqs.(275) and (276) in the vicinity of the solidus line.
As a measure of proximity of the grains to the solidus line at given T0 and p0g , we choose
the diffusion potential deviation

∆M ≡ M(T0 , p0g , c0g ) − M(T0 , p0g , c∗g ) (279)

from its solidus value M(T0 , p0g , c∗g ). All functions appearing in Eqs.(275) and (276) can be
then approximated by linear expansions in the small parameters ∆M, (c0g − c∗g ), (c0L − c∗L )
and (p0L − p0g ). The linearized form of these equations becomes
   
∂(µ2 − µ1 ) 0 ∂(µ2 − µ1 )
∆M = ∗
(cL − cL ) + p0d , (280)
∂cL ∗ ∂p L ∗

     
∂ω ∂g ∂g
(c0g − c∗g ) = (c0L − c∗L ) + p0d
∂cg ∗ ∂cL ∗ ∂pL ∗
+ M(T0 , p0g , c∗g )(c0g − c0L − c∗g + c∗L ) + (c∗g − c∗L )∆M, (281)
17
We will not go into details, but Eqs.(277) and (278) can be given a geometric interpretation of a common
tangent construction, M being the slope of the common tangent to the curves ω(T0 , p0g , c) and g(T0 , p0g , c)
as functions of composition c at fixed T0 and p0g .

52
where the asterisk indicates that the derivatives are taken in the bulk-coexistence state
(∆M = 0). All zero-order terms have canceled out by Eqs.(277) and (278). The equations
obtained are simplified by utilizing the relations
   
∂ω ∂g
= = M(T0 , p0g , c∗g ), (282)
∂cg ∗ ∂cL ∗
 
∂g
= vL∗ , (283)
∂pL ∗
 
∂(µ2 − µ1 )
= v̄L2∗ − v̄L1∗ , (284)
∂pL ∗

where v̄Li = ∂µi /∂pL are partial molar volumes of the components in the liquid phase, with
the asterisk indicating that these volumes refer to the bulk-coexistence state. Using these
relations, Eqs.(280) and (281) become, respectively,
 
∂(µ2 − µ1 )
∆M = (c0L − c∗L ) + (v̄L2∗ − v̄L1∗ )p0d , (285)
∂cL ∗

0 = vL∗ p0d + (c∗g − c∗L )∆M. (286)

Equation (286) immediately gives us the equilibrium disjoining pressure,

c∗L − c∗g
p0d = ∆M, (287)
vL∗

while Eq.(285) can be solved for the deviation of the liquid-layer composition c0L from the
liquidus composition c∗L :

vL2∗ − vL1∗ ∗
2∗ 1∗ 0 1 − (cL − c∗g )
∆M − (vL − vL )p v ∗
c0L − c∗L =   d
=  L  ∆M. (288)
∂(µ2 − µ1 ) ∂(µ2 − µ1 )
∂cL ∗ ∂cL ∗

In many systems, the difference between the partial molar volumes of the components is
small in comparison with vL∗ and the latter equation can be simplified to
1
c0L − c∗L =   ∆M. (289)
∂(µ2 − µ1 )
∂cL ∗

For a single-component system, ∆M loses its significance. The appropriate measure

of proximity of the system to bulk melting is then the difference (T0 − Tm ) between the
temperature of the bicrystal and the bulk melting point Tm at a given pressure p0g . The
equilibrium conditions (275) and (276) reduce to one equation,

ω(T0 , p0g , 0) = g(T0 , p0L , 0), (290)

53
which for bulk equilibrium becomes

ω(Tm , p0g , 0) = g(Tm , p0g , 0). (291)

Linearizing Eq.(290) in the small parameters (T0 − Tm ) and p0d we have

     
∂ω ∂g ∂g
(T0 − Tm ) = (T0 − Tm ) + p0 , (292)
∂T ∗ ∂T ∗ ∂pL ∗ d

from which
s∗L − s∗g Hm
p0d = (T0 − Tm ) = (T0 − Tm ), (293)
vL
∗
Tm vL∗
where Hm = (s∗L − s∗g )Tm is the enthalpy of melting per atom.

12.2 Analysis of fluctuations

To describe fluctuations in the pre-melted GB, we treat the grains as an infinite reservoir
and apply the canonical fluctuation relation (98). The system coupled to the reservoir is
the GB possessing the energy
EGB ≡ EL + Ẽ + Ed . (294)
We will initially consider the full set of fluctuating parameters, which will later be reduced
to a smaller set. Namely, we will initially treat all arguments of the fundamental equations
(247), (250) and (252), except for the fixed GB area, as the fluctuating parameters. In the
energy scheme, the fluctuating parameter set is

Z = (SL , S̃, Sd , VL , N1L , N2L , Ñ2 ). (295)

The thermodynamic forces conjugate to these parameters are

∂ Ẽ ∂Ed ∂Ed ∂ Ẽ
P = (TL , , , −pL + , µ1 , µ2 , ). (296)
∂ S̃ ∂Sd ∂VL ∂ Ñ2
The reversible work of GB formation (99) takes the form
7
X
0
R = EGB − EGB − Pi0 (Zi − Zi0 ), (297)
i=1

where the thermodynamic forces are taken in the state of equilibrium. Using the equilibrium
conditions (262)-(264), we obtain

R = (EL − EL0 ) + (Ẽ − Ẽ0 ) + (Ed − Ed0 )

−T0 (SL − SL0 ) − T0 (S̃ − S̃0 ) − T0 (Sd − Sd0 )
+p0L (VL − VL0 ) − p0d (VL − VL0 )
−µ01 (N1L − N1L0 ) − µ02 (N2L − N2L0 ) − M0 (Ñ2 − Ñ20 ). (298)

54
 Next, we will make some approximations that will reduce the number of fluctuating
parameters. The difference (EL − EL0 ) appearing in Eq.(298) will be approximated by its
linear expansion around equilibrium:

EL − EL0 = T0 (SL − SL0 ) − p0L (VL − VL0 ) + µ01 (N1L − N1L0 ) + µ02 (N2L − N2L0 ). (299)

Likewise, (Ẽ − Ẽ0 ) will also be evaluated in the linear approximation:

Ẽ − Ẽ0 = T0 (S̃ − S̃0 ) + M0 (Ñ2 − Ñ20 ). (300)

The terms related to the disjoining interaction will remain intact. Inserting the linear
expansions (299) and (300) in Eq.(298), several terms cancel out and we are left with

R = (Ed − Ed0 ) − T0 (Sd − Sd0 ) − p0d (VL − VL0 ). (301)

Although Eq.(301) contains variations of three variables, only one fluctuating parameter
is independent. To see this, we apply the Legendre transformation with respect to the
disjoining entropy to obtain the thermodynamic potential

Ψ(T, VL ) = Ed (Sd , VL ) − T Sd . (302)

This thermodynamic potential is called the disjoining potential and its differential form is

dΨ = −Sd dT + pd dVL . (303)

The latter equation suggests another definition of the disjoining pressure:

 
∂Ψ
pd = . (304)
∂VL T

Thus, Eq.(301) can be rewritten in the form

R = Ψ(T0 , VL ) − Ψ(T0 , VL0 ) − p0d (VL − VL0 ), (305)

showing that, at a given temperature T0 , R is solely a function of VL . In other words, VL is

the only independent fluctuating parameter. By neglecting higher-order terms in Eqs.(299)
and (300), we have suppressed fluctuations of all other independent parameters.
The probability distribution of the liquid volume VL is given by Eq.(98), which becomes

   
R Ψ(T0 , VL ) − Ψ(T0 , VL0 ) − p0d (VL − VL0 )
W (VL ) = Wm exp − = Wm exp − . (306)
kT0 kT0

Remembering that the cross-section of the bicrystal is fixed, the finally arrive at the distri-
bution function of the liquid-layer width w:
 
ψ(T0 , w) − ψ(T0 , w0 ) − p0d (w − w0 )
W (w) = Wm exp −A , (307)
kT0

55
where w0 = VL0 /A is the equilibrium width of the layer and ψ = Ψ/A is the disjoining
potential per unit area. Note that this distribution depends on the GB area A: the larger
the area, the sharper is the peak of the distribution around w0 . The disjoining pressure can
be rewritten as pd = (∂ψ/∂w)T . Its equilibrium value p0d appearing in Eq.(307) can be found
from the bulk thermodynamic properties using Eq.(287) for a binary alloy and Eq.(293) for
a single-component system. Knowing p0d and the distribution of w at a given temperature
T0 , one can invert Eq.(307) and extract the disjoining potential ψ(T0 , w).
The GB width distribution (307) in conjunction with Eqs.(287) and (293) for the equi-
librium disjoining pressure answers the goal of this section. For single-component systems,
an equation similar to (307) was derived in previous work [16–18, 68] and utilized for calcu-
lations of disjoining interactions in Ni GBs by atomistic simulations [16–18]. In the present
paper, this distribution has been extended to binary systems.
As a further approximation, we can expand ψ(T0 , w) around the equilibrium width keep-
ing only linear and quadratic terms:
 
0 1 ∂2ψ
ψ(T0 , w) = ψ(T0 , w0 ) + pd (w − w0 ) + (w − w0 )2 , (308)
2 ∂w 2 0
where the second derivative is taken at equilibrium. This results in the Gaussian distribution
  2  
A ∂ ψ 2
W (w) = Wm exp − (w − w0 ) , (309)
2kT0 ∂w 2 0
from which the mean-square fluctuation of the GB width is
kT0
(w − w0 )2 = . (310)
A(∂ ψ/∂w 2 )
2
0

This estimate shows that √ the root-mean-square fluctuation of the liquid-layer width scales
with the GB area as 1/ A.
According to the foregoing derivation, the GB width w appearing in the sharp-interface
model of GB premelting (Fig. 6) has the meaning of the distance between two dividing
surfaces of the solid-liquid interfaces. Such dividing surfaces are characterized by a zero
excess of the total number of particles [Ṽ = 0, Ñ = 0, see excess definition (249)]. It
should be, therefore, understood that the alternate definitions of w adopted in atomistic
simulations [16–18] constitute approximations to the strict definition.

13 Concluding remarks
The thermodynamic theory of fluctuations presented here expands the postulational ba-
sis of the traditional classical thermodynamics by incorporating equilibrium fluctuations.
The main additional element of the expanded theory is the fluctuation postulate expressed
by Eq.(15) [or equivalently, Eq.(19)]. Due to its statistical nature, this postulate signifies a
break with the determinism of classical thermodynamics. The statistical component intro-
duced by this postulate propagates through the entire logical structure of thermodynamics.
The equilibrium state of a system is no longer a state but a distribution over states. Ther-
modynamic properties become distributions around the average values which are measured

56
by macroscopic experiments. Similar expansions of classical thermodynamics were previ-
ously discussed by Tisza and Quay [39] and Callen [26, 27, 38] and can be traced back to
Einstein [40]. The fluctuation postulate formulated in this paper describes fluctuations in an
isolated system and stems from the microcanonical distribution in statistical mechanics. In
this sense, our approach is different from that of Tisza and Quay [39] and Callen [26, 27, 38]
and is more aligned with Einstein’s proposal [40].
The fluctuation postulate relies on the concept of non-equilibrium entropy. Without
a proper definition of non-equilibrium entropy, we cannot talk about its increase during
equilibration or behavior during fluctuations. In some sense, non-equilibrium entropy is
defined already in classical thermodynamics. The thermodynamic equilibrium conditions of
classical thermodynamics are formulated through the entropy maximum principle [24, 26–
28]. In this principle, the actual equilibrium state of the system is compared with possible
virtual states in which the system is imagined to be partitioned into isolated equilibrium
subsystems; their entropies are then summed up and compared with the equilibrium entropy
of the actual system. The expanded thermodynamics goes one step further: it assumes
that such states are actually implemented during spontaneous fluctuations away from the
equilibrium state. The term “quasi-equilibrium” introduced in this paper captures three
defining properties of such fluctuated states: (1) non-equilibrium, (2) real (as opposed to
virtual states of classical thermodynamics), and (3) composed of subsystems that can be
treated as isolated (on a certain time scale) and each described by a fundamental equation.
The quasi-equilibrium states are characterized by a set of internal parameters λi . The
fluctuation postulate specifies the distribution function of these parameters for an isolated
system.
The second law of thermodynamics is re-formulated to include fluctuations. In the refined
formulation, it is only the entropy S̄ averaged over the fluctuations that monotonically in-
creases and reaches a maximum at equilibrium. On shorter time scales, the non-equilibrium
entropy Ŝ incessantly fluctuates around its average value and can never accede a certain
maximum value S (Figs. 2 and 3). Thus, the average entropy S̄ obtained by macroscopic
measurements (performed on the time scale tT D ) is always somewhat smaller than S.
In classical thermodynamics, a central role is played by the fundamental equation which
encapsulates all equilibrium thermodynamic properties of the system. This equation, first
discovered by Gibbs [24], looks exactly like our Eq.(4) but has a subtly different meaning.
In the classical fundamental equation, the entropy S is a static property defined through
the maximum-entropy principle (e.g., Callen’s Postulate II [26, 27] or Tisza’s Postulate Pb1
[28]) and measured experimentally. In the expanded theory, S appearing in Eq.(4) has
the statistical-mechanical meaning of the maximum micro-canonical entropy of an isolated
system. As noted in the previous paragraph, this entropy is slightly larger than the measured
entropy S̄. In other words, the measured entropy S̄ does not exactly satisfy the fundamental
equation (4) of the present theory. This fact, of course, does not create any problems.
It only reflects the possibility of different definitions of the fundamental equation for a
fluctuating system. Since both S and S̄ are well-defined for any given set of conserved
parameters X1 , ..., Xn , S̄(X1 , ..., Xn ) could also be taken as the fundamental equation [26, 27,
37, 38]. We would then have to re-define the temperature, pressure and all other intensities,
and to modify many other equations of the theory. In the end, taking the fundamental
equation S̄(X1 , ..., Xn ) as the point of departure would not lead to a logically consistent

57
thermodynamic fluctuation theory consistent with statistical mechanics.18 The advantage
our Eq.(4) with the entropy defined as above is that it leads to a simple and complete logical
structure fully consistent with statistical mechanics. In the limit of an infinite large system,
the fluctuations vanish and S and S̄ become numerically equal. Accordingly, Eq.(4) becomes
the fundamental equation in the sense of classical thermodynamics.
Similar considerations apply to canonical fluctuations. If the average values X̄1 , ..., X̄m
of the fluctuating extensive parameters X1 , ..., Xm are formally inserted in the fundamental
equation (4), we do not obtain the canonical entropy S̄. However, in the limit of a large
system, all canonical averages do satisfy (4) and the latter becomes the classical fundamental
equation.
For the sake of simplicity, we have only considered a particular class of thermodynamic
systems for which all arguments of the fundamental equation are additive invariants. In
the future, the theory can be generalized to include pseudo-thermodynamic and quasi-
thermodynamic variables [28]. Pseudo-thermodynamic variables are neither additive nor
conserved; an example is furnished by long-range order parameters in atomically ordered
compounds. Quasi-thermodynamic variables can be additive but need not be conserved.
They give rise to additional work terms in the total energy variation. However, the conjugate
thermodynamic forces do not become spatially homogeneous when the system reaches equi-
librium and thus cannot be considered as intensities. Examples include the total magnetic
or electric-dipole moment of the system, the conjugate variables for which are, respectively,
the magnetic and electric field. Elastic strains and stresses in inhomogeneously deformed
crystalline solids belong to the same category. All such cases require a separate treatment.
To facilitate practical applications of the theory, we have presented a regular procedure
for calculations of equilibrium fluctuations of extensive parameters, intensive parameters
and densities in systems with any number of parameters. The proposed formalism has
been demonstrated by deriving the complete set of fluctuation relations for a simple fluid
in three different ensembles. The results are summarized in Tables 1, 2 and 3 and may
present a reference value. It should be noted that these calculations treat all fluctuations in
the Gaussian approximation. Greene and Callen [26, 37] proposed a more general method
that does not rely on the Gaussian approximation and enables calculations of higher-order
correlation moments such as ∆Xi ∆Xj ∆Xk , ∆Xi ∆Xj ∆Xk ∆Xl , etc. Such moments are
rarely important and were not discussed in this paper.
The proposed fluctuation formalism has been applied to solve two problems related to
GBs in binary solid solutions. First, we developed a set of equations relating the fluctua-
tions of excess GB properties, readily accessible by atomistic simulations, to equilibrium GB
characteristics such as the excess elastic modulus, excess heat capacity and others. These
equations offer a computationally efficient approach to computing the entire set of such
properties from a single simulation run. Furthermore, the covariances between the fluctua-
tions generate a set of thermodynamic relations that can be tested by simulations to verify
the fundamental understanding of GB thermodynamics. Implementation of this program is
presently underway.
18
As already mentioned, Greene and Callen [37] resolved these difficulties by making an additional ap-
proximation that essentially identifies the canonical and micro-canonical entropies. This identification is
not satisfying for a rigorous fluctuation theory. The distinction between the two entropies was discussed in
Sec. 9.3

58
 Second, we have extended the sharp interface model of GB pre-melting to binary alloys
and derived the distribution of the fluctuating liquid-layer width w. This distribution has re-
cently been applied to calculate the disjoining potential for several GBs in Cu-Ag alloys [23].
The disjoining potential is a key GB property allowing predictions of different pre-melting
scenaria, such as the continuous melting, thin-to-thick transitions or abrupt melting after
over-heating [69]. In the version of the pre-melting model considered here, we essentially
treat the grains, the liquid layer and the solid-liquid interfaces as parts of the same infinite
reservoir. We only account for fluctuations of w due disjoining interaction. All thermody-
namic relations are in place to develop a more general model that would include fluctuations
of intensive properties of the liquid layer and/or the solid-liquid interface properties. This
generalization could be the subject of future work.
Finally, the approach demonstrated here for GBs can be readily extended to other inter-
faces. It would be interesting to study the fluctuations of excess properties of solid-liquid
and especially solid-solid interfaces to determine their excess elastic and thermal properties
and test some of the basic equations of interface thermodynamics.

Acknowledgements. This work was supported by the National Science Foundation,

Division of Materials Research, the Metals and Metallic Nanostructures Program.

References
[1] D. Frenkel and B. Smit. Understanding molecular simulation: from algorithms to
applications. Academic, San Diego, second edition, (2002).

[2] D. P. Landau and K. Binder. A guide to Monte Carlo simulations in statistical

physics. Cambridge University Press, Cambridge, third edition, (2009).

[3] J. J. Hoyt, Z. T. Trautt and M. Upmanyu. Fluctuations in molecular dynamics

simulations. Mathematics and Computers in Simulation 80 (2010) 1382–1392.

[4] M. Parrinello and A. Rahman. Strain fluctuations and elastic constants. J. Chem.
Phys. 76 (1982) 2662–2666.

[5] A. A. Gusev, M. M. Zehnder and U. W. Suter. Fluctuation formula for elastic

constants. Phys. Rev. B 54 (1996) 1–4.

[6] Z. Zhou and B. Joos. Fluctuation formula for elastic constants of an arbitrary
system. Phys. Rev. B 66 (2002) 054101.

[7] Z. Cui, Y. Sun and J. Qu. Combination method for the calculation of elastic
constants. Phys. Rev. B 75 (2007) 214101.

[8] M. T. Meyers, J. M. Rickman and T. J. Delph. The calculation of elastic

constants from displacement fluctuations. J. Appl. Phys. 98 (2005) 066106.

[9] Y. Mishin. Calculation of open and closed system elastic coefficients for multicompo-
nent solids. Phys. Rev. B 91 (2015) 224107.

59
[10] P. G. Debenedetti. Fluctuation-based computer calculation of partial molar proper-
ties. I. Molecular dynamics simulation of constant volume fluctuations. J. Chem. Phys.
86 (1987) 7126–7137.

[11] P. G. Debenedetti. Fluctuation-based computer calculation of partial molar proper-

ties. II. A numerically accurate method for the determination of partial molar energies
and enthalpies. J. Chem. Phys. 88 (1988) 2681–2684.

[12] R. Davidchack, J. Morris and B. Laird. The anisotropic hard-sphere crystal-melt

interfacial free energy from fluctuations. J. Chem. Phys. 125 (2006) 094710.

[13] M. Amini and B. B. Laird. Crystal-melt interfacial free energy of binary hard
spheres from capillary fluctuations. Phys. Rev. B 78 (2008) 144112.

[14] J. J. Hoyt, M. Asta and A. Karma. Method for computing the anisotropy of the
solid-liquid interfacial free energy. Phys. Rev. Lett. 86 (2001) 5530–5533.

[15] S. M. Foiles and J. J. Hoyt. Computation of grain boundary stiffness and mobility
from boundary fluctuations. Acta Mater. 54 (2006) 3351–3357.

[16] J. J. Hoyt, D. Olmsted, S. Jindal, M. Asta and A. Karma. Method for

computing short-range forces between solid-liquid interfaces driving grain boundary
premelting. Phys. Rev. E 79 (2009) 020601.

[17] H. Song, S. J. Fensin, M. Asta and J. J. Hoyt. A molecular dynamics simulation

of (110) surface premelting in Ni. Scripta Mater. 63 (2010) 128–131.

[18] S. J. Fensin, D. Olmsted, D. Buta, M. Asta, A. Karma and J. J. Hoyt.

Structural disjoining potential for grain-boundary premelting and grain coalescence
from molecular-dynamics simulations. Phys. Rev. E 81 (2010) 031601.

[19] A. Karma, Z. T. Trautt and Y. Mishin. Relationship between equilibrium fluc-

tuations and shear-coupled motion of grain boundaries. Phys. Rev. Lett. 109 (2012)
095501.

[20] A. Adland, A. Karma, R. Spatschek, D. Buta and M. Asta. Phase-field-

crystal study of grain boundary premelting and shearing in bcc iron. Phys. Rev. B 87
(2013) 024110.

[21] Y. Mishin. Calculation of the γ/γ′ interface free energy in the Ni-Al system by the
capillary fluctuation method. Modeling Simul. Mater. Sci. Eng. 22 (2014) 045001.

[22] F. Schmitz, P. Virnau and K. Binder. Determination of the origin and magnitude
of logarithmic finite-scale effects on interfacial tension: Role of interfacial fluctuations
and domain breathing. Phys. Rev. Lett. 112 (2014) 125701.

[23] J. Hickman and Y. Mishin. Atomistic modeling of grain boundary pre-melting in

alloy systems. To be published (2015).

60
[24] J. W. Gibbs. On the equilibrium of heterogeneous substances. In: The collected works
of J. W. Gibbs, volume 1. Yale University Press, New Haven, 1948 55–349.

[25] A. Münster. Classical Thermodynamics. Wiley-Interscience, London - New York -

Sydney - Toronto, (1970).

[26] H. B. Callen. Thermodynamics and an introduciton to the physical theories of equi-

librium thermostatics and irreversible thermodynamics. Wiley, New York - London -
Sydney, (1960).

[27] H. B. Callen. Thermodynamics and an introduciton to thermostatistics. Wiley, New

York, second edition, (1985).

[28] L. Tisza. The thermodynamics of phase equilibrium. Annals of Physics 13 (1961)

1–92.

[29] C. Carathéodory. Untersuchungen Über die Grundlagen der Thermodynamik.

Mathematische Annalen 67 (1909) 355–386.

[30] P. Ehrenfest and T. Ehrenfest. The conceptual foundations of the statistical

approach in mechanics. Dover books on Physics. Dover, Ithaca, New York, (1959).

[31] P. T. Landsberg. Main ideas in the axiomatics of thermodynamics. Pure Appl.

Chem. 22 (1970) 215–227.

[32] J. Kestin. A simple, unified approach to the first and second laws of thermodynamics.
Pure Appl. Chem. 22 (1970) 511–518.

[33] C. Gurney. Principles of classical thermodynamics. Pure Appl. Chem. 22 (1970)

519–526.

[34] K. A. Masavetas. The problem of alternative formalisms for the mathematical foun-
dations of chemical thermodynamics. Mathl. Comput. Modelling 10 (1988) 629–635.

[35] R. J. J. Jongschaap and H. C. Ottinger. Equilibrium thermodynamics —

Callen’s postulational approach. J. Non-Newtonian Fluid Mech. 96 (2001) 5–17.

[36] F. Weinhold. Classical and geometrical theory of chemical and phase thermodynam-
ics. Wiley, Hoboken, New Jersey, (2009).

[37] R. F. Greene and H. B. Callen. On the formalism of thermodynamic fluctuation

theory. Phys. Rev. 83 (1951) 1231–1235.

[38] H. Callen. Thermodynamic fluctuations. American Journal of Physics 33 (1965)

919–922.

[39] L. Tisza and P. M. Quay. The statistical mechanics of phase equilibrium. Annals
of Physics 25 (1963) 49–90.

61
[40] A. Einstein. Theorie der Opaleszenz von homogenen Flüssigkeiten und Flüssigkeit-
enmischungen in der Nähe der kritischen Zustandes. Annalen der Physik 33 (1910)
1275–1298.

[41] J. W. Gibbs. Elementary principles of statistical mechanics. In: The collected works
of J. W. Gibbs, volume 2. Yale University Press, New Haven, 1948 1–207.

[42] Y. D. Rudoi and A. D. Sukhanov. Thermodynamic fluctuations within the Gibbs

and Einstein approaches. Phys. Usp. 43 (2000) 1169–1199.

[43] L. Szilard. Über die Ausdehnung der phänomenologischen Thermodynamik auf die
Schwankungserscheinungen. Zeitschrift für Physik 32 (1925) 753–788.

[44] A. I. Khinchin. Mathematical foundations of statistical mechanics. Dover, New York,

(1949).

[45] D. J. Evans, E. G. D. Cohen and G. P. Morris. Probability of second law

violations in shearing steady states. Phys. Rev. Lett. 71 (1993) 2401.

[46] L. Bertini, A. de Sole, D. Gabrielli, G. Jona-Lasinio and C. Landim. Macro-

scopic fluctuation theory. Rev. Mod. Phys. 87 (2015) 593–363.

[47] M. S. S. Challa and J. H. Hetherington. Gaussian ensemble as an interpolating

ensemble. Phys. Rev. Lett. 60 (1988) 77–80.

[48] A. Hüller. Finite size scaling at first order phase transitions? Z. Phys. B 95 (1994)
63–66.

[49] M. Promberger and A. Hüller. Microcanonical analysis of a finite three-

dimensional Ising system. Z. Phys. B 97 (1995) 341–344.

[50] R. S. Johal, A. Planes and E. Vives. Statistical mechanics in the extended

Gaussian ensemble. Phys. Rev. E 68 (2003) 056113.

[51] L. Velazquez and S. Curilef. On the thermodynamic stability of macrostates with

negative heat capacities. J. Stat. Mech. 2009 (2009) P03027.

[52] L. Velazquez and S. Curilef. Extending canonical Monte Carlo methods. J. Stat.
Mech. 2010 (2010) P02002.

[53] L. Velazquez and S. Curilef. Extending canonical Monte Carlo methods: II. J.
Stat. Mech. 2010 (2010) P04026.

[54] B. Sadigh, P. Erhart, A. Stukowski, A. Caro, E. Martinez and

L. Zepeda-Ruiz. Scalable parallel Monte Carlo algorithm for atomistic simulations of
precipitation in alloys. Phys. Rev. B 85 (2012) 184203.

[55] B. Sadigh and P. Erhart. Calculation of excess free energies of precipitates via
direct thermodynamic integration across phase boundaries. Phys. Rev. B 86 (2012)
134204.

62
[56] H. Nyquist. Thermal agitation of electric charge in conductors. Phys. Rev. 32 (1928).

[57] H. B. Callen and T. A. Welton. Irreversibility and generalized noise. Phys. Rev.
83 (1951) 34–40.

[58] L. D. Landau and E. M. Lifshitz. Statistical Physics, Part I , volume 5 of Course

of Theoretical Physics. Butterworth-Heinemann, Oxford, third edition, (2000).

[59] R. W. Gerling and A. Hüller. First order phase transitions studied in the dynam-
ical ensemble: The q-state Potts model as a test case. Z. Phys. B 90 (1993) 207–214.

[60] A. P. Sutton and R. W. Balluffi. Interfaces in Crystalline Materials. Clarendon

Press, Oxford, (1995).

[61] P. W. Voorhees and W. C. Johnson. The thermodynamics of elastically stressed

crystals. In: H. Ehrenreich and F. Spaepen, eds., Solid State Physics, volume 59.
Elsevier Academic Press, 2004 2–201.

[62] T. Frolov and Y. Mishin. Thermodynamics of coherent interfaces under mechanical

stresses. I. Theory. Phys. Rev. B 85 (2012) 224106.

[63] F. Larché and J. W. Cahn. A linear theory of thermochemical equilibrium of solids

under stress. Acta Metall. 21 (1973) 1051–1063.

[64] F. C. Larché and J. W. Cahn. A nonlinear theory of thermochemical equilibrium

under stress. Acta Metall. 26 (1978) 53–60.

[65] F. C. Larché and J. W. Cahn. The interactions of composition and stress in

crystalline solids. Acta Metall. 33 (1985) 331–367.

[66] R. Courant and F. John. Introduction to calculus and analysis, volume II. Springer
Verlag, New York, (1989).

[67] T. Frolov and Y. Mishin. Thermodynamics of coherent interfaces under mechanical

stresses. II. Application to atomistic simulation of grain boundaries. Phys. Rev. B 85
(2012) 224107.

[68] R. Spatschek, A. Adland and A. Karma. Structural short-range forces between

solid-melt interfaces. Phys. Rev. B 87 (2013) 024109.

[69] Y. Mishin, W. J. Boettinger, J. A. Warren and G. B. McFadden. Thermo-

dynamics of grain boundary premelting in alloys. I. Phase field modeling. Acta Mater.
57 (2009) 3771–3785.

63