Thermal Equilibrium

APAM 1601
Columbia University

Introduction

In this notebook, we introduce the notion of a "thermal equilibrium" and investigate simple properties like mean energy and average particle density using a biased random walk model.

A simple atmosphere

Our starting point is a model for an ideal atmosphere. The gas within the atmosphere is confined, at the top, by a downward pointing gravity and, at the bottom, by the surface of the planet. The gas will be characterized by a "temperature", and the temperature will bias how likely the particle will take a random step upward (or downward). If the motion of a gas particle causes an increase in the overall energy of the atmosphere, then this motion must be restricted by a "transition probability" containing the Boltzman factor, Exp[–ΔU/T], where ΔU is the change of energy of the atmosphere.

Monte Carlo Simulation (Metropolis, et al., 1953)

We will use the famous "Metropolis" algorithm (Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller, "Equation of State Calculations by Fast Computing Machines," J. Chem. Phys. 21 (1953) pp. 1087-1092) to simulate the thermal fluctuations of our model atmosphere. A Monte Carlo simulation uses a random number generator to simulate the random fluctuations which occur in a large collection of particles.

The "Metropolis" algorithm is one of the is one of the most important and widely-used numerical methods in the field of statistical physics.(The story goes that the authors searched in vain for another physicist named Metropolis,so that the paper could have the symmetrical author list Metropolis,Metropolis,Rosenbluth,Rosenbluth,Teller and Teller! Metropolis (1915-1999) was an early pioneer in Monte Carlo methods and group leader at Los Alamos.His name appears first on the paper (cited above),and as a result most people refer to this method as the "Metropolis" Monte Carlo method.However,Marshall Rosenbluth (1927-2003) actually invented the algorithm.Marshall is perhaps one of the greatest pioneers in plasma physics and fusion energy science,and most plasma physicists refer to this algorithm as the Rosenbluth-Teller-Metropolis method.)

Marshall Rosenbluth (1927-2003)

For our atmosphere model, we wish to simulate the random thermal agitation of a large collection particles interacting only with gravity (and the surface of the Earth.) Basically, each particle makes random steps upward or downward. However, these probabilities will no longer be "unbiased" (i.e. 50-50 up or down.) Instead, for each potential step, we check the energy change involved and reject those steps that are thermally unlikely. Our simulation will start at an initial state and continue through a series of time steps towards some sort of "thermal equilibrium" or a "relaxed state". 
How is this done?

We will follow a set of specific rules.
Step 1) Initialize the positions of the particles in the atmosphere.
Step 2) Select a random particle. Make a potential random step (up or down).
Step 3) Calculate the energy change to the entire system if the particle makes the step.
Step 4a) If the energy change is negative, accept the step.
Step 4b) If the energy change is positive, pick a random number between 0 and 1. If the random number is less than Exp[ - (energy change)/kT], then accept the step.
Step 4c) Otherwise, reject the step and leave the position unchanged.
Step 5) Go to Step 2 and repeat (for a very long time....)

The Potential Energy of an "Ideal" Atmosphere

What is the potential energy of our ideal atmosphere? For simplicity, we will neglect all interactions between particles. Essentially, the particles feel only a gravitational force, and their dynamics are independent of other particles.

The total potential energy for N particles will be
	(Subscript[U, tot]) = (∑_(i=1)^N)ϕ((Subscript[x, i]))
where ϕ((Subscript[x, i])) = m g (Subscript[x, i]), is the gravitational potential. Particles are at high potential energy when they are high in the sky and lowering potential energy when they are close to the ground. (The ground will be at x = 0.) In order to simulate the ground, we'll take ϕ = (10^5) (a huge number) whenever x < 0.

Computation of the change in the atmosphere's energy if a particle takes a random step is easy to specify. The change in energy is ΔU= m g Δx. If ΔU > 0, the upward step is accepted if a random number between 0 and 1 is less than Exp[-ΔU / T].

We characterize thermal equilibrium with a temperature. The energy (per degree of freedom, one in our case) of an average particle is (Subscript[k, B])T, where (Subscript[k, B])= 1.3807 × (10^(-23)) J/K is Boltzman's constant.We will take (Subscript[k, B])=1, so temperature is in energy units. With m g= 1, the potential energy at x = 1 is equal to a thermal energy at T = 1.

A Gravitationally-Biased Random Walker

We now have everything needed to define a "gravitationally-biased" random walker used to simulate particle confinement within a collisionless atmosphere. For simplicity, we set m = g = 1.

s[]:=Δx RandomReal[{-1,1}]

step[x_]:=Module[{trialStep,ΔU},trialStep=s[];If[x+trialStep<0,reject+=1;Return[x]];ΔU=trialStep;If[ΔU<0,x+trialStep,If[RandomReal[]<(E^(-((ΔU)/T))),x+trialStep,reject+=1;x]]]

T = 0.1; 
Δx = 0.1;
reject = 0;
NestList[step, 1.0, 10]

reject

reject = 0;
(w1 = NestList[step,1.0,50000];Print["Rejection Rate = ", reject/50000 100.0 PerCent];) // Timing

We generally want to have the reject rate be between 30 and 40%. As Δx increases, then the rejection rate increases. To decrease the rejection rate, the value for Δx must be decreased. Δx can be larger for larger T.

ListPlot[w1[[Range[1,50000,100]]],AxesLabel->{"n/100","x"}]

walkAtmosphere[n_] := NestList[step, 1.0, n] (* initial height is x = 1.0 *)

(Exercise for thought: does the initial condition affect the final thermal equilibrium?)

One thousand walkers…

We take 250,000 steps…

nWalkers = 1000;
reject = 0;
walkList = Table[walkAtmosphere[250],{nWalkers}];
Print["Rejection Rate = ", reject/250000 100.0 PerCent];

walkList // Dimensions (* 1000 walks 500 steps long *)

ListAnimate[Table[ListPlot[walkList[[All,n]], PlotRange -> {0, 1.0}, AxesLabel -> {"walker", "x, n= " <> ToString[n-1]}], {n, 1, 251, 10}],AnimationRunning -> False]

Statistics

BinCounts

BinCounts[walkList[[All,1]], {0.0, 1.5,0.05}]

BinCounts[walkList[[All,100]], {0.0, 1.5,0.05}]

grid = Table[BinCounts[walkList[[All,n]], {0.0, 1.5,0.05}], {n, 1, 251}];

ListDensityPlot[grid,Mesh->False,ColorFunction->Hue,PlotRange->{1,100},FrameLabel->{"Box(x/Δx), "<>ToString[nWalkers]<>" Walkers","n"}]

p1=ListPlot[grid[[15,All]],PlotRange->All,PlotStyle->{Red,PointSize[0.01]}];
p2=ListPlot[grid[[200,All]],PlotRange->All,PlotStyle->{Blue,PointSize[0.01]}];
Show[p1,p2,AxesOrigin->{0,0},AxesLabel->{"Box(x)","Count"}]

Histogram

p3=Histogram[walkList[[All,100]],{0,1.5,1.5/45},PlotLabel->"h[j,n], n= "<>ToString[100],AxesLabel->{"x",""}]

histograms = Table[p1 = Histogram[walkList[[All,n]],{0,1.5,1.5/45},PlotLabel -> "h[j,n], n= " <> ToString[n], AxesLabel -> {"x/Δx",""}];
Show[p1, PlotRange -> {{0.0,1.5},{0, 100}},AxesOrigin -> {0,0}], {n, 1, 251, 10}];

ListAnimate[histograms,AnimationRunning -> False]

Mean and Standard Deviation

avgX=Mean/@Transpose[walkList];
avgPlot1=ListPlot[avgX,AxesLabel->{"n","<x>"}]

varX=Variance/@Transpose[walkList];
ListPlot[varX,AxesLabel->{"n","Var(x)"}]

stdDevX=StandardDeviation/@Transpose[walkList];
ListPlot[stdDevX,AxesLabel->{"n","StdDev(x)"}]

The standard deviation indicates the magnitude of the "fluctuations" in x.

Entropy

Entropy is a measure of the disorder of the distribution of random walkers. For particles with negligible kinetic energy, entropy is defined simply as S(t) = – ∫P(x,t)  Log[P(x,t)]dx, where P(x,t) is the probability of a walker being at cell x at time t. Initially, all of the walkers are located at x = 1 with probability of 100%, The entropy is zero. As the probability become more uniform in x, then the entropy will increase. As the particles accelerate to  the ground, the entropy appears to increase and then decrease. In a non-equilibrium transient state of our simulation, our definition is not very useful. Recall, our estimate of the probability within a gird-box is the ratio, P(x,t) ≈ (# of walkers within grid box x)/(total number of walkers).

Dimensions[grid]

entropyGrid = Table[Plus @@ (If[# > 1,-N[(#/nWalkers) Log[(#/nWalkers)]],0]& /@  grid[[t,All]]), {t, 1, 250}];

ListPlot[entropyGrid,AxesLabel->{"n","Entropy[n], "<>ToString[nWalkers]<>" Walkers"},Joined->True,PlotRange->{0,4}]

Potential Energy

The average potential energy of the particles within our model atmosphere shows relaxation towards lower energy as thermal equilibrium is approached. For our simple model, ⟨U⟩ = (∑_(i=1)^nWalkers)ϕ((Subscript[x, i])) = m g ⟨x⟩.

avgU=avgX;
ListPlot[avgU,AxesLabel->{"n","Avg Potential Energy"}]

In thermal equilibrium, the average potential energy should equal the temperature…

equilU = Mean[avgX[[Range[200,250]]]]

Using the Boltzman distribution for the probability function, the equilibrium value of x is proportional to ∫x Exp[-U(x)/T]dx.

Integrate[x Exp[-x/T], {x, 0, Infinity}]/Integrate[Exp[-x/T], {x, 0, Infinity}]

The mean fluctuation in the potential energy in time is (Sqrt[⟨((((Subscript[U, n]) - (Subscript[U, avg])))^2)⟩]), or

avgδU = Sqrt[Mean[Table[(avgU[[n]]  - equilU)^2, {n, 200, 250}]]]

And this gives a fluctuation about 1/40 the size of ⟨U⟩ (about ± 2.5 %).

Convergence in Monte Carlo Calculations

Perhaps the more important of the practical questions concerning Monte Carlo calculations is one of convergence. How many particles are required to estimate accurate averages? The answer to this question can be found in what we already know about random walks. The mean square distance from the average increases linearly with the number of steps in a random process. If we are trying to estimate the true equilibrium potential energy ((Subscript[U, true])) of a system by averaging over N steps, then the calculated average would be ⟨U⟩ ~ ( N (Subscript[U, true]) + (Sqrt[N ⟨((δU)^2)⟩])) /N. The term (Sqrt[N ⟨((δU)^2)⟩]) ~ ΔU is the deviation of the potential energy due to random diffusion. The value of the average is ⟨U⟩ ~ (Subscript[U, true]) + (Sqrt[⟨((δU)^2)⟩/N]). As N -> ∞, the difference (i.e. the error) between the Monte Carlo estimate of the error and the true equilibrium potential energy decreases as 1/(Sqrt[N]).

For our first calculation, we had 1000 walkers.If we had only 62 walkers (16 times fewer), the error in our estimate of ⟨U⟩ will be four (4) times larger.

T = 0.1; 
Δx = 0.1;

nWalkers = 62;
reject = 0;
walkList = Table[walkAtmosphere[250],{nWalkers}];
Print["Rejection Rate = ", reject/(250 62) 100.0 PerCent];

avgX=Mean/@Transpose[walkList];
avgPlot2=ListPlot[avgX,AxesLabel->{"n","<x>"},PlotRange->All,AxesOrigin->{0,0},PlotStyle->Red]

equilU = Mean[avgX[[Range[200,250]]]]

avgU = avgX;

avgδU = Sqrt[Mean[Table[(avgU[[n]]  - equilU)^2, {n, 200, 250}]]]

This uncertainty (± 8.3 %) is about 4 times larger than with 1000 ideal atmospheric particles, as expected. The reality of Monte Carlo calculations is the need for a large number of realizations with which to find averages!

Show[avgPlot1,avgPlot2,PlotLabel->"62 vs 1000 Walkers"]

Summary

The Rosenbluth-Teller-Metropolis algorithm is a general prescription for finding properties of systems at thermal equilibrium. Particles, or components of the larger system, are allowed to make random (thermal motions) and these motions are tested using a Boltzman probability distribution function. If the motion causes the system energy to increase, then the motion is rejected if the change in energy is too large. In this way, the "simulated" random motions of the system eventually relax to a state of thermal equilibrium. The properties of the thermal equilibrium can be estimated by taking averages of the simulation. The accuracy of the averages (as compared to the "true" equilibrium properties) improves as 1/(Sqrt[N])when the number of particles (and configuration steps) increase.