Rubber band dynamics II: Diffusion

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

“RB dynamics I” studied the dynamic fluctuations of an entropic model for a rubber band:  N segments of length  𝑑=1, fluctuating between pointing forward and backward at random. It studied the random walk of lengths  𝐿 as the segments hopped, both without and with an external parabolic potential stretching the band. Here we shall derive a spatially dependent diffusion equation describing the evolution of the probability distribution of lengths with time, in the limit of large  𝑁.

Current is flow forward minus flow backward. The current past the midpoint  𝐿𝑜=𝐿+1  between two possible lengths  𝐿 and  𝐿+2 is given by the probabilities  𝜌(𝐿)
  and  𝜌(𝐿+2) times the probabilities  𝑃± of flipping forward and backward. For this exercise, we assume N is even.)

Since the sum of the probabilities of being at length  𝐿,  ∑ 𝑝(𝐿)=1, is constant, and our dynamics only shifts  𝐿 locally (by  ±2≪𝑁), we are advised to write our dynamics in terms of the probability current. Note that  𝑝(𝐿)=0 for odd integers  𝐿. Let  𝐽(𝐿𝑜), for odd  𝐿𝑜 (midway between possible lengths of the chain) be the net current from  𝐿−1 to  𝐿+1 per segment flip. In “RB dynamics I”, we gave the probability  𝑃+(𝐿) per flip that a chain of length  𝐿 will grow to  𝐿+2 (contributing to  𝐽 at  𝐿𝑜=𝐿+1), and  𝑃−(𝐿) that a chain of length  𝐿 will shrink to  𝐿−2. (contributing to  𝐽 at  𝐿𝑜=𝐿−1).

(a) Our rubber band ensemble at time  𝑡 has probability  𝑝(𝐿) of having length  𝐿. Argue that the probability current J(Lo) of our rubber band ensemble growing past the (odd) length  𝐿𝑜 is
               𝐽(𝐿𝑜)Δ𝑡=𝑝(𝐿𝑜−1)𝑃+(𝐿𝑜−1)−𝑝(𝐿𝑜+1)𝑃−(𝐿𝑜+1)
where  Δ𝑡 is the time for one segment flip.

Here by convention we set  Δ𝑡=1/𝑁, so a sweep that flips every segment on average once takes one unit of time.

Your answer here (or in a separate writeup).

In taking the continuum limit as  𝑁→∞ (figure below), let us keep the total unfolded length fixed. To do so, we use  𝑥=𝐿/𝑁. Also, the harmonic stretching force  𝐹=𝛼𝐿=𝛼𝑁𝑥, so we change variables to  𝑎=𝑁𝛼. Finally, the probability  𝑝(𝐿) represents the probability density  𝜌(𝑥) between  𝐿−1 and  𝐿+1 (i.e.,  𝑥−1/𝑁 and  𝑥+1/𝑁), so we substitute  𝜌(𝑥)(2/𝑁) for  𝑝(𝐿).

Changing to continuum variables. In going from the microscopic description to the continuum limit, we change all lengths by a factor of  𝑁, we change the negative “spring” constant  𝛼 to  𝑎=𝑁𝛼, and we change from probabilities  𝑝(𝐿) to probability densities  𝜌(𝑥)=𝑝(𝐿)𝑁/2.

(b) Substituting 𝐿=𝑁𝑥, 𝛼=𝑎/𝑁, and 𝑃±, show that
        𝐽(𝐿)=[𝑒^(𝑎𝑥) (1−𝑥+1/𝑁)𝜌(𝑥−1/𝑁) − 𝑒^(−𝑎𝑥) (1+𝑥+1/𝑁)𝜌(𝑥+1/𝑁)] / [𝑒^(𝑎𝑥)+𝑒^(−𝑎𝑥)]

(* From 'RB Dynamics I', code also given below *)
energy[L_,α_] := -(1/2) α L^2;
Pplus[L_,N_,α_] := (N-L)/(2 N) (1/(1+Exp[energy[L+2, α]-energy[L,α]]))
(* Flips of up spins to down, from length L to L-2 *)
Pminus[L_,N_,α_] := (N+L)/(2N)(1/(1+Exp[energy[L-2,α]-energy[L,α]]))

J[Lo_,N_,α_]:= (1/Δt)( ...)
Δt = ...;
JofX[x_,N_,a_] = J[Lo,N,α] /.p[L_]->.../N/. {Lo->...,...}//Simplify

(c) What is the net current  𝐽 in the limit  𝑁→∞ (holding  𝑎 and  𝑥 constant)? Use the fact that  (𝑒^(𝑎𝑥)−𝑒^(−𝑎𝑥))/(𝑒^(𝑎𝑥)+𝑒^(−𝑎𝑥))=tanh(𝑎𝑥) to remove the exponentials from your answer, and show your work. Argue that the end of the rubber band as  𝑁→∞ has a position-dependent velocity 
                      𝑣(𝑥)=tanh(𝑎𝑥)−𝑥.
In ‘RB Dynamics I’(e),  you found that for large a the length quickly moved to
a final off-center position. Find a  numerical solution for this final length at  𝑎=0.75 and  𝑎=1.5 for  𝑁→∞. Derive from  𝑣(𝑥) the critical value 𝑎𝑐 when the rubber band equilibrium length splits away from the origin.

JNinf = ExpToTrig[Normal[Series[...,{N,∞,...}]]]//Simplify
vHB[x_,a_] = ...;

Your answer here.

Xmin[a_]:=  x/. FindRoot[...==0,{x,1}]
Print["Final L for a=1.5 is N times ",...] 
Print["Final L for a=0.75 is N times ",...]

(d) Using your random-walk simulation from “RB Dynamics I”, make a histogram of lengths for  𝑁=100 at  𝑎=0.75 and  𝑁=1000 at  𝑎=1.5, adding points until you get good histograms. When necessary, drop the transient first part of the trajectory, while the rubber band moves from zero to the new minimum. Are your histograms concentrated near the predicted value you found in part (c)?

We supply the complete code here. You may wish to substitute your answer from RB Dynamics I, or use this to help solve the previous problem.

energy[L_,α_] := -(1/2) α L^2;
(* Flips of down spins to up, from length L to L+2. The rate of flipping from L to right is Pplus[L] *)
Pplus[L_,N_,α_] := (N-L)/(2 N) (1/(1+Exp[energy[L+2, α]-energy[L,α]]))
(* Flips of up spins to down, from length L to L-2 *)
Pminus[L_,N_,α_] := (N+L)/(2N)(1/(1+Exp[energy[L-2,α]-energy[L,α]]))
HeatBathFlip[L_,N_,α_] := Block[{r = Random[]},If[r<Pplus[L,N,α], L+2, If[r<Pplus[L,N,α]+Pminus[L,N,α], L-2, L]]]

n = 100;
a = 0.75;
α = a/n;
...
trajL = NestList[HeatBathFlip[#,n,α]&, ...];
trajectory = ...;
trajPlot[a] = ListPlot[...]
histPlot[a] = Histogram[...]

n = 1000;
a = ...
...
histPlot[a] = Histogram[Drop[trajL,10n], {Table[m,{m,-n,n,2}]}, "PDF"]

Your answer here.

(e) Starting at  𝑥=0, launch a trajectory for  𝑁=1000 and  𝑎=1.5, and examine how it flows to its final value once it deviates from the local fixed point at zero (say, in the first 20000 segment flips). Compare the flow to that predicted by your equation for the velocity as a function of length in part (c). (That is, numerically solve  𝑑𝑥/𝑑𝑡=𝑣(𝑥) from part (c), starting from a small positive or negative value of  𝑥, and rescale it from  𝑥 and 𝑎 to 𝐿 and  𝛼.) The time spent in the vicinity of  𝐿=0 in the random walk will depend on the fluctuations. Adjust the theory curve right and left to make a good comparison.

a = 1.5;
xSol[a] = x/.NDSolve[{x'[t]==..., x[0]==0.001},x,{t,0,40}][[1]];

Plot[xSol[a][t],{t,...,...}]

(* We've stored the trajectory plots as trajPlot[a],so that we can vary the time lag and sign of the theory to match it without regenerating new curves. *)

a = 1.5;
n = 1000;
dt = ...;
theoryPlot = Plot[-...,{t,-dt,20-dt}, PlotStyle->Green];
Show[trajPlot[a],theoryPlot]

The current in the limit  𝑁→∞ is the answer that the continuum limit supplies. Thermodynamics and other continuum theories often ignore the fluctuations in the system. We can study the statistical mechanics of the fluctuations by studying the leading corrections in  1/𝑁.

(f) Now find the first correction in  1/𝑁 to the net current. Write the current to this order in the form
				𝐽(𝑥)=(𝑣(𝑥)+𝑣1(𝑥)/𝑁) 𝜌(𝑥) − 𝐷(𝑥) ∂𝜌/∂𝑥.
Calculate  𝑣1(𝑥). Show that
				𝐷(𝑥)=(1−𝑥 tanh(𝑎𝑥))/𝑁.

Clear[a]
J2ndOrder = ExpToTrig[Normal[Series[...]]]//Simplify

NewTerms = J2ndOrder - JNinf // Simplify
Diff[x_,a_]= ...
v1HB[x_,a_] = ...

The corrections  𝑣1(𝑥) make tiny corrections to the fixed point of part(d) and the velocity curve of part(e) above -- unimportant for large  𝑁. But the term proportional to  ∂𝜌/∂𝑥, the typical domain of statistical mechanics, dominates many of the properties.)

(g) Write the emergent forced diffusion equation governing our entropic rubber band. (To simplify things, leave it in terms of  𝑣,  𝐷,  𝜌, and their derivatives. Ignore the terms involving  𝑣1.)

j[x_,t_]:= ...
D[ρ[x,t],t] ==...

(h) Approximate your differential equation to linear order in 𝑥 about zero (again ignoring  𝑣1), and solve your diffusion equation for the stationary distribution at  𝑎=0.75. (Hint: Solving for the distribution that makes the current equal to zero is easier. Try a Gaussian.) Compare to a histogram of equilibrated values of your random walk simulation.

Clear[diff,v,a,σSol]
diff[a_] = Normal[Series[..., {x,0,1}]]
v[a_] = ...
current[a_, ρ_] := ...
ρStationary[x_,σ_] = (1/Sqrt[...]) Exp[-...];
current[a,ρStationary[x,σ]]//Simplify
σSol[a_,N_] = ...
n = 100;
stationaryPlot = Plot[ρStationary[x,n σSol[0.75,n]], {x,-n,n}];
Show[histPlot[0.75], stationaryPlot]