Foundations of Computational Materials Science/Introduction and Basic Concepts of Classical Statistical Mechanics

Variables in many-particle systems

The state of a N particles system is completely characterized by a set of 2×3 coordinates, the positions and momenta of the particles, listed as

${\begin{aligned}\left\{q_{x}^{\left(1\right)},q_{y}^{\left(1\right)},q_{z}^{\left(1\right)},q_{x}^{\left(2\right)},q_{y}^{\left(2\right)},q_{z}^{\left(2\right)},...,q_{x}^{\left(N\right)},q_{y}^{\left(N\right)},q_{z}^{\left(N\right)}\right\}:=&{\vec {q}}\\\left\{p_{x}^{\left(1\right)},p_{y}^{\left(1\right)},p_{z}^{\left(1\right)},p_{x}^{\left(2\right)},p_{y}^{\left(2\right)},p_{z}^{\left(2\right)},...,p_{x}^{\left(N\right)},p_{y}^{\left(N\right)},p_{z}^{\left(N\right)}\right\}:=&{\vec {p}},\\\end{aligned}}$

where $q_{i}^{\left(j\right)}$ and $p_{i}^{\left(j\right)}$ are the the position and the impulse of the ith particle in the j coordinate respectively. A mapping in the indices is widely used to make the notations easier as follows: $\left({{q}_{1}},{{q}_{2}},{{q}_{3}},{{q}_{4}},{{q}_{5}},{{q}_{6}},...,{{q}_{3N-2}},{{q}_{3N-1}},{{q}_{3N}}\right):=\left(q_{x}^{\left(1\right)},q_{y}^{\left(1\right)},q_{z}^{\left(1\right)},q_{x}^{\left(2\right)},q_{y}^{\left(2\right)},q_{z}^{\left(2\right)},...,q_{x}^{\left(N\right)},q_{y}^{\left(N\right)},q_{z}^{\left(N\right)}\right),$

and for the impulse as well.

Both $q_{i}^{\left(j\right)}$ and $p_{i}^{\left(j\right)}$ are a function of the time t, and the equation of their time derivative (which is called the equation of motion) is given by the Hamilton's equations:

${{\dot {q}}_{k}}={\frac {\partial H\left({\vec {q}},{\vec {p}}\right)}{\partial {{p}_{k}}}}\quad {{\dot {p}}_{k}}=-{\frac {\partial H\left({\vec {q}},{\vec {p}}\right)}{\partial {{q}_{k}}}}$

where $H\left({\vec {q}},{\vec {p}}\right)$ is the Hamiltonian function of the system and called energy and is given by

$H\left({\vec {q}},{\vec {p}}\right)=\sum \limits _{k=1}^{3N}{{\frac {p_{k}^{2}}{2{{m}_{k}}}}+V\left({\vec {q}}\right)},$

where $m_{k}$ is the mass of the kth particle and $V\left({\vec {q}}\right)$ is the potential (external and internal).

Trajectory

We can see that the N particle system and it's evolution is described by $6N$ coordinates. It can be envisaged as N points moving in a 6-dimensional space (3 space coordinates and 3 impulse coordinates), or as a single point moving in a $6N$-dimensional space, called phase space. In this case, we can take a look at the trajectory of the system, i.e. the image of the mapping from the time to the $6N$ dimensional space. This is a $1D$ curve embedded in a $6N$-dimensional space, and it has the following properties:

1. starting from two different points of the phase field (investigated the same 6N particle system with different initial conditions), the two curves can never intersect each other. (Although, they can be as close to each other, as they cannot be resolved at a numeric precision.)
2. a curve cannot intersect itself. (Although, at a given numeric precision it can be impossible to determine which direction the one curve continues in from the 'quasi' intersection.)
3. the mapping could take the same value in the $6N$-dimensional space, but in this case, from this point on, the curve must form a closed loop and the motion is periodic.

The aim of statistical mechanics

• When dealing with macroscopic systems, we are in general not interested in the complete microscopic information
• There are in general many microsystems which correspond to the same macroscopic state of the system
• We try to obtain, starting from the dynamics of individual microsystems, a statistical description of the properties of macrosystems. On the way, we try to lose as much of the detailed (but mostly irrelevant) information on the microsystems dynamics as possible
• Deriving macro-equations from micro-dynamics

Statistical ensemble

To obtain a statistical description of a macrosystem, we consider not only one microsystem but many of them corresponding to a set of macro-variables (a statistical ensemble). This ensemble is characterized by a probability to find the system in any given microstate that is consistent with the macro-variables. This probability is given by the N-particle density function:

${{\rho }_{N}}\left({{q}_{1}},{{q}_{2}},...,{{q}_{3N}},{{p}_{1}},{{p}_{2}},...,{{p}_{3N}}\right)d{{q}_{1}}d{{q}_{2}}...d{{q}_{3N}}d{{p}_{1}}d{{p}_{2}}...d{{p}_{3N}}$

We can envisage this function as a density of points in the $6N$ dimensional phase space. Each point corresponds to a microsystem. We require normalization (the system must be in some microstate):

$\int _{\Omega }{{{\rho }_{N}}\left({{q}_{1}},{{q}_{2}},...,{{q}_{3N}},{{p}_{1}},{{p}_{2}},...,{{p}_{3N}}\right)d{{q}_{1}}d{{q}_{2}}...d{{q}_{3N}}d{{p}_{1}}d{{p}_{2}}...d{{p}_{3N}}}=1,$ where $\Omega$ is the 6N dimensional phase space where the particle could ever be at.

Since the evolution is deterministic, the ensemble is fully defined by its initial density function. This can be thought of as a statistical rule for constructing initital conditions in multiple simulations.

Example

Let us consider a set of simulations where particles are initially placed with equal probability inside a cube of edge length L and with initial momenta in x, y, and z equi-distributed between 0 and P. This corresponds to the initial density function

${{\rho }_{N}}\left({{q}_{1}},{{q}_{2}},...,{{q}_{3N}},{{p}_{1}},{{p}_{2}},...,{{p}_{3N}}\right)={{\left({\frac {1}{L}}\right)}^{3N}}{{\left({\frac {1}{P}}\right)}^{3N}}$

Equation of Motion (EOM) for the phase space density

Consider the motion of all phase space points originally contained in a volume ${{G}_{0}}$. At time t this evolves into a volume ${{G}_{t}}$. Because the trajectories of two microsystems cannot intersect, the number of phase space points contained in ${{G}_{t}}$ is conserved, therefore it's total time derivative is 0:

${\frac {d}{dt}}\int _{{G}_{t}}{{{\rho }_{N}}\left({\vec {q}},{\vec {p}}\right){{d}^{3N}}q{{d}^{3N}}p}=0$

Let's look at the consequence! Move the time derivation into the integral:

${\begin{aligned}&{\frac {d}{dt}}\int _{{G}_{t}}{{{\rho }_{N}}\left({\vec {q}},{\vec {p}}\right){{d}^{3N}}q{{d}^{3N}}p}=0=\\&=\int _{{G}_{t}}{\left[{\frac {\partial {{\rho }_{N}}}{\partial t}}+\sum \limits _{k=1}^{3N}{\left({\frac {\partial {{\rho }_{N}}}{\partial {{q}_{k}}}}{{\dot {q}}_{k}}+{\frac {\partial {{\rho }_{N}}}{\partial {{p}_{k}}}}{{\dot {p}}_{k}}\right)+{{\rho }_{N}}\sum \limits _{k=1}^{3N}{\left({\frac {\partial {{\dot {q}}_{k}}}{\partial {{q}_{k}}}}+{\frac {\partial {{\dot {p}}_{k}}}{\partial {{p}_{k}}}}\right)}}\right]{{d}^{3N}}q{{d}^{3N}}p}\\\end{aligned}}$

Since it is true for arbitrary volume $G_{t}$, it can be only true if the integrand itself is 0:

${\frac {\partial {{\rho }_{N}}}{\partial t}}+\sum \limits _{k=1}^{3N}{\left({\frac {\partial {{\rho }_{N}}}{\partial {{q}_{k}}}}{{\dot {q}}_{k}}+{\frac {\partial {{\rho }_{N}}}{\partial {{p}_{k}}}}{{\dot {p}}_{k}}\right)}+{{\rho }_{N}}\sum \limits _{k=1}^{3N}{\left({\frac {\partial {{\dot {q}}_{k}}}{\partial {{q}_{k}}}}+{\frac {\partial {{\dot {p}}_{k}}}{\partial {{p}_{k}}}}\right)}=0$

Substituting into ${{\dot {q}}_{k}}$ and ${{\dot {p}}_{k}}$ from the Hamilton's equation

${{\dot {q}}_{k}}={\frac {\partial H}{\partial {{p}_{k}}}}\quad {{\dot {p}}_{k}}=-{\frac {\partial H}{\partial {{q}_{k}}}},$

we get

${\frac {\partial {{\rho }_{N}}}{\partial t}}+\sum \limits _{k=1}^{3N}{\left({\frac {\partial {{\rho }_{N}}}{\partial {{q}_{k}}}}{\frac {\partial H}{\partial {{p}_{k}}}}-{\frac {\partial {{\rho }_{N}}}{\partial {{p}_{k}}}}{\frac {\partial H}{\partial {{q}_{k}}}}\right)+}{{\rho }_{N}}\sum \limits _{k=1}^{3N}{\left({\frac {{{\partial }^{2}}H}{\partial {{q}_{k}}\partial {{p}_{k}}}}-{\frac {{{\partial }^{2}}H}{\partial {{p}_{k}}\partial {{q}_{k}}}}\right)}=0$

Using the Young's theorem on second derivatives, the right term disappear:

${{\rho }_{N}}\sum \limits _{k=1}^{3N}{\left({\frac {{{\partial }^{2}}H}{\partial {{q}_{k}}\partial {{p}_{k}}}}-{\frac {{{\partial }^{2}}H}{\partial {{q}_{k}}\partial {{p}_{k}}}}\right)}=0$

What is left, the EOM of the phase space density, and it has the structure of the continuity equation for an incompressible fluid:

${\frac {\partial {{\rho }_{N}}}{\partial t}}+\sum \limits _{k=1}^{3N}{\left({\frac {\partial {{\rho }_{N}}}{\partial {{q}_{k}}}}{\frac {\partial H}{\partial {{p}_{k}}}}-{\frac {\partial {{\rho }_{N}}}{\partial {{p}_{k}}}}{\frac {\partial H}{\partial {{q}_{k}}}}\right)}=0$

It follows that the phase space density in the vicinity of a given system is conserved in time (Liouville's theorem). A widely used form of the EOM is

${\frac {\partial {{\rho }_{N}}}{\partial t}}=-\sum \limits _{k=1}^{3N}{\left({\frac {\partial {{\rho }_{N}}}{\partial {{q}_{k}}}}{\frac {\partial H}{\partial {{p}_{k}}}}-{\frac {\partial {{\rho }_{N}}}{\partial {{p}_{k}}}}{\frac {\partial H}{\partial {{q}_{k}}}}\right)}=:\left\{H\left({\vec {q}},{\vec {p}}\right),{{\rho }_{N}}\right\},$

where on the right side the Poisson bracket is introduced. The notation of Poisson bracket has the advantage that the analogy with quantum mechanics can be more easily seen later on.

The recurrence theorem

Consider a phase space of finite volume G and let a system be at time $t=0$ be contained within a sub-volume g (note: g may be arbitrarily small). Then the system will, in the course of time, return with probability 1 into g.

Proof

1. We consider a phase space filled with a uniform density (microcanonical ensemble) and calculate the volume flow of points out of g. The rate r of this flow is constant in time. Hence the cumulative volume that has flown out of g at time t is $rt$.
2. Let f denote the fraction of this volume that has not returned to g until time t.
3. The phase space density is incompressible, hence the inequality $frt<G-g<G$ must be fullfilled. As $t\to \infty$ it follows $f\to 0$.