Rubber band dynamics III: Free energy

(Sethna, “Entropy, Order Parameters, and Complexity”, ex. XXX)
© 2024, James P. Sethna, all rights reserved.

Statistical mechanics is a complete theory for the static properties of Hamiltonian systems: the probability of a snapshot of the system having any particular configuration. It constrains the dynamics of the system (entropy cannot decrease, ...) but different microscopic physics or simulation methods can change how a system evolves in time.

Here we study the statics and two kinds of dynamics in the entropic rubber band model, introduced in Exercise 5.2 in the microcanonical ensemble, and analyzed in Exercise 6.16 in the fixed-force ensemble. In ‘RB Dynamics I’ we added a parabolic potential energy to the model, and found a transition between a state with one equilibrium length at zero and a state with two equilibrium lengths.

We start by analyzing the static properties of the rubber band model in an ensemble fixing the external force on the random chain, and with the external parabolic potential.

We reformulate our model in terms of segment orientations  𝐬. Each of the  𝑁 segments of the rubber band has length one and can point in one of two directions  𝑠_𝑖 = ±1, with the rubber band length  𝐿=∑ 𝑠_𝑖.  𝐹 is the external force on the tip of the rubber band, and the external potential is  (1/2) 𝛼 𝐿^2, so
						H[L]=-F L-(((L^2) α)/2)
			H[s]=-F (∑_(i=1)^N)(Subscript[s, i])-(1/2) α ((((∑_(i=1)^N)(Subscript[s, i])))^2)
			          =-F (∑_(i=1)^N)(Subscript[s, i])-(1/2) α ((∑_(i=1)^N)(Subscript[s, i]))((∑_(j=1)^N)(Subscript[s, j]))
				    =-F (∑_(i=1)^N)(Subscript[s, i])-(1/2) α (∑_(j=1)^N)((∑_(i=1)^N)(Subscript[s, i]) (Subscript[s, j]))- (1/2) α N
where  𝐬={𝑠_1,…,𝑠_(2^N)} runs over all  (2^N) possible segment orientations. Here the last formula for  H connects our rubber band problem to the well-studied infinite-range Ising model with  𝐽/N = 𝛼 and  𝐹 = 𝐻.

Our ensemble fixes the force 𝐹 and the coupling  𝛼. The partition function sums the Boltzmann weight over all possible segment orientation patterns,
			Z[F,α]==∑ (e^(-H[𝐬]))==∑(Subscript[Z, F,α])[L]
 where we measure energies in units of kB T and entropy in nats (kB = 1). Thus we shall change F and α up and down rather than changing the temperature T down and up to explore the behavior.

What is this last decomposition into  𝑍(𝐿)? Since  H(𝐬) depends on the spins only through their sum, we can count the number of segment configurations  Ω(𝐿)=𝑁!/((N+L)/2)!((N-L)/2)! and weigh them by  exp(−H(𝐿)):
		𝑍(𝐿) = Ω(𝐿)exp(−H(𝐿)) = exp(−(H(𝐿)−𝑆(𝐿))).
 where  𝑆(𝐿) = log(Ω(𝐿)) is the microcanonical entropy we studied in Exercise 5.2. Instead of using Stirling’s formula to approximate the entropy, we will study the exact  𝑍(𝐿) and  F(𝐿) numerically.

The separation  𝑍=∑𝑍(𝐿) allows us to find the probability distribution of lengths at fixed force. Just as we studied the free energy density for the ideal gas in Section 6.7, we can use  𝑍(𝐿) to define a free energy density  Free(𝐿) for the rubber band at fixed force and coupling.

(a) As in Exercise 6.17, give the formula for  Free(𝐿)=−  log𝑍(𝐿) in terms of  𝐿,  𝛼, and  𝑆(𝐿). Write the probability  𝑝(𝐿) = 𝜌(𝐿)Δ𝐿 = 2𝜌(𝐿) of the equilibrium rubber band being of length 𝐿, in three ways. First, write it in terms of  𝑍(𝐿) and  𝑍. Then write it in terms of  Free(𝐿) and  𝑍. And finally, write it in terms of the Boltzmann-like weights  exp(−Free(𝐿’)), for all the different lengths  𝐿’.

Your answer here (or in a separate writeup)

Our ensemble fixes the force 𝐹 and the coupling  𝛼. The partition function sums the Boltzmann weight over all possible segment orientation patterns,
			Z[F,α]==∑ (e^(-H[𝐬]))==∑(Subscript[Z, F,α])[L]
 where we set  𝑘𝐵 𝑇 = 𝑘𝐵 = 1 for simplicity (measuring energy in units of  𝑘𝐵 𝑇 and entropy in ‘nats’, where  𝑘nat=1).

(b) Plot  Free(𝐿) and  𝜌(𝐿) for  𝐹=0,  𝑁=100, and with  𝛼 𝑁 = 𝑎 = 0,  0.25, ...,  1.5. (Remember that  𝜌(𝐿)= p(L) Δ L = 𝑝(𝐿)/2.) Confirm that the critical value  𝛼𝑐 at which  𝜌(𝑥) splits away from the origin is close to the continuum value of the spring constant  𝐾 in Exercise 5.12.

Clear[Ω,S,Energy, Free]
Energy[L_,α_,F_] := -... - F L
Ω[L_,Nn_]:= Nn!/(((Nn+L)/2.)!((...)!)) //N
Z[L_,α_, F_, Nn_]:= ... Exp[-Energy[...]]
Ztot[α_,F_,Nn_]:= Sum[...,{L,-Nn,Nn,2}]
S[L_,Nn_] := Log[...]
Free[L_,α_, F_, Nn_] :=...
ρ[L_,α_,F_,Nn_]:= (1/2.) .../...

Nn = 100;
F := 0
as = Table[a, {a,0,1.5,0.25}];
curves = Table[Free[L,a/Nn,F, Nn],{a,as}];
Plot[curves,{L,-Nn,Nn}]
Print["Varying α: ", as/Nn]
curves = Table[ρ[...,{a,as}];
Plot[curves,{L,-Nn,Nn}, PlotRange->All]

Your answer here (or in a separate writeup)

Note that the free energy near the transition is quantitatively similar to that of the quartic potential  (Subscript[f, 0])+(1/2) 𝑎 𝐿^2 + 𝑔 𝐿^4  as 𝑎(𝑇) passes through zero. Exercise 9.5 discusses Landau’s approach to the Ising phase transition using this quartic polynomial (eqn 9.18). He posits a quartic free energy density as a function of magnetization at fixed temperature and external field. See also Exercises 12.5 and 12.26 for other mean-field approaches to the Ising model.

(c) Plot  Free(𝐿) and  𝜌(𝐿) for  𝑎=1.25,  𝑁=100, and with a few interesting values for the force  𝐹.  (Note that for small values of F there are two stable minima. We call the higher energy minimum metastable.) At what value  𝐹𝑐 does the metastable minimum become unstable? (A rough answer for  𝐹𝑐
is fine. But if you want a precise answer, calculate the spring force  𝑓(𝐿) needed for part (d), and see where it last crosses zero as the local minimum of  Free(𝐿,𝐹) disappears at  𝐹𝑐.)

Nn = 100
Clear[F]
a = ...;
α = ...;
Fs = Table[F,...];
curves = Table[Free[...],{F,Fs}];
plotF = Plot[...]
Print["Varying F: ", Fs]

In Exercise 2, we extracted a prediction for the evolution law of the length from heat-bath dynamics. But this is not the only choice. In later chapters, we shall often assume ‘gradient’ dynamics: that the velocity is a mobility 𝛾 times minus the (variational) derivative of the free energy with respect to the “order parameter” (in this case L, see also Section 2.3). Gradient dynamics says that the tip of the rubber band evolves with the law 
				d𝐿/d𝑡 = 𝑣_gradient = 𝛾 𝑓(𝐿) = − 𝛾 𝑑(Free)/𝑑𝐿.
Here the force  𝑓(𝐿) is the force exerted by the spring when it is not in its equilibrium position. It is partly due to the external force  𝛼𝐿+𝐹 and partly due to the entropic spring force.
Let us consider the case where there is no force  𝐹 from the external world.

(d) Numerically compute the force  𝑓(𝐿) (either by finite differences or by symbolic differentiation) for  𝐹=0,  𝛼=1.25/𝑁 and  𝑁=1000. Does it go to zero at the equilibrium lengths? Compare it to the velocity of the tip of the rubber band given by heat-bath dynamics,  𝑣_HB(𝐿) = 𝑁 tanh(𝛼 𝐿) − 𝐿, derived in ‘RB Dynamics II’, by plotting
𝛾𝑓(𝐿) and  𝑣_HB(𝐿) on the same graph. Can you find a constant mobility  𝛾 that makes these two agree everywhere? ( 𝛾 is proportional to  𝑁.) Can you find a constant mobility that allows them to agree near the positive and negative equilibrium lengths? (Focus on matching slopes; a rough estimate is fine. The fixed point shifts quite a bit between  𝑁=100 for  𝑣_gradient and  𝑁=∞ for  𝑣HB; we reduce this by using  𝑁=1000.)

Clear[L,α,F,Nn]
f[L_,α_,F_,Nn_] = -D[...];

(* Mathematica's derivative goes crazy for N>100. In Python, we can't take an analytic derivative, but the discrete derivative works to N=1000. *)
Nn=100;
a=1.25;
α=a/Nn;
F=0;
vHB[L_,α_,Nn_] := ...
vHBPlot = Plot[...,{L,-Nn,Nn}];
γ=...;
Plot[{γ f[...],vHB[...]},{L,-0.95Nn,0.95Nn}]
Print["γ overall = ", γ, " does / does not match shape"]
γClose =...;
Plot[{γClose ...,...},{L,...}]
Print["γ near fixed point = ", γClose, " can / cannot match slope"]

So, we can match gradient to heat-bath dynamics locally near equilibrium by a suitable choice of mobility. This is reassuring. But they disagree in general! Is one or the other wrong? Or are they both consistent, possible dynamics that yield the same equilibrium behavior? 
The heat-bath algorithm is not an accurate representation of real rubber bands! Had we written a diffusion equation for the (efficient, but somewhat unphysical) Metropolis algorithm (Exercise 8.6, or the (grossly unphysical) Wolff algorithm (Exercise 8.8), we would have yet a different (rather strange) prediction for the velocity.

What do we need to check to see if gradient dynamics and heat-bath dynamics are both OK? Let us add fluctuations to answer this question. Again, we can compare two stochastic dynamics.

The tradition in the field is to extend gradient dynamics to Langevin dynamics by adding noise. They assume a constant  𝛾, and white noise corresponding to a fixed diffusion constant  𝐷 (see Exercises 6.18, 6.19, and 10.7). By fixing  𝐷/𝛾 = 𝑘𝐵 𝑇, they guarantee that the ensemble generated at late times is the equilibrium thermal ensemble given by the Boltzmann distribution.

Feynman, in his Lectures on Physics (vol. I sec. 43.5 near the end), derives the Einstein relation 𝐷/𝛾 = 𝑘𝐵 𝑇. He notes that the current from diffusion must cancel the current from the force due to the free energy in order for the system to be in equilibrium. He then uses the fact that the equilibrium density is given by the Boltzmann distribution. Let us consider a general free energy G(x) with equilibrium probability density 𝜌(x)exp(−G(x)/𝑘𝐵 𝑇) / 𝑍, diffusion current −𝐷 𝜌’(x), and force-driven current  𝛾𝑓(x)𝜌(x).

(e) Derive the Einstein relation,  𝐷/𝛾 = 𝑘𝐵 𝑇 by balancing the currents and using the equilibrium probability density. (It should be easier to do this on the fly than to look up Feynman’s argument, but his discussion is worth reading.)

Your answer here (or in a separate writeup)

So, do both gradient dynamics and heat-bath dynamics pass the Einstein relation test? Since Langevin dynamics uses the Einstein relation to set the noise from the damping, it certainly passes. But what about heat-bath dynamics?

In ‘RB Dynamics II’, in addition to finding  𝑣_HB(𝐿), we used the microscopic heat-bath dynamics to derive the spatially dependent diffusion constant  𝐷_HB(𝐿) = 𝑁−𝐿 tanh(𝛼𝐿). The Einstein relation then implies a spatially dependent mobility  𝛾HB(𝐿).

(f) Check numerically if  𝛾_HB(𝐿) does yield the heat-bath  𝑣_HB(𝐿), by plotting the latter along with  𝛾_HB(𝐿) 𝑓(𝐿) for  𝛼=1.25/𝑁 and  𝑁=100 (or 1000 if feasible). Is our heat-bath diffusion equation consistent with free energies and the Einstein relation? (Remember, for us kB T = 1). Discuss.

Nn = 100;
α = ...;
F = 0;
DHB[L_,α_,Nn_] := ...
γHB[L_,α_,Nn_] := ...
Plot[{...},{L,-0.95Nn,0.95Nn}]

Your discussion here (or in a separate writeup)