140B Statistical Physics

Part 3. Phase Transitions

Liquid-Gas Transition

Phases and Phase Transitions

States of Matter

Grid[{{"","Fluid",SpanFromLeft,"Solid"},{"","Gas","Liquid",SpanFromAbove},
Block[{atom},
atom[r_]:=Disk[r,Scaled[0.05]];
SeedRandom[17348];
Prepend[Graphics[{Line[{{0,0},{0,5},{5,5},{5,0},{0,0}}],atom/@#},PlotRangePadding->0.1,PlotRange->{{0,5},{0,5}},ImageSize->{50,50}]&/@{{{1,1},{3,4}},(Tuples[Range[4],2]+RandomVariate[NormalDistribution[0,0.4],{16,2}]),Tuples[Range[4],2]},""]],
{"Atom arrangement","Random","Random","Periodic"},{"Traslation symmetry","Preserved","Preserved","Broken"},{"Correlation","Negligible","Short\[Hyphen]range","Long\[Hyphen]range"},{"Density","Low","High","High"},{"Compressibility","High","Low","Low"},{"Fluidity","✓","✓","×"},{"Sound modes","1","1","3"}},Dividers->{{False,{True},False},{3->True}},ItemStyle->"Text"]

The key difference lies in the density correlation

C(r)=⟨n(r)n(0)⟩,

which quantifies how the presence of an atom at the origin 0 affects the likelihood of finding another atom at some point r away from the origin.

Block[{C},
C[A_,ξ_,b_][x_]:=Exp[A Cos[2π x] Exp[-x/ξ]-b];
Row[Plot[#1[x],{x,0,7.5},PlotStyle->Black,Filling->Axis,PlotRange->{0,3},PlotLabel->Style[#2,"Text"],FrameLabel->{"\!\(r/a\)","\!\(C(r)\)"},ImagePadding->{{#3,2},{42,0}},AspectRatio->0.8,ImageSize->{Automatic,140}]&@@@{{C[3,0.2,1.5],"Gas",42},{C[3,1,0],"Liquid",2},{C[3,∞,0],"Solid",2}}]]

Gas: atoms are far apart, almost independent, no regular arrangement. Correlation is negligible (very weak and short-ranged).

Liquid: atoms are close together, sliding around each other, no regular arrangement. Correlation is strong but short-ranged.

Solid: atoms are closely packed, arranging in regular patterns (crystal lattices). Correlation is strong and long-ranged.

Different microscopic order of atoms leads to different macroscopic properties of the matter.

Phase Diagram

Block[{lt=0.8,col,bdy},
col=<|"g"->White,"l"->Blue,"s"->Green|>;
bdy=<|"sg"->BSplineCurve[{{0,0},{0.1,0},{0.6,0.3},{1,1}}],"lg"->BSplineCurve[{{1,1},{2.5,1.5},{3,3}}],"sl"->Line[{{1,1},{0.4,4}}]|>;
Graphics[{{Lighter[col["s"],lt],FilledCurve[{bdy["sg"],bdy["sl"],Line[{{0.4,4},{0,4}}]}]},
{Lighter[col["l"],lt],FilledCurve[{bdy["lg"],Line[{{3.3,4},{0.4,4}}]}]},
Block[{dx=0.1,h=1,s=0.3,w=0.4},Translate[{EdgeForm[None],Table[{Blend[Lighter[col[#],lt]&/@{"l","g"},(x+1)/2],Triangle[{{0,0},{s+w x,h},{s+w(x+dx),h}}]},{x,-1,1-dx,dx}]},bdy["lg"][[1,-1]]]],
{Opacity[0.2],AbsoluteThickness[1],Line[{{0,3},{3,3},{3,0}}],Line[{{{0,2},{4,2}},{{0.8,0},{0.8,2}},{{2.5,2},{2.5,0}}}],Line[{{0,1},{1,1},{1,0}}]},
{AbsoluteThickness[1.5],bdy["sg"],bdy["lg"],bdy["sl"]},{AbsolutePointSize[5],Point[{bdy["lg"][[1,1]],bdy["lg"][[1,-1]]}]},
{Text[Style["Ice\n(solid)",LineSpacing->{0, 12}],{0.5,1.3}],Text[Style["Water\n(liquid)",LineSpacing->{0,12}],{1.8,2.5}],Text[Style["Vapor\n(gas)",LineSpacing->{0,12}],{2.5,0.8}],
Text[Style["Supercritical\nfluid",12,LineSpacing->{0,10}],{3,3.5}],Text[Style["Triple\npoint",10,LineSpacing->{0,8}],{1,1},{-1,1.2}],
Text[Style["Critical\npoint",10,LineSpacing->{0,8}],{3,3},{-1,1.2}]}},PlotRange->{{0,4},{0,4}},Frame->True,FrameTicks->{{{{1,"0.6"},{2,"101"},{3,"\!\(2.2×10\^4\)"}},{{2,"1 atm"}}},{{{0.8,"0"},{1,""},{1.2,Rotate["0.01",-π/4],{0,0}},{2.5,Rotate["100",-π/4]},{3,Rotate["374",-π/4]}},None}},FrameLabel->{"\!\( T\) [°C]","\!\( P\) [kPa]"},ImageSize->250]]

Each macrostate of a system is labeled by a set of coordinates g=((Subscript[g, 1]),(Subscript[g, 2]),…) in the phase diagram. For example, g=(T,P) in the above phase diagram.

Phase transition happens across the phase boundary, where the free energy density

f(g)=(lim_(V→∞))((F(g))/V)

becomes non-analytic with respect to g in the thermodynamic limit, i.e. its derivatives become singular or discontinuous.

Phase: Two macrostates g and g' belongs to the same phase as long as there exists a continuous path from g to g' without encountering any phase transition.

Liquid and gas are in the same phase (the fluid phase).

Even though liquid and gas are separated by a phase transition, it does not make them different phases of matter, because the phase transition can be circumvented through the supercritical fluid regime.

This is the only logically consistent way of defining phases: 
they should be unified by the absence of phase transitions, 
other than separated by the presence of phase transitions.

Row[Graphics[{{Lighter[Orange,0.8],FilledCurve[{#1,Line[{{0,1},{-1,1},{-1,-1},{0,-1}}]}]},{Lighter[Blue,0.8],FilledCurve[{#1,Line[{{0,1},{1,1},{1,-1},{0,-1}}]}]},#1,{Point[{0,#}&/@{-0.5,0.5}],Text["\!\( g\)",{0.1,-0.5},{-1,0}],Text["\!\( g\^′\)",{0.1,0.5},{-1,0}]},{Opacity[0.2],Arrowheads[{-0.1,0.1}],Red,Arrow[{0,#}&/@{-0.5,0.5}](*,Green,Arrow[BSplineCurve@{{0,-0.5},{0.9,-0.2},{0.9,0},{0.9,0.2},{0,0.5}}]*)}},PlotRangePadding->0,Frame->True,FrameTicks->None,ImageSize->80]&@@@{{BSplineCurve[#+{-0.4,0}&/@{{0,-1},{0,-0.3},{0.5,-0.2},{1,0},{0.5,0.2},{0,0.3},{0,1}}]},{BSplineCurve[#+{-0.4,0}&/@{{0,-1},{0,-0.3},{0.2,0},{0.4,0},{1.2,0},{0.4,0},{0.2,0},{0,0.3},{0,1}}]}},Spacer[1]]

Solid is in a different phase from liquid and gas.

There can be many different solid phases characterized by different ways of breaking the translation and rotation symmetry of the space.

Lattice Gas Model

Model Hamiltonian

Assuming the space is partitioned into many small cells (or lattice sites), each labeled by an index, denoted as i, j.

Block[{L=4},SeedRandom[1];Graphics[{{LightGray,AbsoluteThickness[1],Rotate[Translate[Line[{#,0}&/@{-L,L}],{0,#}]&/@Range[-L,L],# π/2,{0,0}]&/@{0,1}},{Black,Disk[#,0.4]&/@RandomChoice[Range[-L+1/2,L-1/2],{10,2}]},{White,Text["i",{-1/2,1/2}]},{Black,Text["j",{1/2,1/2}]}},ImageSize->140]]

The microstate of the system is described by a collection of binary random variables n={(Subscript[n, i])}, where

describes the number of atoms on site i, assuming each site can not host more than one atom (the hard-core condition).

g>0 is the attractive interaction strength between atoms: it provides an energy reward if atoms are neighboring to each other.

⟨i j⟩ denotes the summation over pairs of neighboring sites i and j.

The energy -g (Subscript[n, i])(Subscript[n, j]) under different scenarios:

({{(Subscript[n, i]), (Subscript[n, j]), -g (Subscript[n, i])(Subscript[n, j]),  }, {0, 0, 0, (no energy reward)}, {0, 1, 0, }, {1, 0, 0, }, {1, 1, -g, (attractive energy reward)}})

μ∈ℝ is the chemical potential, s.t. each atom experiences a potential energy of -μ to occupy a site.

It is introduced to tune the total particle number

To lower the energy, N tends to increase with increasing  μ.

Lattice Structure

Block[{col1,col2,site,a=0.6,z={0.4,0.3}},
{col1,col2}={Gray,Black};
site=Block[{r=0.12},
{EdgeForm[Directive[col2,AbsoluteThickness[1.2]]],White,Disk[{0,0},r]}];
Grid[{{"1D","2D",SpanFromLeft,"3D"},
Graphics[#,PlotRange->1.5,ImageSize->{60,60}]&/@{{{col1,Line[{#,0}&/@{-1.4,1.4}]},{Translate[site,{#,0}]&/@Range[-1,1]}},{{col1,Rotate[Translate[Line[{#,0}&/@{-1.4,1.4}],{0,#}]&/@Range[-1,1],# π/2,{0,0}]&/@{0,1}},{Translate[site,#]&/@Tuples[Range[-1,1],2]}},{{col1,Rotate[{Line[{#,0}&/@{-1.5,1.5}],Translate[Line[{#,0}&/@{-1,1}],{0,Sqrt[3]/2 #}]&/@{-1,1}},# 2π/3,{0,0}]&/@Range[3]},{Translate[site,#]&/@Cases[Tuples[Range[-2,2],2].{{1,0},{-1/2,Sqrt[3]/2}},r_/;Norm[r]<=1]}},{{col1,Translate[Rotate[Translate[Line[{#,0}&/@{-1,1}],{0,#}]&/@{-a,a},# π/2]&/@{0,1},z a #]&/@{-1,1},Translate[Line[z #&/@{-1,1}],RotationMatrix[# π/2].{a,a}]&/@Range[4]},{Translate[site,#.{{1,0},{0,1},z}]&/@Tuples[{-a,a},3]}}},
{"chain","sqaure","triangle","cubic"},{"\!\(q=2\)","\!\(q=4\)","\!\(q=6\)","\!\(q=6\)"}},ItemStyle->"InlineFormula",Dividers->{{False,{True},False},{False,True}}]]

Coordination number q: the number of nearest-neighbor sites around each given site, reflecting how connected or coordinated each site is within the lattice.

Volume: Number of sites on the lattice.

Number of nearest-neighbor bonds (links) on the lattice.

Probability Distribution

Z=(∑_n)(e^(-β E(n))).

High temperature limit (T→∞,  β->0): p(n)=1/Z reduces to a uniform distribution of n. In this case, each site will be occupied or empty with 1/2 to 1/2 probability => Supercritical fluid (an uniform mixture of liquid and gas).

Block[{L=30,cols={White,Black},n},
SeedRandom[0];
n=RandomChoice[{0.5,0.5}->{0,1},{L,L}];
ArrayPlot[n,ColorRules->Thread[{0,1}->cols],PlotRange->{0,1},PlotRangePadding->0,PlotLabel->Style[StringTemplate["\!\(N/V=`1`\)"][NumberForm[Total@Flatten@n/L^2,{3,2}]],Black,14],ImageSize->108]]

Low temperature limit (T→0,  β->∞): p(n) will be dominated by the lowest-energy configuration n.

A over-simplified estimation: assuming the same particle number on every site, i.e. ∀i:(Subscript[n, i])=n, the energy function becomes

where V is the lattice volume and q is the coordination number (assuming square lattice with q=4).

Derive Eq. ().

E(n)=-g(∑_(⟨i j⟩))(Subscript[n, i])(Subscript[n, j])-μ(∑_i)(Subscript[n, i])
=-g(∑_(⟨i j⟩))(n^2)-μ(∑_i)n
=-g (n^2)(∑_(⟨i j⟩))1-μ n(∑_i)1.

Number of sites on the lattice is

(∑_i)1=V.

Number of bonds on the lattice is

(∑_(⟨i j⟩))1=(q/2)V.

E(n)=-(q/2)g (n^2)V-μ n V
=-V n ((q/2) g n+μ).

Block[{q=4},
Plot[Evaluate@Table[-n(q n/2+μ),{n,0,1}],{μ,-4,0},PlotStyle->{Gray,Black},PlotRange->{-3,3},FrameLabel->{"\!\(μ/g\)","\!\(E(n)/g V\)"},GridLines->{{-2},None},Epilog->{{Gray,Text["\!\(n=0\)",##]&@@@{{{0,0},{1,-1}},{{-4,0},{-1,1}}}},{Black,Text["\!\(n=1\)",##,{1,-0.4}]&@@@{{{0,-2},{1,1}},{{-4,2},{-1,-1}}}}},Axes->False,ImageSize->180]]

When  μ<-2 g, the configuration n=0 is lower in energy. => Gas (in the dilute limit). Atoms are far apart from each other — so dilute that we can’t see any atom on a finite lattice.

Block[{L=30,cols={White,Black},n},
SeedRandom[0];
n=ConstantArray[0,{L,L}];
ArrayPlot[n,ColorRules->Thread[{0,1}->cols],PlotRange->{0,1},PlotRangePadding->0,PlotLabel->Style[StringTemplate["\!\(N/V=`1`\)"][NumberForm[Total@Flatten@n/L^2,{3,2}]],Black,14],ImageSize->108]]

When  μ>-2 g, the configuration n=1 is lower in energy. => Liquid (in the dense limit). Atoms are closely packed, next to each other — so dense that we can’t see any bubble on a finite lattice.

Block[{L=30,cols={White,Black},n},
SeedRandom[0];
n=ConstantArray[1,{L,L}];
ArrayPlot[n,ColorRules->Thread[{0,1}->cols],PlotRange->{0,1},PlotRangePadding->0,PlotLabel->Style[StringTemplate["\!\(N/V=`1`\)"][NumberForm[Total@Flatten@n/L^2,{3,2}]],Black,14],ImageSize->108]]

Markov Chain Monte Carlo

Monte Carlo Sampling

start with an arbitrary initial sample (x^((0)))\[Distributed](Subscript[p, 0])((x^((0)))),

generate a new sample (x^((1))) following the transition probability p((x^((1)))|(x^((0)))) (to be designed later)

(x^((0)))→(x^((1)))\[Distributed]p((x^((1)))|(x^((0)))),

iterate the process to construct a chain of samples

(x^((0)))→(x^((1)))->(x^((2)))->…→(x^((n)))->….

Markov Chain

The chain of samples Eq. () generated by the Markov process is a Markov chain, along which the sample probability evolves as

(Subscript[p, n])((x^((n))))=(∑_(x^((n-1))))p((x^((n)))|(x^((n-1))))(Subscript[p, n-1])((x^((n-1)))).

The evolution is uniquely determined by the transition probability p(x'|x).

One may start with any initial distribution (Subscript[p, 0])(x). Under the Markov process, the distribution (Subscript[p, n])(x) might evolve and relax to a stationary distribution p(x)

(lim_(n→∞)) (Subscript[p, n])(x)=p(x).

Detailed Balance

Argument: detailed balance implies that, for every pair of states , the flow of probability is perfectly balanced between the forward sampling x→x' and the reversed sampling x'→x, such that the probability will not evolve in time, i.e. stationary.

Alternatively, Eq. () can be written as

((p(x'|x))/(p(x|x')))=((p(x'))/(p(x))).

Note: for x'=x, Eq. () is trivially satisfied, thus one only need to check the detailed balance for the case of x'≠x.

Metropolis-Hastings Algorithm

Given the current sample x.

Propose a new sample (OverTilde[x]) by drawing from the proposal probability q((OverTilde[x])|x).

Calculate the acceptance probability

Accept or Reject:

with probability α((OverTilde[x]),x), accept the new sample and set x'=(OverTilde[x]),

with probability 1-α((OverTilde[x]),x), reject the new sample and set x'=x.

Return x' as the updated sample.

p(x'|x)=q(x'|x)α(x',x),

Show that the transition probability in Eq. () satisfies the detailed balance condition in Eq. ().

((p(x'|x))/(p(x|x')))=((q(x'|x)α(x',x))/(q(x|x')α(x,x')))
=((q(x'|x))/(q(x|x')))((min(1,((p(x')q(x|x'))/(p(x)q(x'|x)))))/(min(1,((p(x)q(x'|x))/(p(x')q(x|x'))))))
=((q(x'|x))/(q(x|x')))((min(1,a(x',x)))/(min(1,(1/(a(x',x))))))

a(x',x)=((p(x')q(x|x'))/(p(x)q(x'|x))).

if a(x',x)>1,

min(1,a(x',x))=1, 
min(1,(1/(a(x',x))))=(1/(a(x',x))),

therefore

((p(x'|x))/(p(x|x')))=((q(x'|x))/(q(x|x')))(1/(1/a(x',x)))=((q(x'|x))/(q(x|x')))a(x',x).

if a(x',x)<1,

min(1,a(x',x))=a(x',x), 
min(1,(1/(a(x',x))))=1,

therefore

((p(x'|x))/(p(x|x')))=((q(x'|x))/(q(x|x')))((a(x',x))/1)=((q(x'|x))/(q(x|x')))a(x',x).

((p(x'|x))/(p(x|x')))=((q(x'|x))/(q(x|x')))a(x',x)
=((q(x'|x))/(q(x|x')))((p(x')q(x|x'))/(p(x)q(x'|x)))
=((p(x'))/(p(x))),

q(x|(OverTilde[x]))=q((OverTilde[x])|x),

if the proposed sample has a lower (or equal) energy, i.e. E((OverTilde[x]))<=E(x), the proposal is always accepted,

if the proposed sample has a higher energy, i.e. E((OverTilde[x]))>E(x), the proposal is accepted with probability (e^(E(x)-E((OverTilde[x])))) that decays exponentially with the energy difference.

Monte Carlo Simulation

E(n)=-g(∑_(⟨i j⟩))(Subscript[n, i])(Subscript[n, j])-μ(∑_i)(Subscript[n, i]).

Start with (arbitrary) initial configuration n.

Randomly choose a site i.

Propose a new configuration n' by flipping the occupation number on site i, i.e.

(Subscript[n, i]^′)=1-(Subscript[n, i]),
(Subscript[n, j]^′)=(Subscript[n, j])  (for all j≠i),

Compute the energy difference

Δ E=E(n')-E(n).

If the energy lowers, Δ E<0, accept the update n->n',

If the energy raises, Δ E>0, accept the update n->n' with probability

α=(e^(-β Δ E)),

and reject with probability 1-α (if rejected, configuration will remain in n).

Repeat to collect a series of configurations n.

DynamicModule[{L=50,cols={White,Black},rate=0,n,update,Tμ,bg},
SeedRandom[0];
n=RandomChoice[{0.5,0.5}->{0,1},{L,L}];
update[{T_,μ_}][n_]:=Block[{i,js,δE},
i=RandomInteger[{1,L},{2}];
js=Mod[i+#,L,1]&/@{{1,0},{-1,0},{0,1},{0,-1}};
δE=(1-2Extract[n,i])(-Total[Extract[n,#]&/@js]-μ);
RandomChoice[{1-#,#}&@Min[Exp[-δE/T],1]->{n,ReplacePart[n,{i->1-Extract[n,i]}]}]];
Tμ={1.,-2};
n=Nest[update[Tμ],n,10000];
bg=Block[{μc=-2,Tc=1/(2 Log[1+Sqrt[2]]),m,lt=0.5},
m[{T_,μ_}]:=FixedPoint[Tanh[(μ-μc+Tc #)/T]&,0.,100];Graphics[{{Thread@{Lighter[#,lt]&/@cols,Rectangle[{0,-2},{0.35,#}]&/@{-3,-1}},
First@ContourPlot[(m[{T,μ}]+1)/2,{T,0.1,1.2},{μ,-3,-1},PlotRange->{0,1},Contours->Range[-0.05,1.05,0.1],ContourStyle->None,ColorFunction->(Blend[Lighter[#,lt]&/@cols,#]&),ColorFunctionScaling->False]},{AbsoluteThickness[2],Line[{#,μc}&/@{0,Tc}]},{AbsolutePointSize[5],Point[{Tc,μc}]}},PlotRange->{{0,1.2},{-3,-1}},Frame->True,FrameLabel->{"\!\( T/g\)","\!\(μ/g\)"},AspectRatio->1,ImageSize->180]];
Column[{Row@{Style["Update rate: ","Text"],SetterBar[Dynamic[rate],{0,1,10,200,5000}]},
Row[
{LocatorPane[Dynamic[Tμ],
Column[{Style[Dynamic[StringTemplate["\!\((T,μ)/g=(`1`,`2`)\)"]@@(NumberForm[#,{3,2}]&/@N[Tμ])],"InlineFormula"],bg},Alignment->Right],{{0.1,-3},{1.2,-1}}],
Dynamic[Refresh[n=Nest[update[Tμ],n,rate];ArrayPlot[n,ColorRules->Thread[{0,1}->cols],PlotRange->{0,1},PlotRangePadding->0,PlotLabel->Style[StringTemplate["\!\(N/V=`1`\)"][NumberForm[Total@Flatten@n/L^2,{3,2}]],Black,14],ImageSize->180],UpdateInterval->0.05],TrackedSymbols:>{n,Tμ}]},Spacer[1]]}]]

Magnetic Transition

Ising Model

Model Hamiltonian

The microscopic state is describes by a set of Ising spins s={(Subscript[s, i])}, where

describes the spin state on the lattice site i.

J>0 is the coupling constant, representing the interaction strength between neighboring spins.

⟨i j⟩ denotes the summation over pairs of neighboring sites i and j.

The energy -J (Subscript[s, i])(Subscript[s, j]) under different scenarios:

({{spins, (Subscript[s, i]), (Subscript[s, j]), -J (Subscript[s, i])(Subscript[s, j]), energy}, {aligned, +1, +1, -J, gain}, {anti-aligned, +1, -1, +J, cost}, {anti-aligned, -1, +1, +J, cost}, {aligned, -1, -1, -J, gain}})

h∈ℝ is the external magnetic field, s.t. the spin gains energy h by aligning with the external magnetic field, and costs energy h to anti-align.

It is introduced to tune the total magnetization

The magnetization M is energetically favored to align with the external magnetic field h (i.e. to be of the same sign).

Ising Symmetry

Mapping between on-site binary variables.

Grid[{{"Lattice gas","\!\(n\_i=0\)","\!\(n\_i=1\)"},Prepend[Graphics[{{LightGray,AbsoluteThickness[1],Line[{{-1,-1},{-1,1},{1,1},{1,-1},{-1,-1}}]},{#,Disk[{0,0},0.8]}},ImageSize->16]&/@{White,Black},SpanFromAbove],{SpanFromAbove,"empty","occupied"},{"Ising","\!\(s\_i=+1\)","\!\(s\_i=-1\)"},
Prepend[Graphics[{Arrowheads[{{0.01,1,Graphics[Line[Offset/@{{-5,3},{0,0},{-5,-3}}]]}}],<|1->Red,-1->Blue|>[#],Arrow[#{{0,-0.8},{0,0.9}}]},PlotRange->1,ImageSize->16]&/@{1,-1},SpanFromAbove],{SpanFromAbove,"spin up","spin down"}},ItemStyle->"InlineFormula",Dividers->{{2->True},{4->True}},Alignment->{Center,Center}]

Using Eq. (), the lattice gas model

E(n)=-g(∑_(⟨i j⟩))(Subscript[n, i])(Subscript[n, j])-μ(∑_i)(Subscript[n, i]),

is equivalent to the Ising model (up to a constant energy shift)

E(s)=-J(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j])-h(∑_i)(Subscript[s, i]).

with the following relations

Derive the relations in Eq. ().

E=-g(∑_(⟨i j⟩))((1+(Subscript[s, i]))/2)((1+(Subscript[s, j]))/2)-μ(∑_i)((1+(Subscript[s, i]))/2)
=-(g/4)(∑_(⟨i j⟩))(1+(Subscript[s, i])+(Subscript[s, j])+(Subscript[s, i])(Subscript[s, j]))-(μ/2)(∑_i)(1+(Subscript[s, i]))
=-(g/4)((∑_(⟨i j⟩))1+(∑_(⟨i j⟩))((Subscript[s, i])+(Subscript[s, j]))+(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j]))-(μ/2)((∑_i)1+(∑_i)(Subscript[s, i])).

(∑_i)1=V, (∑_(⟨i j⟩))1=(q/2)V.

(∑_(⟨i j⟩))(Subscript[s, i])=(1/2)(∑_i)(∑_(j<->i))(Subscript[s, i])=(1/2)(∑_i)(Subscript[s, i])(∑_(j<->i))1=(q/2)(∑_i)(Subscript[s, i]),

E=-(g/4)(V+q(∑_i)(Subscript[s, i])+(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j]))-(μ/2)(V+(∑_i)(Subscript[s, i]))
=-(g/4)(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j])-(((q g+2μ)/4))(∑_i)(Subscript[s, i])-(((g+2μ)V)/4).

The Ising model has an explicit (Subscript[ℤ, 2]) Ising symmetry at h=0 (in the absence of external magnetic field).

The symmetry acts by flipping all Ising spins together,

The symmetry transformations a (Subscript[ℤ, 2]) group — the two-fold cyclic group containing two elements: identity (do nothing) and spin flip, s.t.

flip+flip=identity.

The energy function (at h=0) is invariant under the symmetry transformation

E(s)=-J(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j]).

Probability Distribution

Z=(∑_s)(e^(-β E(s))).

Combinatorial complexity: With V spins, there are (2^V) possible configurations. This exponential growth in configurations makes exact summation in Eq. () intractable for large systems, especially in higher dimensions.

Correlation complexity: The correlated nature of spins rooted in the interactions among them, where the state of one spin can influence the probable states of others. This preclude simplifying the analysis to isolated spins, necessitating a collective analysis.

d - spatial dimension of lattice,

h - external magnetic field.

Mean Field Theory

General Idea

Key assumption: neglect the correlation between Ising spins, assuming each spin only feels an average interaction effect from other spins, hence the name “mean field”.

Introduce average magnetization density m,

Block[{pts={{0,0},{1,0},{0,1},{-1,0},{0,-1}},cross,w=1.5},cross[mf_:False]:={{Gray,AbsoluteThickness[1],Rotate[Line[{#,0}&/@{-1,1}],# π/2,{0,0}]&/@{0,1}},{White,Disk[#,0.2]&/@pts},Text@@@Thread@{If[mf,{"\!\(s\_0\)"}~Join~ConstantArray["\!\( m\)",4],StringTemplate["\!\(s\_`1`\)"]/@Range[0,4]],pts}};
Graphics[{Translate[cross[#1],{#2 w,0}]&@@@{{False,-1},{True,1}},Text["\[Implies]",{0,0}]},ImageSize->200]]

This enables two crucial simplifications:

Independence: Each spin (Subscript[s, 0]) is treated as if it is independent of the others (Subscript[s, i≠0]), under the influence of the mean field m only.

Factorization: The joint probability distribution p(s) is approximated by a product of independent distributions for each spin.

Key steps:

Probability Factorization: propose a factorized probability model (Subscript[p, MF])(s) to approximate the ground truth p(s).

Objective Function: formulate an objective function that measure the divergence between (Subscript[p, MF])(s) and p(s).

Parameter Optimization: optimize model parameters in (Subscript[p, MF])(s) to minimize its distance to p(s), such that (Subscript[p, MF])(s) approximates p(s) as much as possible (within the model capacity).

Model Prediction: use the model (Subscript[p, MF])(s) to predict whatever thermodynamic properties of interest.

Probability Factorization

(Subscript[p, MF])(s)=(∏_i)(Subscript[p, MF])((Subscript[s, i])).

Normalization is automatically satisfied

(∑_((Subscript[s, i])=±1))(Subscript[p, MF])((Subscript[s, i]))=1.

The model parameter m has a clear physical meaning — the expectation value of (Subscript[s, i]) (i.e. the magnetization density)

⟨(Subscript[s, i])⟩=(∑_((Subscript[s, i])=±1))(Subscript[s, i])(Subscript[p, MF])((Subscript[s, i]))=m.

Show that Eq. () is the unique solution of Eq. ().

(Subscript[p, MF])((Subscript[s, i]))=({{{, ({{q, (Subscript[s, i])=+1,}, {1-q, (Subscript[s, i])=-1,}})}})

⟨(Subscript[s, i])⟩=(∑_((Subscript[s, i])=±1))(Subscript[s, i])(Subscript[p, MF])((Subscript[s, i]))
=(+1)q+(-1)(1-q)
=2q-1.

2q-1=m,

q=((1+m)/2),

1-q=((1-m)/2).

Objective Function

Kullback-Leibler (KL) Divergence

The term “divergence” is used instead of “distance” because (Subscript[D, KL]) is not symmetric in general, i.e. (Subscript[D, KL])(p||q)!=(Subscript[D, KL])(q||p), therefore not qualified as a distance.

Non-negativity of KL divergence: for any two distributions p(x) and q(x),

(Subscript[D, KL])(p||q)≥0.

The equality is achieved iff p(x) and q(x) are identical distributions, i.e.

(Subscript[D, KL])(p||q)=0⇔∀x:p(x)=q(x).

Prove Eq. () that the KL divergence is non-negative.

Plot[{Log[x],x-1},{x,0,4},Epilog->{Blue,Text["\!\(log x\)",{3.5,Log[3.5]},{0,1}],Red,Text["\!\(x-1\)",{2.5,2.5},{1,1}]},FrameLabel->{"\!\( x\)",None},ImageSize->160]

(Subscript[D, KL])(p(x)||q(x))=-(∑_(x∈𝒳))p(x)log ((q(x))/(p(x)))
≥-(∑_(x∈𝒳))p(x)(((q(x))/(p(x)))-1)
=-(∑_(x∈𝒳))q(x)+(∑_(x∈𝒳))p(x)
=-1+1=0.

Reality p(x): The world is full of randomness. Let p(x) describe the probability of event x happening in reality.

Belief q(x): When we observe reality, our brain unconsciously construct a (different) probability model q(x) in our mind, representing our belief about the world.

Surprise: Sometimes, we witness an event x that, according to our belief q(x), was highly unlikely, yet it actually occurs in reality. This discrepancy results in surprise:

surprise=-log q(x).

If an event x happens:

({{q(x), -log q(x), Level of suprise}, {1, 0bit, Not at all.}, {1/2, 1bit, As if flipping a coins, landing heads up.}, {1/1024, 10bit, As if flipping  10 coins, all landing heads up!}})

Relative surprise: Now, I observe the world with my friend. We hold different beliefs q(x) (me) and q'(x) (my friend). Thus, even if we witness the same event x, we might experience different levels of surprise. My surprise relative to my friend:

relative surprise=(-log q(x))-(-log q'(x))
=log ((q'(x))/(q(x))).

Compare belief with reality: It turns out that my friend is  omniscient, that their belief equals the reality, i.e. q'(x)=p(x), then relative to them, I would feel more surprised.

This is because

when I am more surprised than my friend, it indicates that I underestimated the likelihood of an event occurring,

when I am less surprised, it means I overestimated it.

Consequently, I tend to be more surprised than my friend more frequently, leading to a higher average relative surprise for me.

(Subscript[D, KL])(p||q)=(∑_(x∈𝒳))p(x)log ((p(x))/(q(x))).

The KL divergence (Subscript[D, KL]) precisely measures this relative, comparative level of surprise, quantifying the deviation of my belief q(x) from reality p(x).

({{x, Biden, Trump}, {(Subscript[q, Pro-Biden])(x), 0.999, 0.001}, {(Subscript[q, Neutral])(x), 0.5, 0.5}, {(Subscript[q, Pro-Trump])(x), 0.001, 0.999}})

({{(Subscript[D, KL])(p||q), , belief q(x), , }, {[bit], , (Subscript[q, Pro-Biden]), (Subscript[q, Neutral]), (Subscript[q, Pro-Trump])}, {reality p(x), (Subscript[p, Biden wins]), 0, 0.99, 9.94}, {, (Subscript[p, Election ties]), 3.98, 0, 3.98}, {, (Subscript[p, Trump wins]), 9.94, 0.99, 0}})

Block[{ps},
ps={#,1-#}&/@{0.001,0.5,0.999};
TableForm@Table[p1.Log2[p1/p2],{p1,ps},{p2,ps}]]

Code.

Variational Free Energy

In practice, it does not matter which distribution is treated as belief and which is reality. Both objectives are fine:

(Subscript[D, KL])((Subscript[p, MF])||p)    or    (Subscript[D, KL])(p||(Subscript[p, MF])).

However, (Subscript[D, KL])((Subscript[p, MF])||p) is much easier to evaluate analytically, because sampling s from the uncorrelated (factorized) distribution (Subscript[p, MF]) is tractable.

So we choose the objective function to be

By minimizing the objective function, we expect to bring (Subscript[p, MF])(s) close to p(s) as much as possible.

p(s)=(1/Z)(e^(-β E(s)))=(1/Z)exp((1/T)(J(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j])+h(∑_i)(Subscript[s, i]))),
(Subscript[p, MF])(s)=(∏_i)(Subscript[p, MF])((Subscript[s, i]))=(∏_i)((1+m (Subscript[s, i]))/2),

(Subscript[D, KL])((Subscript[p, MF])||p)=β F+log Z,

F is the variational free energy, given by

V - lattice volume

q - coordination number (assuming square lattice with q=4)

Derive Eq. () by substituting Eq. () into Eq. ().

(Subscript[D, KL])((Subscript[p, MF])||p)=(∑_s)(Subscript[p, MF])(s)log(((Subscript[p, MF])(s))/(p(s)))
=(∑_s)(Subscript[p, MF])(s)(log(Subscript[p, MF])(s)-log p(s))
=(Underscript[(Underscript[(∑_s)(Subscript[p, MF])(s)(-log p(s)), ⏟]), (Ι)])+(Underscript[(Underscript[(∑_s)(Subscript[p, MF])(s)log (Subscript[p, MF])(s), ⏟]), (ΙΙ)])

Evaluate the first term:

(Ι)=(∑_s)(Subscript[p, MF])(s)(-log p(s))
=(∑_s)(Subscript[p, MF])(s)(-log((1/Z)(e^(-β E(s)))))
=(∑_s)(Subscript[p, MF])(s)(-log((1/Z))+β E(s))
=(∑_s)(Subscript[p, MF])(s)(log Z+β (-J(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j])-h(∑_i)(Subscript[s, i])))
=log Z(∑_s)(Subscript[p, MF])(s)+β(∑_s)(Subscript[p, MF])(s)(-J(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j])-h(∑_i)(Subscript[s, i])).

By Eq. () and Eq. (), we can derive

(∑_s)(Subscript[p, MF])(s)=1,
(∑_s)(Subscript[p, MF])(s) (Subscript[s, i])=m,
(∑_s)(Subscript[p, MF])(s) (Subscript[s, i])(Subscript[s, j])=({{{, ({{(m^2), i≠j,}, {1, i=j.}})}})

Their detailed derivations are as follows

Normalization

(∑_s)(Subscript[p, MF])(s)=(∑_s)(∏_i)(Subscript[p, MF])((Subscript[s, i]))
=(∏_i)(∑_(Subscript[s, i]))(Subscript[p, MF])((Subscript[s, i]))
=(∏_i)1=1.

1st moment

(∑_s)(Subscript[p, MF])(s) (Subscript[s, i])=(∑_s)(∏_j)(Subscript[p, MF])((Subscript[s, j]))(Subscript[s, i])
=((∑_(Subscript[s, i]))(Subscript[p, MF])((Subscript[s, i]))(Subscript[s, i]))(∏_(j≠i))((∑_(Subscript[s, j]))(Subscript[p, MF])((Subscript[s, j])))
=m(∏_(j≠i))1=m.

2nd moment (assuming i≠j)

(∑_s)(Subscript[p, MF])(s) (Subscript[s, i])(Subscript[s, j])=(∑_s)(∏_k)(Subscript[p, MF])((Subscript[s, k]))(Subscript[s, i])(Subscript[s, j])
=((∑_(Subscript[s, i]))(Subscript[p, MF])((Subscript[s, i]))(Subscript[s, i]))((∑_(Subscript[s, j]))(Subscript[p, MF])((Subscript[s, j]))(Subscript[s, j]))(∏_(k≠i,j))((∑_(Subscript[s, k]))(Subscript[p, MF])((Subscript[s, k])))
=(m^2)(∏_(k≠i,j))1=(m^2).

Therefore,

(Ι)=log Z+β (-J(∑_(⟨i j⟩))(m^2)-h(∑_i)m)
=log Z+β (-J(q V/2)(m^2)-h V m)
=log Z+β (-((q J)/2) (m^2)-h m)V.

If we define

E=(-((q J)/2) (m^2)-h m)V,

the above equation can be written as

(Ι)=log Z+β E.

Evaluate the second term:

(ΙΙ)=(∑_s)(Subscript[p, MF])(s)log (Subscript[p, MF])(s)
=(∑_s)(∏_j)(Subscript[p, MF])((Subscript[s, j]))log((∏_i)(Subscript[p, MF])((Subscript[s, i])))
=(∑_s)(∏_j)(Subscript[p, MF])((Subscript[s, j]))((∑_i)log (Subscript[p, MF])((Subscript[s, i])))
=(∑_i)(∑_s)(∏_j)(Subscript[p, MF])((Subscript[s, j]))log (Subscript[p, MF])((Subscript[s, i]))
=(∑_i)((∑_(Subscript[s, i]))(Subscript[p, MF])((Subscript[s, i]))log (Subscript[p, MF])((Subscript[s, i])))(∏_(j≠i))((∑_(Subscript[s, j]))(Subscript[p, MF])((Subscript[s, j])))
=(∑_i)((∑_((Subscript[s, i])=±1))((1+m (Subscript[s, i]))/2)log ((1+m (Subscript[s, i]))/2))(∏_(j≠i))1
=(∑_i)(((1+m)/2)log ((1+m)/2)+((1-m)/2)log ((1-m)/2))
=(((1+m)/2)log ((1+m)/2)+((1-m)/2)log ((1-m)/2))V.

If we define

S=-(((1+m)/2)log ((1+m)/2)+((1-m)/2)log ((1-m)/2))V,

the second term is simply (ΙΙ)=-S.

(Subscript[D, KL])((Subscript[p, MF])||p)=(Ι)+(ΙΙ)
=log Z+β E-S
=log Z+β (E-T S)
=log Z+β F,

Note that (Subscript[F, true]):=-T log Z is the true free energy of the system, Eq. () can be written as

Since (Subscript[F, ture]) (or Z) is independent of the model parameter m, it can be dropped from the objective function without affecting the optimization of m.

=> To minimize (Subscript[D, KL])((Subscript[p, MF])||p), we only need to minimize F.

Given (Subscript[D, KL])((Subscript[p, MF])||p)≥0,

i.e. the variational free energy F serves as an upper bound of the true free energy.

Plot[{1,2,3Exp[-t]+2},{t,0,5},PlotStyle->{Black,Directive[Red,Dashed],Red},Filling->{3->{1,Lighter[Red,0.8]}},FrameTicks->None,FrameLabel->{"Optimization steps",None},Epilog->{Text["\!\(-T log Z\)",{1,1},{0,1}],{Red,Text["\!\( F\)",{1,3.2},{0,-1}],Text["\!\(F\_★\)",{5,2},{1,-1}]},{{Arrowheads[{-0.05,0.05}],Arrow[{{#,1},{#,3 Exp[-#]+2}}&@1]},{{Lighter[Red,0.8],Table[Translate[#,{0.618,1}r],{r,0.1CirclePoints[32]}]},#}&@Text["\!\(T D\_KL\)",{1,2},{-1.1,0}]}},PlotRange->{0,5},ImageSize->180]

DynamicModule[{J=1,f},
f[m_,T_,h_]:=(-2J m^2-h m)-T(-(1+m)/2 Log[(1+m)/2]-(1-m)/2 Log[(1-m)/2]);
Manipulate[Plot[f[m,T,h],{m,-1,1},PlotLabel->StringTemplate["\!\((T,h)/J=(`1`,`2`)\)"][NumberForm[T,{3,2}],NumberForm[h,{3,2}]],PlotStyle->Black,Axes->False,FrameLabel->{"\!\( m\)","\!\(F/J V\)"},ImagePadding->{{55,5}, {40, 5}},ImageSize->200],{{T,3.,"\!\(T/J\)"},0.01,8,ImageSize->Small},{{h,0.,"\!\(h/J\)"},-2,2,ImageSize->Small},LabelStyle->"InlineFormula",Paneled->False]]

Parameter Optimization

-h-q J m+T arctanh m=0.

Verify Eq. ().

D[(-q J/2 m^2-h m)-T(-(1+m)/2 Log[(1+m)/2]-(1-m)/2 Log[(1-m)/2]),m]//FullSimplify

Manipulate[Block[{q=4,J=1,ms},
ms=Sort@DeleteDuplicates[(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,#,100]},AccuracyGoal->3,PrecisionGoal->3])&/@{-1.,1.},Abs[#1-#2]<0.01&];
Plot[{m,Tanh[(q J m+h)/T]},{m,-1,1},Epilog->{{AbsoluteThickness[1],Opacity[0.3],Circle[#{1,1},Offset@{4,4}]&/@ms}},PlotLabel->StringTemplate["\!\(m∈`1`\)"][NumberForm[Chop[ms,0.01],{3,2}]],FrameLabel->{"\!\( m\)",None}]],{{T,3.,"\!\(T/J\)"},0.01,8,ImageSize->Small},{{h,0.,"\!\(h/J\)"},-2,2,ImageSize->Small},LabelStyle->"InlineFormula",Paneled->False]

DynamicModule[{q=4,J=1,m,mval,legend,bg,T0=4.,h0=0.,mode=0},
mval[T_,h_,m0_:None]:=If[T==0.,Sign[h],(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,If[m0===None,N@Sign[h],m0],100]},AccuracyGoal->3,PrecisionGoal->3])];
legend=Graphics[Raster[List/@Range[0,1,0.05],{{0,-1.05},{1,1.05}},ColorFunction->ColorData["BlueRedSplit"]],PlotLabel->Style["\!\( m\_★\)",Black,"InlineFormula"],PlotRangePadding->0,Frame->True,FrameTicks->{{None,None},{None,None}},BaselinePosition->Scaled[0.32],AspectRatio->10,ImagePadding->{{2,2},{2,0}},ImageSize->{Automatic,112}];
bg=First@ListContourPlot[Flatten[Table[{T,h,mval[T,h]},{T,Sort@Join[Range[0,8,0.2],Range[3.81,4.19,0.01]]},{h,Sort@Join[Range[-2,2,0.2],Range[-0.1,0.1,0.01]]}],1],Contours->21,ContourStyle->None];
Row[{LocatorPane[Dynamic[{T0,h0}],Graphics[{bg,{{Black,Line[{#,0}&/@{0,q J}]},{Black,Point[{q J,0}]}},Dynamic@{Opacity[0.2],Switch[mode,0,Line[{#,h0}&/@{0,8}],1,Line[{T0,#}&/@{-2,2}]]}},Frame->True,FrameLabel->{"\!\(T/J\)","\!\(h/J\)"},AspectRatio->1,ImageSize->180],{{0,-2},{8,2}},Appearance->None,BaselinePosition->Scaled[0.55]],legend,
Column[{Toggler[Dynamic[mode],{0->Dynamic[StringTemplate["\!\( h/J=`1`\)"][NumberForm[h0,{3,2}]]],1->Dynamic[StringTemplate["\!\( T/J=`1`\)"][NumberForm[T0,{3,2}]]]},BaseStyle->"InlineFormula"],
Dynamic@Switch[mode,0,Plot[If[Abs[h0]<0.01,{mval[T,0.,1.],mval[T,0.,-1.]},mval[T,h0]],{T,0,8},PlotPoints->21,MaxRecursion->3,FrameLabel->{"\!\(T/J\)",None},PlotStyle->Black,PlotRange->{{0,8},{-1,1}},PlotRangePadding->{0.2,0.1},ImageSize->180],1,Plot[mval[T0,h],{h,-2,2},Exclusions->If[T0>=q J,{},{0}],PlotPoints->21,MaxRecursion->3,FrameLabel->{"\!\(h/J\)",None},PlotStyle->Black,PlotRange->{{-2,2},{-1,1}},PlotRangePadding->{0.1,0.1},ImageSize->180]]},Alignment->Right]}]]

h=0: there is a phase transition, with the critical temperature

T≥(Subscript[T, c]): Only one stable solution m=0. => Ising spins are disordered, resulting in zero magnetization on average.

T<(Subscript[T, c]): m=0 is unstable, two stable solutions emerge at m≠0 (as T→0, m→±1). => Ising spins are ordered (spontaneously align in the same direction), resulting in finite magnetization.

m is also called the order parameter — a local observable that indicates the ordering of Ising spins.

h≠0: there is no sharp transition, only a crossover.

T≫q J: m∼h/T varies smoothly.

T≲q J: The anti-aligned solution (h m<0) will become meta-stable or unstable. The aligned solution (h m>0) is always stable.

Magnetization curve shows how magnetization (Subscript[m, ★]) responding to external magnetic field h.

T>(Subscript[T, c]): Paramagnetic (spin disordered) phase. Magnetization curve is continuous with finite slope near h→0.

T<(Subscript[T, c]): Ferromagnetic (spin ordered) phase. Magnetization curve is discontinuous across h→0.

T=(Subscript[T, c]): Critical point of the magnetic transition (order-disorder transition). Magnetization curve is continuous but with divergent slope near h→0.

Model Prediction

DynamicModule[{q=4,J=1,mval,f,legend,bg,T0=4.,h0=0.,mode=0},
mval[T_,h_,m0_:None]:=Block[{m},If[T==0.,Sign[h],(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,If[m0===None,N@Sign[h],m0],100]},AccuracyGoal->3,PrecisionGoal->3])]];
f[T_,h_]:=Block[{m},
m=If[h==0.,If[T==0.,1.,mval[T,0.,1.]],mval[T,h]];
(-2J m^2-h m)-T If[Abs[m]==1.,0,(-(1+m)/2 Log[(1+m)/2]-(1-m)/2 Log[(1-m)/2])]];
legend=Graphics[Raster[List/@Range[0,1,0.05],{{0,-1.05},{1,1.05}},ColorFunction->ColorData["DeepSeaColors"]],PlotLabel->Style["   \!\( F\_★/J\)",Black,"InlineFormula"],PlotRangePadding->0,Frame->True,FrameTicks->{{None,None},{None,None}},BaselinePosition->Scaled[0.3],AspectRatio->10,ImagePadding->{{2,5},{2,0}},ImageSize->{Automatic,112}];
bg=First@ListContourPlot[Flatten[Table[{T,h,f[T,h]},{T,0,8,0.2},{h,-2,2,0.1}],1],Contours->Range[-6,-2,0.2],ContourStyle->None,ColorFunction->"DeepSeaColors"];
Row[{LocatorPane[Dynamic[{T0,h0}],Graphics[{bg,{{Black,Line[{#,0}&/@{0,q J}]},{Black,Point[{q J,0}]}},Dynamic@{Opacity[0.2],Switch[mode,0,Line[{#,h0}&/@{0,8}],1,Line[{T0,#}&/@{-2,2}]]}},Frame->True,FrameLabel->{"\!\(T/J\)","\!\(h/J\)"},AspectRatio->1,ImageSize->180],{{0,-2},{8,2}},Appearance->None,BaselinePosition->Scaled[0.55]],legend,
Column[{Toggler[Dynamic[mode],{0->Dynamic[StringTemplate["\!\( h/J=`1`\)"][NumberForm[h0,{3,2}]]],1->Dynamic[StringTemplate["\!\( T/J=`1`\)"][NumberForm[T0,{3,2}]]]},BaseStyle->"InlineFormula"],
Dynamic@Switch[mode,0,Plot[f[T,h0],{T,0,8},Exclusions->None,PlotPoints->21,MaxRecursion->3,FrameLabel->{"\!\(T/J\)",None},PlotStyle->Black,PlotRange->{{0,8},{-6,-2}},PlotRangePadding->{0.2,0.2},ImageSize->170],1,Plot[f[T0,h],{h,-2,2},Exclusions->None,PlotPoints->21,MaxRecursion->3,FrameLabel->{"\!\(h/J\)",None},PlotStyle->Black,PlotRange->{{-2,2},{-6,-2}},PlotRangePadding->{0.1,0.2},ImageSize->170]]},Alignment->Right]}]
]

d(Subscript[F, ★])=-(Subscript[P, ★])dV-(Subscript[S, ★])dT-(Subscript[M, ★])dh,

Entropy

(Subscript[S, ★])=-(Subscript[(((∂(Subscript[F, ★]))/(∂T))), V,h]).

Block[{q=4,J=1,mval,f,fT},
mval[T_,h_,m0_:None]:=Block[{m},If[T==0.,Sign[h],(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,If[m0===None,N@Sign[h],m0],100]},AccuracyGoal->3,PrecisionGoal->3])]];
f[T_,h_]:=Block[{m},
m=If[h==0.,If[T==0.,1.,mval[T,0.,1.]],mval[T,h]];
(-2J m^2-h m)-T If[Abs[m]==1.,0,(-(1+m)/2 Log[(1+m)/2]-(1-m)/2 Log[(1-m)/2])]];
fT=Interpolation[Table[{T,f[T,0.]},{T,0,8,0.02}]];
Plot[-fT'[T],{T,0,8},PlotStyle->Black,GridLines->{{q J},None},FrameLabel->{"\!\(T/J\)","\!\( S\)"},Epilog->{Text["\!\(h=0\)",{8,0},{1,-0.5}]}]]

Specific heat

(Subscript[c, V])=(T/V)(Subscript[(((∂(Subscript[S, ★]))/(∂T))), V,h])=-(T/V)(Subscript[((((∂^2)(Subscript[F, ★]))/(∂(T^2)))), V,h]).

Block[{q=4,J=1,mval,f,fT},
mval[T_,h_,m0_:None]:=Block[{m},If[T==0.,Sign[h],(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,If[m0===None,N@Sign[h],m0],100]},AccuracyGoal->3,PrecisionGoal->3])]];
f[T_,h_]:=Block[{m},
m=If[h==0.,If[T==0.,1.,mval[T,0.,1.]],mval[T,h]];
(-2J m^2-h m)-T If[Abs[m]==1.,0,(-(1+m)/2 Log[(1+m)/2]-(1-m)/2 Log[(1-m)/2])]];
fT=Interpolation[Table[{T,f[T,0.]},{T,0,8,0.02}]];
Plot[-T fT''[T],{T,0,8},PlotStyle->Black,GridLines->{{q J},None},FrameLabel->{"\!\(T/J\)","\!\(c\_V\)"},Epilog->{Text["\!\(h=0\)",{8,1.5},{1,1}]}]]

Mean field theory: (Subscript[c, V]) becomes discontinuous at T=(Subscript[T, c]).

Beyond mean field: (Subscript[c, V]) diverges at T=(Subscript[T, c]) as

with a critical exponent α.

Magnetization

(Subscript[M, ★])=-(Subscript[(((∂(Subscript[F, ★]))/(∂h))), V,T]).

Magnetization density (Subscript[m, ★])=(Subscript[M, ★])/V.

Block[{q=4,J=1,mval},
mval[T_,h_,m0_:None]:=Block[{m},If[T==0.,If[m0===None,Sign[h],Sign[m0]],(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,If[m0===None,N@Sign[h],m0],100]},AccuracyGoal->3,PrecisionGoal->3])]];
ListLinePlot[Table[{T,Chop@mval[T,0,1.]},{T,0,8,0.02}],GridLines->{{q J},None},PlotStyle->Black,PlotRange->{{0,8},{0,1}},PlotRangePadding->{0.2,0.05},FrameLabel->{"\!\(T/J\)","\!\(m\_★\)"},Epilog->{Text["\!\(h→0\_+\)",{8,1},{1,1}]}]]

Mean field theory: (Subscript[m, ★]) drops to zero following (Subscript[m, ★])∼((((Subscript[T, c])-T))^(1/2)) as T→(Subscript[T, c]) (from below).

Beyond mean field: (Subscript[m, ★]) drops to zero as T→(Subscript[T, c]) following

with a critical exponent β.

Magnetic susceptibility: in the paramagnetic phase, (Subscript[m, ★]) responses to h linearly (for small h), a behavior known as linear response. The linear coefficient is defined as susceptibility

χ=((∂(Subscript[m, ★]))/(∂h))=-(1/V)(Subscript[((((∂^2)(Subscript[F, ★]))/(∂(h^2)))), V,T]).

Block[{q=4,J=1,mval},
mval[T_,h_,m0_:None]:=Block[{m},If[T==0.,Sign[h],(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,If[m0===None,N@Sign[h],m0],100]},AccuracyGoal->3,PrecisionGoal->3])]];
Block[{δh=0.00001},ListLinePlot[Table[{T,(mval[T,δh]-mval[T,-δh])/(2δh)},{T,q J,8,0.02}],GridLines->{{q J},None},PlotStyle->Black,PlotRange->{{0,8},{0,5}},PlotRangePadding->{0.2,0.2},FrameLabel->{"\!\(T/J\)","\!\(χ\)"},Epilog->{Text["\!\(h=0\)",{8,5},{1,1}],{LightGray,Rectangle[{0,0},{q J,5}]},Text["\!\(χ→∞\)",{2,5},{0,1}]}]]]

Mean field theory: χ diverges as χ∼1/(T-(Subscript[T, c])) as T→(Subscript[T, c]) (from above).

Beyond mean field: χ diverges as T→(Subscript[T, c]) as

with a critical exponent γ.

At the critical point T=(Subscript[T, c]), χ→∞ at h=0, indicating that (Subscript[m, ★]) is no longer linear in h.

Block[{q=4,J=1,mval},
mval[T_,h_,m0_:None]:=Block[{m},If[T==0.,If[m0===None,Sign[h],Sign[m0]],(m/.FindRoot[Tanh[(q J m+h)/T]==m,{m,FixedPoint[Tanh[(q J #+h)/T]&,If[m0===None,N@Sign[h],m0],100]},AccuracyGoal->3,PrecisionGoal->3])]];
ListLinePlot[Table[{h,mval[q J,h]},{h,2Range[-1,1,0.1]^3}],GridLines->{{q J},None},PlotStyle->Black,PlotRange->{{-2,2},{-1,1}},PlotRangePadding->{0.2,0.05},FrameLabel->{"\!\(h/J\)","\!\(m\_★\)"},Epilog->{Text["\!\(T=T\_c\)",{2,-1},{1,-1}]}]]

Mean field theory: (Subscript[m, ★]) responses to h as (Subscript[m, ★])∼(h^(1/3)) at T=(Subscript[T, c]).

Beyond mean field: (Subscript[m, ★]) responses to h as

with a critical exponent δ.

Landau Theory of Phase Transition

Spontaneous Symmetry Breaking

Grid[{{"Name","\!\(ℤ\_2≅O(1)\)","\!\(O(2)\)","\!\(O(3)\)"},{"Meaning","2\[Hyphen]fold reflection","Circular rotation","Spherical rotation"},Prepend[Graphics[#,PlotRange->1,ImageSize->{80,80}]&/@{{{AbsolutePointSize[5],Red,Point[{0,#}&/@{-0.7,0.7}]},{AbsoluteThickness[1],Dashed,Line[{0,#}&/@{-0.7,0.7}]},{AbsoluteThickness[1.2],Line[{#,0}&/@{-1,1}]}},{{AbsoluteThickness[1.2],Circle[{0,0},0.9]},{Arrowheads[0.15],Red,Arrow[ReIm@{0,0.9Exp[I π/3]}]},Point[{0,0}]},{{AbsoluteThickness[1.2],Circle[{0,0},0.9],AbsoluteThickness[1],Circle[{0,0},{0.9,0.2},{-π,0}],{Dashed,Circle[{0,0},{0.9,0.2},{0,π}]},Circle[{0,0},{0.7,0.9},{-π/2,π/2}]},{Arrowheads[0.15],Red,Arrow[{0.7,0.9}#&/@ReIm@{0,Exp[I π/4]}]},Point[{0,0}]}},""],
{Style["Order\nparameter",LineSpacing->{0,12}],"\!\( m\)",Style["\!\(m=(({{(Subscript[m, 1])}, {(Subscript[m, 2])}}))\)",SpanMaxSize->∞],Style["\!\(m=(({{(Subscript[m, 1])}, {(Subscript[m, 2])}, {(Subscript[m, 3])}}))\)",SpanMaxSize->∞]},{"Transform","\!\(m→±m\)","\!\(m→E\^\(I θ σ\^y\)m\)","\!\(m→E\^\(I\/2 θ·σ\)m\)"}},Dividers->{2->True},ItemStyle->{Automatic,{"Text","Text",{"InlineFormula"}}}]

For example, if the microscopic energy model E(s)=E(-s) is invariant under the Ising symmetry s→-s, the probability distribution p(s)∝(e^(-β E(s))) will inherit the same symmetry, hence any macroscopic properties ⟨X⟩ must also be symmetric,

∀X:⟨X(s)⟩
=(∑_s)X(s)p(s)(Overscript[=, s->-s])(∑_s)X(-s)p(-s)
=(∑_s)X(-s)p(s)=⟨X(-s)⟩.

This would imply the magnetization density m:=⟨m(s)⟩ to vanish, because it is odd under the Ising symmetry,

m(s)=(1/V)(∑_i)(Subscript[s, i])=-m(-s),

thus ⟨m(s)⟩=-⟨m(-s)⟩=-⟨m(s)⟩ => m=0.

Block[{lt=0.8},
Graphics[{{Lighter[Red,lt],Rectangle[{-1,-1},{1,1}]},{AbsoluteThickness[10],Lighter[Purple,lt],Line[{#,0}&/@{-1,0}],White,Line[{#,0}&/@{0,1}]},Text["Symm.",{0.5,0}],Text["(\!\(m≠0\))",{-0.5,0},{0,1}],Text["(\!\(m"<>#1<>"0\))",{0,#2}]&@@@{{">",0.45},{"<",-0.45}},{Red,Text["ESB",{0,#}]&/@{-0.7,0.7}},{Purple,Line[{#,0}&/@{-1,0}],Text["SSB",{-0.5,0},{0,-0.5}]},{Point[{0,0}]}},PlotRange->1,PlotRangeClipping->True,Frame->True,FrameTicks->{{{{0,0}},None},{{{0,0}},None}},FrameLabel->{"\!\(T-T\_c\)","\!\( h\)"},ImageSize->160]]

Explicit symmetry breaking (ESB): When h≠0, the external magnetic field explicitly breaks the Ising symmetry, s.t. E(s) is no longer invariant under s→-s. If the cause doesn’t have the symmetry, so doesn’t its effect.

Spontaneous symmetry breaking (SSB): When h=0 and T<(Subscript[T, c]), while the underlying microscopic model E(s) has the Ising symmetry, the system settles into a macroscopic state that does not exhibit this symmetry.

SSB presents a paradox where a symmetric cause leads to asymmetric effects. There are several interpretations:

Collective choice: At low temperature, the Ising spins collectively choose a particular symmetry breaking order to lower the free energy. SSB is a collective effect driven by interaction.

Historical contingency: At the critical point, the Ising spins fluctuates strongly. Once a small, random fluctuation biases towards a symmetry breaking pattern, the bias can be amplified by interaction and gets frozen in the system at low temperature. SSB is a historical contingency driven by critical fluctuations.

Ergodicity breaking: Ergodicity — the assumption that a system can sample all feasible microstate in its equilibrium distribution — can break down when the distribution splits into disconnected modes (peaks). Once settled in one mode, the system can hardly explore other modes via local fluctuations. The symmetry appears to be broken despite that the distribution is still symmetric.

Manipulate[Block[{a0=1,b0=1/2,Tc=1,f,Z,p},
f[m_]:=a0(T-Tc) m^2+b0 m^4;
Off[General::munfl];
Z=NIntegrate[Exp[-f[m]/T],{m,-2,2}];
p[m_]:=Exp[-f[m]/T]/Z;
Plot[p[m],{m,-2,2},PlotStyle->Black,Filling->Axis,PlotRange->{0,All},PlotRangePadding->{0.1,Scaled[.05]},PlotLabel->StringTemplate["\!\(T/T\_c=`1`\)"][NumberForm[T,{3,2}]],FrameLabel->{"\!\( m\)","\!\(p(m)\)"},ImagePadding->{{40, 2}, {40,2}},ImageSize->180]],{{T,1.,"\!\(T/T\_c\)"},0.01,2,ImageSize->Small},LabelStyle->"InlineFormula",Paneled->False]

SSB is an asymmetric illusion in a symmetric probability distribution due to mode collapse.

Classification of Phases

Grid[{{"Phase","Order","Broken symmetry"},{"Ferromagnet","Spin alignment","Spin rotational symmetry"},{"Liquid crystal\n(nematic)","Molecule alignment","Spatial rotation symmetry"},{"Crystalline solid","Atom arrangement\n(periodic lattice)","Translation and rotation symmetries"},{"Superfluid","Phase coherence\nof particles","\!\(U(1)\) symmetry\n(particle number conservation)"},{"Electricmagnetic field","Phase coherence\nof strings","Magnetic 1\[Hyphen]form symmetry"}},Dividers->{{},{2->True}},ItemStyle->{Directive["Text",TextAlignment->Center],{1->Directive[14,Bold,FontFamily->"CMU Concrete"]}},Background->{None,{{White,LightGray}}}]

Landau Free Energy

Given a symmetry G of the system, an order parameter is a local observable that transforms non-trivially under the symmetry action.

For example, G=(Subscript[ℤ, 2]) in Ising model is the Ising (spin-flip) symmetry. Given m=⟨(Subscript[s, i])⟩,

(Subscript[ℤ, 2]):m→-m,

thus m qualifies as an order parameter, as it can indicate symmetry breaking:

m({{{, ({{=0, (Subscript[ℤ, 2]) symmetric,}, {≠0, (Subscript[ℤ, 2]) symmetry breaking.}})}})

Landau theory employs a Taylor series expansion of the free energy (density) in terms of the order parameter, without considering the microscopic details of the system.

For example, expanding the variational free energy F(V,T,h;m) in Eq. () with respect to the order parameter m gives

f:=(F/V)=-T log 2-h m+((T-(Subscript[T, c]))/2)(m^2)+((T (m^4))/12)+….

Verify Eq. ().

Block[{f},
f[m_,T_,h_]:=(-q J/2 m^2-h m)-T(-(1+m)/2 Log[(1+m)/2]-(1-m)/2 Log[(1-m)/2]);
Series[f[m,T,h],{m,0,4}]]

Landau argues that, in the absence of external field (i.e. h=0), the free energy density f should take the following general form

Manipulate[Plot[a m^2+b m^4,{m,-1,1},PlotRange->{{-1,1},{-0.5,1}},PlotLabel->StringTemplate["\!\(a=`1`,b=`2`\)"][NumberForm[a,{3,2}],NumberForm[b,{3,2}]],PlotStyle->Black,FrameLabel->{"\!\( m\)","\!\( f(T;m)\)"},ImageSize->200],{{a,-1.,"\!\( a\)"},-1,1,ImageSize->Small},{{b,2.,"\!\( b\)"},0.5,2,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Symmetry argument: With h=0, the system is (Subscript[ℤ, 2]) symmetric => f(T;m) must be an even function of m as

f(T;m)=f(T;-m),

So the expansion Eq. () must only contain even powers of m.

To include an external magnetic field h, simply add -h m to the free energy density:

Symmetry argument: h≠0 explicitly breaks the (Subscript[ℤ, 2]) symmetry, such that f(T,h;m) is no longer an even function of m, therefore a linear term -h m is allowed.

Phase Diagram

The (m^0) term (Subscript[f, 0])(T) is unimportant (unrelated to m).

The (m^2) coefficient a(T), called the tuning parameter, tunes the phase transition:

Phases along h=0:

Grid[{{Style["Tuning\nparameter",LineSpacing->{0,12}],Style["Order\nparameter",LineSpacing->{0,12}],"Phase","\!\(ℤ\_2\) symmetry"},{"\!\(a(T)>0\)","\!\(m\_★=0\)","Disorder","Preserved"},
{"\!\(a(T)<0\)","\!\(m\_★≠0\)","Ordered","Broken (spontaneously)"}},ItemStyle->"InlineFormula",Dividers->{None,{2->True}},Alignment->Left]

The (m^4) coefficient b(T) should be positive to ensure the stability (if not, (m^6) term must be considered). We may assume

DynamicModule[{a=-0.6,h=0.4,b=1,pd,legend},
legend=Graphics[Raster[List/@Range[0,1,0.05],{{0,-1.05},{1,1.05}},ColorFunction->ColorData["BlueRedSplit"]],PlotLabel->Style["\!\( m\_★\)",Black,"InlineFormula"],PlotRangePadding->0,Frame->True,FrameTicks->{{None,None},{None,None}},BaselinePosition->Scaled[0.32],AspectRatio->10,ImagePadding->{{0,2},{2,0}},ImageSize->{Automatic,112}];
pd=ContourPlot[-Sign[h]Root[2a #+4b #^3+Abs[h]&,1],{a,-1,1},{h,-0.6,0.6},Contours->Range[-1,1,0.05],ContourStyle->None,Exclusions->None,Epilog->{{Opacity[0.2],Table[Polygon@Append[Table[{a,s(-2a/3)^(3/2)/b^(1/2)},{a,-1,0,0.1}],{-1,0}],{s,{1,-1}}]},Line[{#,0}&/@{-1,0}],Point[{0,0}]},FrameLabel->{"\!\(a/b\)","\!\(h/b\)"},ImageSize->180];
Row[{Row[{LocatorPane[Dynamic[{a,h}],pd,{{-1,-0.6},{1,0.6}},BaselinePosition->Scaled[0.55]],legend}],Dynamic[Plot[a m^2+b m^4-h m,{m,-1,1},PlotRange->{{-1,1},{-0.8,1}},PlotStyle->Black,Epilog->{Translate[{{Point[{0,0}]},{Opacity[0.2],AbsoluteThickness[1],InfiniteLine[{{0,0},{0,1}}]},Text["\!\(m\_★\)",{0,0},{0,-1}]},{#,a #^2+b #^4-h #}]&[-Sign[h]Root[2a #+4b #^3+Abs[h]&,1]]},PlotLabel->StringTemplate["\!\((a,h)/b=(`1`,`2`)\)"][
NumberForm[a,{3,2}],NumberForm[h,{3,2}]],FrameLabel->{"\!\( m\)","\!\( f(m)\)"},AspectRatio->1,ImageSize->{Automatic,190}]]},Spacer[1]]]

Critical Behaviors

the stable solution of order parameter

which exhibits (Subscript[m, ★])∼((((Subscript[T, c])-T))^(1/2)) behavior near the critical point;

Block[{a0=1,b0=1,Tc=1},Plot[Piecewise[{{Sqrt[a0/(2b0)]Sqrt[Tc-T],T<Tc}},0],{T,0.9Tc,1.1Tc},PlotStyle->Black,Exclusions->None,FrameLabel->{"\!\(T/T\_c\)","\!\(m\_★\)"}]]

the optimal free energy density

Derive Eq. () and Eq. ().

f(m)=(Subscript[f, 0])+a (m^2)+b (m^4).

((∂f)/(∂m))=2a m+4 b (m^3)=2m (a+2b (m^2))=0,

m=0,
m=±(Sqrt[((-a)/(2b))])   (for a<0).

When a≥0, m=0 is the stable solution.

When a<0, m=±(Sqrt[((-a)/(2b))]) are stable solutions.

(Subscript[m, ★])=({{{, ({{±(Sqrt[((-a)/(2b))]), a<0,}, {0, a>=0.}})}})

(Subscript[f, ★]):=f((Subscript[m, ★]))
=(Subscript[f, 0])+a (Subscript[m, ★]^2)+b (Subscript[m, ★]^4).

When a>=0, (Subscript[m, ★])=0, so (Subscript[f, ★])=(Subscript[f, 0]).

When a<0,

(Subscript[f, ★])=(Subscript[f, 0])+a (((-a)/(2b)))+b (((((-a)/(2b))))^2)
=(Subscript[f, 0])-((a^2)/(2b))+((a^2)/(4b))
=(Subscript[f, 0])-((a^2)/(4b)).

(Subscript[f, ★])=({{{, ({{(Subscript[f, 0])-((a^2)/(4b)), a<0,}, {(Subscript[f, 0]), a>=0.}})}})

(Subscript[c, V])=-T(((∂^2)(Subscript[f, ★]))/(∂(T^2)))=({{{, ({{((Subscript[a, 0]^2)/(2(Subscript[b, 0])))(Subscript[T, c]), T<(Subscript[T, c]),}, {0, T>=(Subscript[T, c]),}})}})

Block[{a0=1,b0=1,Tc=1},Plot[Piecewise[{{a0^2/(2b0)Tc,T<Tc}},0],{T,0.9Tc,1.1Tc},PlotStyle->Black,Exclusions->None,FrameLabel->{"\!\(T/T\_c\)","\!\(c\_V\)"}]]

Derive Eq. ().

f=(Subscript[f, 0])-h m+a (m^2)+b (m^4).

((∂f)/(∂m))=-h+2 a m+4 b (m^3)=0,

h=2 a m+4 b (m^3).

Manipulate[Block[{a0=1,b0=1,Tc=1,h,ma,mc},
h[m_]:=2a0(T-Tc)m+4b0 m^3;
ma=Sqrt[a0(Tc-T)/(2b0)];mc=Sqrt[a0(Tc-T)/(6b0)];
Show[ParametricPlot[{h[m],m},{m,#2,#3},PlotStyle->#1]&@@@If[T<Tc,{{LightGray,-mc,mc},{Gray,-mc,-ma},{Gray,ma,mc},{Black,-0.5,-ma},{Black,ma,0.5}},{{Black,-0.5,0.5}}],PlotRange->{{-0.2,0.2},{-0.5,0.5}},Epilog->{If[T<Tc,Block[{m,mb},mb=First@DeleteCases[m/.NSolve[h[mc]==h[m]],mc];
{Arrowheads[{{0.01,0.5,Graphics[{AbsoluteThickness[1],Dashing[{}],Line[Offset/@{{-4,3},{0,0},{-4,-3}}]}]}}],Gray,Dashed,Arrow[{h[mc],#}&/@{mc,mb}],Arrow[{h[-mc],#}&/@{-mc,-mb}]}],{}]},Axes->True,AxesOrigin->{0,0},FrameLabel->{"\!\( h\)","\!\( m\)"},AspectRatio->0.618,ImageSize->180]],{{T,1,"\!\( T\)"},0.7,1.5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

In the disorder phase (T>(Subscript[T, c])), the magnetic susceptibility at h=0 is given by

χ=(Subscript[(((∂(Subscript[m, ★]))/(∂h))), h=0])=(1/(2 (Subscript[a, 0])(T-(Subscript[T, c])))),

which diverges as χ∼(((T-(Subscript[T, c])))^(-1)) as T→(Subscript[T, c]) from above, known as the Curie-Weiss law.

At the critical point T=(Subscript[T, c]), a(T)=0 vanishes, then Eq. () reduces to

h=4 (Subscript[b, 0])(m^3),

indicating that (Subscript[m, ★]) responses to h as (Subscript[m, ★])∼(h^(1/3)).

({{Exponent, Definition, Condition, Prediction}, {α, (Subscript[c, V])∼((|T-(Subscript[T, c])|)^(-α)),  , 0}, {β, (Subscript[m, ★])∼((((Subscript[T, c])-T))^β), T<(Subscript[T, c]), 1/2}, {γ, χ∼(((T-(Subscript[T, c])))^(-γ)), T>(Subscript[T, c]), 1}, {δ, (Subscript[m, ★])∼(h^(1/δ)), T=(Subscript[T, c]), 3}})

Limitation of Landau Theory

Limitations of Landau theory:

Absence of critical fluctuations: Fluctuations are long-ranged and highly correlated near the critical point. Ignoring them leads to inaccurate predictions of critical exponents.

({{Exponent, Definition, D=2, D=3, D=4, Mean field}, {α, (Subscript[c, V])∼((|T-(Subscript[T, c])|)^(-α)), 0, 0.110, 0, 0}, {β, (Subscript[m, ★])∼((((Subscript[T, c])-T))^β), 1/8, 0.336, 1/2, 1/2}, {γ, χ∼(((T-(Subscript[T, c])))^(-γ)), 7/4, 1.237, 1, 1}, {δ, (Subscript[m, ★])∼(h^(1/δ)), 15, 4.790, 3, 3}})

Ignorance of dimensionality: Dimensionality of the system significantly influences critical phenomena, which is not account for in mean field theories.

Approaches beyond Landau theory:

Field theory: Promote the order parameter m to a fluctuating field m(x), e.g.

F=∫(d^D)x ((1/2)((((∂_μ)m))^2)+a (m^2)+b (m^4)+…),

such that the spatial dimension D can also be incorporated.

Numerical simulations: Use computer to simulate the system at large scales (e.g. Monte Carlo simulation).

Renormalization group: A powerful framework to analyze how physics systems behaves at different length scales, from which critical exponents can be more accurately estimated.

Conformal bootstrap: A non-perturbative approach leveraging conformal symmetries to classify universality classes and solve for critical exponents.

Renormalization Group

Overview

An elegant conceptual framework to understand phases and phase transitions,

A powerful computational approach to identify critical points and estimate critical exponents.

How to identify representative features? - Representation learning.

How to establish the effective theory from scale to scale? - Coarse graining.

What is the usage of tracking the flow of effective theories? - Renormalization group flow.

Representation Learning

Multivariate Probabilities

Positivity

p(x,y,…)≥0.

Normalization

(∑_(x∈𝒳))(∑_(y∈𝒴))… p(x,y,…)=1.

p(x,y): joint probability. The probability of observing X=x and Y=y jointly.

p(x): marginal probability. The probability of observing X=x regardless of Y (before any observation is made on Y).

p(x|y): conditional probability. The probability of observing X=x after observing Y=y.

The conditional probability is simply proportional to the joint probability, and is normalized by the marginal probability for each given condition separately.

Joint and marginal distributions:

({{p(x,y), ,  , y,  , p(x)}, {, , 1, 2, 3, }, { , 1, 0, (1/2), 0, (1/2)}, {x, 2, 0, 0, (1/6), (1/6)}, { , 3, (1/6), 0, (1/6), (1/3)}, {p(y), , (1/6), (1/2), (1/3),  }})

Conditional distributions:

p(x|y)=p(x,y)/p(y)

({{ ,  , p(x|Y=1), p(x|Y=2), p(x|Y=3)}, { , 1, 0, 1, 0}, {x, 2, 0, 0, (1/2)}, { , 3, 1, 0, (1/2)}})

p(y|x)=p(x,y)/p(x)

({{ ,  , y,  }, { , 1, 2, 3}, {p(y|X=1), 0, 1, 0}, {p(y|X=2), 0, 0, 1}, {p(y|X=3), (1/2), 0, (1/2)}})

Action and Probability

p(x)=(1/Z)(e^(-β E(x)))

This concept of “action” parallels that in classical mechanics, where the principle of least action simply states that the most-likely configuration of x is the one that minimizes the action S(x) — a straightforward interpretation of Eq. ().

log Z is often dropped in describing the action, as it is just a constant shift, not affecting how S(x) depends on x, and Z can always be recovered by normalizing p(x).

Interestingly, S(x) can also be interpreted as the entropy reduction, or the information gain, associated with observing the specific state x, because the average entropy of the distribution p(x) is precisely given by

S=(∑_x)p(x)S(x)=-(∑_x)p(x)log p(x).

What a nice coincidence in the use of the symbol “S”!

S(x)(Underscript[(Overscript[⇌, statistical physics]), machine learning])p(x).

Statistical physics starts with the theoretical model of action (or energy function) S(x)=β E(x) and derive the probability distribution p(x) to characterize the equilibrium thermodynamic properties.

Machine learning starts with the empirical data distribution p(x) and infer or learn a underlying action model S(x) that describe or generate the observed data distribution.

Effective Action

p(x,z)=(e^(-S(x,z))),

p(x)=(∑_z)p(x,z)=(∑_z)(e^(-S(x,z))),

S(x)=-log p(x)=-log(∑_z)(e^(-S(x,z))).

({{p(x,z), (Overscript[→, marginal]), p(x)}, {↓,  , ↓}, {S(x,z), (Underscript[→, effective]), S(x)}})

Consider an Ising model of three Ising variables (Subscript[x, 0]),(Subscript[x, 1]),(Subscript[z, 0]) (grouped into x=((Subscript[x, 0]),(Subscript[x, 1])) and z=((Subscript[z, 0]))). Suppose that they interact as

Block[{node},
node[txt_,pos_]:=Translate[{{EdgeForm[Black],White,Disk[{0,0},Offset@{8,8}]},Text["\!\("<>txt<>"\)"]},pos];
Graphics[{{Line[{{-1,0},{0,1},{1,0}}],Text["\!\("<>#1<>"\)",##2]&@@@{{" w",{-0.5,0.5},{1,-0.5}},{" w",{0.5,0.5},{-1,-0.5}}}},node@@@{{"x\_0",{-1,0}},{"x\_1",{1,0}},{"z\_0",{0,1}}}},ImageSize->{Automatic,50}]]

with w=β J>0 being the dimensionless coupling strength (or call the weight in machine learning).

The joint action is

There is no direct interaction between (Subscript[x, 0]),(Subscript[x, 1]).

However, (Subscript[z, 0]) can mediate an effective interaction between (Subscript[x, 0])and (Subscript[x, 1]),

S(x)=-log (∑_z)(e^(-S(x,z)))
=-log (∑_((Subscript[z, 0])=±1))(e^(w ((Subscript[x, 0])+(Subscript[x, 1]))(Subscript[z, 0])))
=-log (2 cosh(w ((Subscript[x, 0])+(Subscript[x, 1])))).

Take w=1 for example, the values of S(x) for different configurations of x are as follows

Block[{w=1.},
Style[Grid[ArrayFlatten[{{"",{{"\!\(x\_1\)",SpanFromLeft},{"\!\(+1\)","\!\(-1\)"}}},{{{"\!\(x\_0\)","\!\(+1\)"},{SpanFromAbove,"\!\(-1\)"}},Table[-Log[2Cosh[w x0+w x1]],{x0,{1,-1}},{x1,{1,-1}}]}}],Alignment->{Center,Center},Dividers->{3->True,3->True}],"InlineFormula"]]

The action is lower if (Subscript[x, 0]) and (Subscript[x, 1]) are aligned => meaning that there is an effective ferromagnetic interaction between (Subscript[x, 0]) and (Subscript[x, 1]), even if there was no direct interaction between them in the original model.

In fact, for Ising variables (Subscript[x, 0]),(Subscript[x, 1]), Eq. () can be rewritten as

Show that Eq. () is consistent with Eq. ().

for (Subscript[x, 0])=(Subscript[x, 1])=±1,

S(x)=-log (2 cosh(2w))=-(Subscript[w, eff])+const.

for (Subscript[x, 0])=-(Subscript[x, 1])=±1,

S(x)=-log 2=+(Subscript[w, eff])+const.

log cosh(2 w)=2(Subscript[w, eff]),

Plot[Log[Cosh[2 w]]/2,{w,0,2},PlotStyle->Black,FrameLabel->{"\!\( w\)","\!\(w\_eff\)"}]

In conclusion, the latent spin (Subscript[z, 0]) can mediate an effective interaction between (Subscript[x, 0]) and (Subscript[x, 1]) with an effective interaction strength (Subscript[w, eff]) that increases with w.

Block[{node},
node[txt_,pos_]:=Translate[{{EdgeForm[Black],White,Disk[{0,0},Offset@{8,8}]},Text["\!\("<>txt<>"\)"]},pos];
Graphics[{Translate[{{Line[{{-1,0},{0,1},{1,0}}],Text["\!\("<>#1<>"\)",##2]&@@@{{" w",{-0.5,0.5},{1,-0.5}},{" w",{0.5,0.5},{-1,-0.5}}}},node@@@{{"x\_0",{-1,0}},{"x\_1",{1,0}},{"z\_0",{0,1}}}},{-2,0}],Text["\[Implies]",{0,0}],Translate[{{Dashed,Line[{{-1,0},{1,0}}],Text["\!\("<>#1<>"\)",##2]&@@@{{" w\_eff",{0,0},{0,-1}}}},node@@@{{"x\_0",{-1,0}},{"x\_1",{1,0}}}},{2,0}]},ImageSize->{Automatic,50}]]

Boltzmann Machine

Random variables:

x={(Subscript[x, i])} - input variables (describing a sample of data, such as an image).

z={(Subscript[z, i])} - latent variables (introduced to mediate effective interactions/correlations among input variables).

Consider (Subscript[x, i]),(Subscript[z, i])=0,1 as binary variables (like occupation numbers in a lattice gas model), such that each sample of x can represent a black-white image (as particles arranged into a pattern on the lattice).

Block[{img},
img=Reverse@Round@ImageData[ColorNegate@First@ResourceData["MNIST", "TrainingData"][[1]]];Graphics[{Text["\!\(x=\)",{0,14},{1.5,0}],Raster[img,{{0,0},Dimensions@img}],{Opacity[0.5],AbsoluteThickness[1],Gray,{Table[Line[{{0,y},{#1,y}}],{y,0,#2}],Table[Line[{{x,0},{x,#1}}],{x,0,#1}]}&@@Dimensions@img}},ImageSize->150]]

(Subscript[p, θ])(x) - the model distribution to optimize.

(Subscript[S, θ])(x,z) - a joint action parametrized by a set of parameters θ (such as the coupling strengths). For example,

(Subscript[S, θ])(x,z)=-(∑_(i,j))(Subscript[W, i j])(Subscript[x, i])(Subscript[z, j])-(∑_i)(Subscript[b, i]^x)(Subscript[x, i])-(∑_j)(Subscript[b, j]^z)(Subscript[z, j])
=-x·W·z-(b^x)·x-(b^z)·z.

Graph[UndirectedEdge@@@Tuples[{StringTemplate["\!\(x\_`1`\)"]/@Range[0,4],StringTemplate["\!\(z\_`1`\)"]/@Range[0,3]}],GraphLayout->"BipartiteEmbedding",EdgeStyle->Directive[Opacity[1],AbsoluteThickness[1.2],Black],VertexShapeFunction->({{EdgeForm[Directive[Opacity[1],Black]],White,Disk[#1,Offset@{8,8}]},{Black,Text[Style[#2,"InlineFormula"],#1]}}&),AspectRatio->1.5,ImageSize->100]

The weights (Subscript[w, i j]) and biases (Subscript[b, i]^x),(Subscript[b, j]^z) are the model parameters, collectively denoted as θ=(W,(b^x),(b^z)).

The setting that the x and z only interact with each other, but not within themselves, is to make the training more efficient by Gibbs sampling. This special variant is called the restricted Boltzmann machine (RBM).

Sample generation: generate new samples x by drawing from the distribution

(Subscript[p, θ])(x)=(1/(Subscript[Z, θ]))(∑_z)(e^(-(Subscript[S, θ])(x,z))).

Other task-specific applications

Anomaly detection

Data imputation

Feature extraction (representation learning)

Classification

Demonstration

SeedRandom[42];
data=;
Grid[Partition[data,4]]

Step 1: write a code to define the restricted Boltzmann machine (RBM) model using PyTorch package.

There are 28×28=784 input variables in x.

We choose to use 32 latent variables for z.

import torch
import numpy
torch.manual_seed(12) # set random seed
# define the Boltzmann machine model
class RBM(torch.nn.Module):
    def __init__(self, nx=784, nz=32):
        super().__init__()
        self.nx = nx
        self.nz = nz
        self.w = torch.nn.Parameter(5*torch.randn((nx, nz)))  # adjusted weight initialization
        self.bx = torch.nn.Parameter(torch.zeros(nx))  # input variables bias
        self.bz = torch.nn.Parameter(torch.zeros(nz))  # latent variables bias

    def x_to_z(self, x, beta=1):
        logit_z = torch.nn.functional.linear(x, self.w.T, self.bz)
        return torch.bernoulli(torch.sigmoid(beta*logit_z))

    def z_to_x(self, z, beta=1):
        logit_x = torch.nn.functional.linear(z, self.w, self.bx)
        return torch.bernoulli(torch.sigmoid(beta*logit_x))

    def forward(self, x, steps=5, beta=1):
        for _ in range(steps):
            z = self.x_to_z(x, beta=beta)
            x = self.z_to_x(z, beta=beta)
        return x
    
    def eff_action(self, x):
        bxx = x.matmul(self.bx)
        logit_z = torch.nn.functional.linear(x, self.w.T, self.bz)
        logZ = -torch.nn.functional.logsigmoid(-logit_z).sum(-1)
        return (- bxx - logZ).mean()
    
    def loss(self, x, steps=5):
        x1 = self(x, steps=steps)
        return self.eff_action(x) - self.eff_action(x1)

Step 2: setup the model, prepare the data, and connect the model to an optimizer.

rbm = RBM()
data = torch.tensor(numpy.array(<*ImageData/@data*>), dtype=torch.float).view(-1, 28*28)
optimizer = torch.optim.Adam(rbm.parameters(), lr=0.1, weight_decay=0.002)

Step 3: train the model on the data for 1000 steps.

loss_sum = 0.
for k in range(1000):
	optimizer.zero_grad()
	loss = rbm.loss(data, steps=5)
	loss.backward()
	optimizer.step()
	loss_sum += loss.item()
	if (k+1)%100==0:
		print(f'step {k+1:4}: loss = {loss_sum/100:6.2f}')
		loss_sum = 0.

The loss function (Subscript[ℒ, θ]) decreases significantly, then gradually converges to the minimum within the error range of stochastic fluctuations.

Start with an initial input x, follow the Markov chain Monte Carlo (MCMC) approach to obtain a sequence of new samples

x→x'→x''→…

according to the model distribution (Subscript[p, θ])(x).

Recall that every action (Subscript[S, θ])(x,z)=β (Subscript[E, θ])(x,z) can be viewed as an rescaling of its energy function by an overall scale β. Tuning this temperature parameter T=1/β controls the strength of thermal fluctuations in the probability model.

DynamicModule[{eval,sample},
eval=ExternalEvaluate[First@ExternalSessions[],#]&;
eval["sample=data.clone()"];
sample[β_]:=Image/@eval[StringTemplate["sample=rbm(sample,beta=`1`);sample.view(-1,28,28).tolist()"][β]];
Manipulate[Style[Dynamic[Column[{Deploy@Row[{Style[StringTemplate["\!\(T=`1`\)"][NumberForm[T,{3,2}]],"InlineFormula"],EventHandler[Style["reset","Text",Underlined],{"MouseClicked":>eval["sample=data.clone()"]}]},Spacer[150]],Grid[Partition[Row/@Thread@{data,sample[1/T]},4]]}]],DynamicUpdating->Dynamic[updating]],{{T,0.1},0.01,2.5,ImageSize->Small},{{updating,True},{False,True}},Paneled->False,LabelStyle->"InlineFormula",SaveDefinitions->True]]

Each digits is a (meta)-stable phase of the lattice gas model.

Raising and lowering the temperature can anneal the model to its stable phase.

(Subscript[S, θ])(x,z)=-x·W·z-(b^x)·x-(b^z)·z.

(b^x) - the chemical potential for the lattice gas.

W - the pattern of x interacting with each of the latent variables (there are 32 of them) in z.

Block[{eval,mplt},
eval=ExternalEvaluate[First@ExternalSessions[],#]&;
mplt=MatrixPlot[#,FrameTicks->None,ImageSize->40,PlotRangePadding->None]&;
Style[Column[{Row[{"\!\(b\^x=\) ",mplt@eval["rbm.bx.view(28,28).tolist()"]}],
Row[{"\!\(W=\) ",Grid[Partition[mplt/@Transpose[eval["rbm.w.view(28,28,-1).tolist()"],{2,3,1}],8]]}]}],"InlineFormula"]]

(Subscript[p, θ])(z|x)=((e^(-(Subscript[S, θ])(x,z)))/((∑_x)(e^(-(Subscript[S, θ])(x,z))))).

z serves as an encoding of x in the latent (feature) space.

By projecting the high-dimensional latent space to a two-dimensional plane, we can visualize the arrangement of hand-written digits x in the latent space: we expect similar digits to be close to each other.

Block[{eval,code,pos},
eval=ExternalEvaluate[First@ExternalSessions[],#]&;
code=eval["rbm.x_to_z(data, beta=10).tolist()"];
pos=DimensionReduce[code].{{1,0},{0,-1}};
Graphics[{Translate[{AbsoluteThickness[1],Line[Offset/@ReIm@{0,15Exp[# I π/2]}]&/@{0,1},Text["\!\( z\)",{0,0},{-1.2,-0.5}]},{-6.5,-5}],Inset[Tooltip[ColorReplace[ColorNegate@#1,White],#1],#2,Automatic,{0.5,0.5}]&@@@Thread@{data,pos}},PlotRangePadding->Scaled[0.05],Frame->True,FrameTicks->None,ImageSize->240]]

Coarse Graining

Problem Setup

Block[{site,sublatt,cb,fine, coarse},
site[col_,pos_]:={EdgeForm[col],Lighter[col,0.8],Disk[pos,0.25]};
sublatt=GroupBy[Function[{x,y},{{x,y},Which[OddQ[x]&&OddQ[y],Red,EvenQ[x]&&EvenQ[y],Green,True,Blue]}]@@@Tuples[Range[-4,4],2],Last->First];
sublatt[Red]=DeleteCases[sublatt[Red],{x_,y_}/;EvenQ[(x+y)/2]];
cb={Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],#]&/@sublatt[Red]};
Graphics[{Line[Cases[Tuples[sublatt[Green],2],{p1_,p2_}/;Norm[p1-p2]==2]],
Table[site[Black,#]&/@sublatt[c],{c,{Green}}],Text["\!\( i\)",{0,0},{-1.5,1}],Text["\!\( j\)",{2,0},{-1.5,1}]},PlotRangePadding->0,ImageSize->{Automatic,120}]]

Degree of freedom: s={(Subscript[s, i])}, each Ising spin (Subscript[s, i])=±1 is a binary variable on site i.

The energy model (without external magnetic field)

E(s)=-J(∑_(⟨i j⟩))(Subscript[s, i])(Subscript[s, j]).

Or the action model

with w=β J=J/T being the dimensionless coupling constant.

Small w: weak coupling, high temperature;

Large w: strong coupling, low temperature.

Probability distribution (Boltzmann distribution)

Block Partition: partitioning spins into blocks.

Representation Learning: emergence of coarse-grained (higher-level) spins.

Decimation: summing over (marginalize) fine-grained (lower-level) spins to obtain an effective theory for coarse-grained (higher-level) spins.

Block Partition

Block[{site,sublatt,cb,fine, coarse},
site[col_,pos_]:={EdgeForm[col],Lighter[col,0.8],Disk[pos,0.25]};
sublatt=GroupBy[Function[{x,y},{{x,y},Which[OddQ[x]&&OddQ[y],Red,EvenQ[x]&&EvenQ[y],Green,True,Blue]}]@@@Tuples[Range[-4,4],2],Last->First];
sublatt[Red]=DeleteCases[sublatt[Red],{x_,y_}/;EvenQ[(x+y)/2]];
cb={Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],#]&/@sublatt[Red]};
Graphics[{cb,Line[Cases[Tuples[sublatt[Green],2],{p1_,p2_}/;Norm[p1-p2]==2]],
site[Black,#]&/@sublatt[Green],
Text["\!\( i\)",{0,0},{-1.5,1}],Text["\!\( j\)",{2,0},{-1.5,1}],Text["\!\( I\)",{1,1}]},PlotRangePadding->0,ImageSize->{Automatic,120}]]

I - block index (larger scale).

i,j - site index (smaller scale).

(∑_(⟨i j⟩))=(∑_I)(∑_(⟨i j⟩∈I)).

p(s)∝(e^(-S(s)))=(∏_I)(e^(-(Subscript[S, I])(s)))∝(∏_I)(Subscript[p, I])(s).

Note: Eq. () does not mean the distribution can be factored into independent distributions in each block, because the spins are shared between blocks — not independent variables.

Representation Learning

Block[{site,sublatt,cb,fine, coarse},
site[col_,pos_]:={EdgeForm[col],Lighter[col,0.8],Disk[pos,0.25]};Graphics[{Translate[{{Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],{1,1}]},Line[{{0,0},{2,0},{2,2},{0,2},{0,0}}],
site[Blue,#]&/@{{0,0},{2,0},{2,2},{0,2}},Text["\!\( i\)",{0,0},{0,1.2}],Text["\!\( j\)",{2,0},{0,1.2}],Text["\!\( I\)",{1,1}]},{-2,0}],
Text["→",{1,1}],
Translate[{{Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],{1,1}]},Line[{{1,1},#}&/@{{0,0},{2,0},{2,2},{0,2}}],
site[Blue,#]&/@{{0,0},{2,0},{2,2},{0,2}},site[Orange,#]&/@{{1,1}},Text["\!\( i\)",{0,0},{0,1.2}],Text["\!\( I\)",{1,1},{0,-1}]},{2,0}]},ImageSize->90]]

Spin degrees of freedom:

s={(Subscript[s, i])} with (Subscript[s, i])=±1 - the original fine-grained spin (in blue),

s'={(Subscript[s, I]^′)} with (Subscript[s, I]^′)=±1 - the emergent coarse-grained spin (in orange), as a representative of the fine-grained spins.

The local action

(Subscript[S, I])(s)=-w(∑_(⟨i j⟩∈I))(Subscript[s, i])(Subscript[s, j]) ,

defines a local distribution as a target distribution for a generative model to learn

in which (Subscript[s, i]) interact with each other.

Introduce a generative model with the following action

where

(Subscript[s, i]) no longer interacts with each other directly, but they all interact with an emergent representative spin (Subscript[s, I]^′) at the block center. — the architecture of a restricted Boltzmann machine.

(OverTilde[w]) is the model parameter to be optimized. It controls how strong the representative spin s' interacts with the underlying spins s, described by the conditional distribution

(Subscript[(OverTilde[p]), I])(s'|s)=((e^(-(Subscript[(OverTilde[S]), I])(s,s')))/((∑_s)(e^(-(Subscript[(OverTilde[S]), I])(s,s'))))).

Manipulate[
Column[{Style[StringTemplate["\!\(w\&~=`1`\)"][NumberForm[w,{3,2}]],"InlineFormula"],Grid[Prepend[{Style[Row[<|-1->Style["-",Blue],+1->Style["+",Red]|>/@#],"CodeText",Bold],Block[{p1,p2},{p1,p2}=#/Total[#]&@Exp[w Total[#]{-1,1}];
Graphics[{{Blue,Text[Style["-","CodeText",Bold],{0,0.5},{1.2,0}],Rectangle[{0,0},{p1,1}]},{Red,Text[Style["+","CodeText",Bold],{1,0.5},{-1.2,0}],Rectangle[{p1,0},{1,1}]},{White,Translate[Rectangle[{-0.01,0},{0.01,1}],{p1-0.005,0}]}},AspectRatio->0.18,ImageSize->{Automatic,10}]]}&/@SortBy[Tuples[{-1,1},4],Total],{"\!\( s\)","\!\(p\&~\_I(s\^′|s)\)"}],ItemStyle->"InlineFormula",Dividers->{{2->True},{2->True,3->True,7->True,13->True,17->True}}]}],{{w,0.5},0,1,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Given (Subscript[(OverTilde[S]), I])(s,s'), the model distribution (Subscript[(OverTilde[p]), I])(s) of the original spins s is

Comparing the target distribution (Subscript[p, I])(s) (parametrized by w) and the model distribution (Subscript[(OverTilde[p]), I])(s) (parametrized by (OverTilde[w])) for all configurations of s:

Manipulate[Block[{s,ptgt,pmdl,bar,surprise},
s=Tuples[{-1,1},4];
ptgt=#/Total[#]&@Exp[w0 #.RotateLeft[#]&/@s];
pmdl=#/Total[#]&@Sum[Exp[w1 Total[#]sp&/@s],{sp,{-1,1}}];
bar[p_]:=Graphics[{{GrayLevel[0.6],Rectangle[{0,-1},{p,1}]},{Opacity[0],Rectangle[{p,-1},{0.5,1}]},{Text[Style[NumberForm[Round[p,0.001],{4,3}],12],{0.5,0},{1,0}]}},PlotRangePadding->0,AspectRatio->0.15,ImageSize->{Automatic,12}];
surprise=ptgt Log[ptgt/pmdl];
Grid[Join[{{"",StringTemplate["\!\(w=`1`\)"][NumberForm[w0,{3,2}]],StringTemplate["\!\(w\&~=`1`\)"][NumberForm[w1,{3,2}]]},{"\!\( s\)","\!\(p\_I(s)\)","\!\(p\&~\_I(s)\)"}},{Style[Row[<|-1->Style["-",Blue],+1->Style["+",Red]|>/@#1],"CodeText",Bold],bar[#2],bar[#3]}&@@@Thread@{s,ptgt,pmdl}],ItemStyle->"InlineFormula",Dividers->{{2->True,3->True},{3->True}},
Background->{None,Join[{White,White},Blend[{Green,White,Purple},Rescale[#,{-1,1}]]&/@surprise]}]
],{{w0,0.72,"\!\( w\)"},0,2.5,ImageSize->Small},{{w1,1.18,"\!\(w\&~\)"},0,2.5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Loss function: the KL divergence

which is a function of (OverTilde[w]) given w. For each w, find the optimal (OverTilde[w]) to minimize the KL divergence ℒ.

Block[{s,ptgt,pmdl,L,rec},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
L[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]]|>&@@FindMinimum[L[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{w,0.01,2.5,0.01}];
Column[{ListLinePlot[Prepend[{#["w"],#["wt"]}&/@rec,{0,0}],PlotStyle->Black,FrameLabel->{"\!\( w\)","\!\(w\&~\)"},ImagePadding->{{48, 5}, {2, 2}},ImageSize->180],
ListLinePlot[Prepend[{#["w"],#["L"]/Log[2]}&/@rec,{0,0}],PlotStyle->Red,FrameLabel->{"\!\( w\)","\!\(ℒ\) [bit]"},
ImagePadding->{{48, 5}, {40, 2}},ImageSize->180]},Spacings->0]
]

The optimization is not perfect given the non-vanishing KL divergence, indicating that the model distribution (Subscript[(OverTilde[p]), I])(s) approximates rather than precisely matches the target distribution (Subscript[p, I])(s).

However, the small residual KL divergence (<0.1 bit) suggests that the approximation is sufficiently accurate.

The approximation can be improved as more representative (latent) spins are introduced, see Ref. [].

As a result, we obtain the optimal model parameter (OverTilde[w])(w) as a function of the original parameter w.

Block[{site,sublatt,cb,fine, coarse},
site[col_,pos_]:={EdgeForm[col],Lighter[col,0.8],Disk[pos,0.25]};Graphics[{Translate[{{Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],{1,1}]},Line[{{0,0},{2,0},{2,2},{0,2},{0,0}}],Text["\!\( w\)",{1,2},{0,-0.5}],
site[Blue,#]&/@{{0,0},{2,0},{2,2},{0,2}},Text["\!\( i\)",{0,0},{0,1.2}],Text["\!\( j\)",{2,0},{0,1.2}],Text["\!\( I\)",{1,1}]},{-2,0}],
Text["→",{1,1}],
Translate[{{Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],{1,1}]},Line[{{1,1},#}&/@{{0,0},{2,0},{2,2},{0,2}}],
Text["\!\( w\&~\)",{1.5,0.5},{-0.5,-0.5}],
site[Blue,#]&/@{{0,0},{2,0},{2,2},{0,2}},site[Orange,#]&/@{{1,1}},Text["\!\( i\)",{0,0},{0,1.2}],Text["\!\( I\)",{1,1},{0,-1}]},{2,0}]},ImageSize->100]]

The optimization result tells us how to treat the local action (Subscript[S, I])(s) as an effective action of (Subscript[(OverTilde[S]), I])(s,s') approximately.

Wanda Hou, Yi-Zhuang You. Machine Learning Renormalization Group for Statistical Physics. arXiv:2306.11054.

Decimation

Block[{site,sublatt,cb,fine, coarse},
site[col_,pos_]:={EdgeForm[col],Lighter[col,0.8],Disk[pos,0.25]};
sublatt=GroupBy[Function[{x,y},{{x,y},Which[OddQ[x]&&OddQ[y],"IR",EvenQ[x]&&EvenQ[y],"VB",True,"UV"]}]@@@Tuples[Range[-4,4],2],Last->First];
sublatt["IR"]=DeleteCases[sublatt["IR"],{x_,y_}/;EvenQ[(x+y)/2]];
cb={Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],#]&/@sublatt["IR"]};
Graphics[{Translate[{cb,Line[Cases[Tuples[sublatt["VB"],2],{p1_,p2_}/;Norm[p1-p2]==2]],
site[Blue,#]&/@sublatt["VB"],Text["\!\(s\_i\)",{0,0},{-1.5,0.8}]},{-5.5,0}],Text["→",{0,0}],
Translate[{cb,Line[Cases[Tuples[{sublatt["IR"],sublatt["VB"]}],{p1_,p2_}/;Norm[p1-p2]==Sqrt[2]]],
site[Blue,#]&/@sublatt["VB"],site[Orange,#]&/@sublatt["IR"],Text["\!\(s\_i\)",{0,0},{-1.5,0.8}],Text["\!\(s\_I\%′\)",{1,1},{-2,0}]},{5.5,0}]},PlotRangePadding->0,ImageSize->{Automatic,120}]]

If we marginalize the emergent coarse-grained spins s', we fall back to the original action for the fine-grained spins s approximately

(∑_(s'))(e^(-(OverTilde[S])(s,s')))≃(e^(-S(s))).

However, if we marginalize the fine-grained spins s, we will obtain a new effective action S'(s') for the coarse-grained spins s'

Block[{site,sublatt,cb,fine, coarse},
site[col_,pos_]:={EdgeForm[col],Lighter[col,0.8],Disk[pos,0.25]};
sublatt=GroupBy[Function[{x,y},{{x,y},Which[OddQ[x]&&OddQ[y],"IR",EvenQ[x]&&EvenQ[y],"VB",True,"UV"]}]@@@Tuples[Range[-4,4],2],Last->First];
sublatt["IR"]=DeleteCases[sublatt["IR"],{x_,y_}/;EvenQ[(x+y)/2]];
cb={Opacity[0.1],Translate[Rectangle[{-1.3,-1.3},{1.3,1.3},RoundingRadius->0.3],#]&/@sublatt["IR"]};
Graphics[{Translate[{Line[Cases[Tuples[{sublatt["IR"],sublatt["VB"]}],{p1_,p2_}/;Norm[p1-p2]==Sqrt[2]]],
site[Blue,#]&/@sublatt["VB"],site[Orange,#]&/@sublatt["IR"],Text["\!\(s\_i\)",{0,0},{-2,0}],Text["\!\(s\_I\%′\)",{-1,-1},{-2,0}],Text["\!\(s\_J\%′\)",{1,1},{-2,0}],Text["\!\(w\&~\)",{-0.5,-0.5},{1,-0.5}]},{-5.5,0}],Text["→",{0,0}],
Translate[{Line[Cases[Tuples[{sublatt["IR"],sublatt["VB"]}],{p1_,p2_}/;Norm[p1-p2]==Sqrt[2]]],
site[Orange,#]&/@sublatt["IR"],Text["\!\(s\_I\%′\)",{-1,-1},{-2,0}],Text["\!\(s\_J\%′\)",{1,1},{-2,0}],Text["\!\(w\^′\)",{0,0},{0.5,-0.5}]},{5.5,0}]},PlotRangePadding->0,ImageSize->{Automatic,120}]]

The fine-grained spins s mediate effective interactions between the emergent coarse-grained spins s',

with an effective coupling w' as a function of the original coupling w by (recall Eq. ())

where (OverTilde[w])(w) is given by the following curve.

Block[{s,ptgt,pmdl,L,rec},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
L[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]]|>&@@FindMinimum[L[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{w,0.01,2.5,0.01}];
ListLinePlot[Prepend[{#["w"],#["wt"]}&/@rec,{0,0}],PlotStyle->Black,FrameLabel->{"\!\( w\)","\!\(w\&~(w)\)"},ImagePadding->{{48, 5}, {40, 2}},ImageSize->180]
]

If we rotate the lattice by 45° and rescale the lattice spacing by (Sqrt[2]), the effective action S'(s') in Eq. () maps to the original action S(s) in Eq. (), with the coupling modified by w->w' in Eq. ().

({{p(s),  , (OverTilde[p])(s,s'), →, p'(s')}, {↓,  , ↑,  , ↓}, {S(s),  , (OverTilde[S])(s,s'),  , S'(s')}, {=(∑_I)(Subscript[S, I])(s), →, =(∑_I)(Subscript[(OverTilde[S]), I])(s,s'),  , …}})

Every step is exact except for the representation learning, where the KL divergence is non-vanishing (i.e. the learning is not perfect).

Under coarse graining,

the lattice up-scales by (Sqrt[2]), i.e. the lattice spacing rescaled by ℓ->ℓ'=(Sqrt[2])ℓ,

the action restores the same form (this doesn’t need to be so in general, but let us assume it for now),

the coupling constant is updated by w→w' (more generally, it is the generative model that gets updated).

Renormalization Group Flow

RG Flow Equation

(w^((0)))->(w^((1)))->(w^((2)))->(w^((3)))->…,

ℛ denotes the RG transformation, mapping the coupling constant w from one scale to another.

Its specific form in Eq. () is derived for 2D Ising model under the given RG scheme (i.e., block size, number of representative spins, parameter space dimension, etc.). ℛ is expected to be different for different models and RG schemes.

Block[{s,ptgt,pmdl,L,rec,wfun,wc},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
L[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]],"w1"->Log[Cosh[2 #2[[1,2]]]]/2|>&@@FindMinimum[L[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{w,0.01,2.5,0.01}];
wfun=Interpolation[Prepend[{#["w"],#["w1"]}&/@rec,{0,0}]];
wc=w/.FindRoot[wfun[w]==w,{w,0.4}];
Plot[{wfun[w],w},{w,0,1},PlotStyle->{Black,Directive[AbsoluteThickness[1],Dashed,Black]},Frame->True,FrameLabel->{"\!\( w\^(k)\)","\!\(w\^(k+1)\)"},PlotRange->{{0,1},{0,1}},GridLines->{{wc},{wc}},Epilog->{AbsoluteThickness[1],Arrowheads[0.04],Red,Arrow[Flatten[{{#1,#2},{#2,#2}}&@@@Thread[{Most[#],Rest[#]}]&@NestList[wfun,wc+0.001,18],1]],
Blue,Arrow[Flatten[{{#1,#2},{#2,#2}}&@@@Thread[{Most[#],Rest[#]}]&@NestList[wfun,wc-0.001,22],1]],Black,Text[StringTemplate["\!\(w\_c≈`1`\)"][NumberForm[wc,3]],{wc,wc},{-1,1}]},AspectRatio->Automatic,ImageSize->220]
]

w=β J=J/T.

Small w (weak coupling, high temperature) - disordered phase. Under coarse graining:

Block[{s,ptgt,pmdl,L,rec,wfun,wc,bern,upsample,d=64},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
L[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]],"w1"->Log[Cosh[2 #2[[1,2]]]]/2|>&@@If[w<=5,FindMinimum[L[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{0,{wt->2w-0.3462}}],{w,Join[Range[0.01,2.5,0.01],Range[3,10],Range[20,90,10],10^Range[2,5]]}];
wfun=Interpolation[Prepend[{#["w"],#["w1"]}&/@rec,{0,0}]];
wc=w/.FindRoot[wfun[w]==w,{w,0.4}];
bern[h_]:=Sign[Tanh[h]-RandomReal[{-1,1},Dimensions[h]]];
upsample[x_,{w0_,w1_}]:=Block[{zs,xs,xraw},
zs=bern[w0(x+RotateLeft[x,#])&/@{{0,-1},{-1,0},{0,1},{1,0}}];
xs=bern[w1(zs+RotateLeft[zs])];
xraw=Flatten[{{xs[[1]],xs[[2]]},{xs[[4]],xs[[3]]}},{{3,1},{4,2}}];
If[Length[xraw]<=d,xraw,xraw[[;;d,;;d]]]];
SeedRandom[3];
Row[Framed[Image[#,ImageSize->d],FrameMargins->None,FrameStyle->AbsoluteThickness[1]]&/@Reverse[Cases[FoldList[upsample,{{1}},Partition[ArcCosh[Exp[2#]]/2&/@Reverse@NestList[wfun,wc-0.03,22],2]],x_/;Length[x]==d]],"→"]
]

Spins are randomly oriented due to strong thermal fluctuations. Their orientations averages out, allowing the representative spin to fluctuate more freely, weakening spin correlation at larger scales.

Large w (strong coupling, low temperature) - ordered phase. Under coarse graining:

Block[{s,ptgt,pmdl,L,rec,wfun,wc,bern,upsample,d=64},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
L[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]],"w1"->Log[Cosh[2 #2[[1,2]]]]/2|>&@@If[w<=5,FindMinimum[L[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{0,{wt->2w-0.3462}}],{w,Join[Range[0.01,2.5,0.01],Range[3,10],Range[20,90,10],10^Range[2,5]]}];
wfun=Interpolation[Prepend[{#["w"],#["w1"]}&/@rec,{0,0}]];
wc=w/.FindRoot[wfun[w]==w,{w,0.4}];
bern[h_]:=Sign[Tanh[h]-RandomReal[{-1,1},Dimensions[h]]];
upsample[x_,{w0_,w1_}]:=Block[{zs,xs,xraw},
zs=bern[w0(x+RotateLeft[x,#])&/@{{0,-1},{-1,0},{0,1},{1,0}}];
xs=bern[w1(zs+RotateLeft[zs])];
xraw=Flatten[{{xs[[1]],xs[[2]]},{xs[[4]],xs[[3]]}},{{3,1},{4,2}}];
If[Length[xraw]<=d,xraw,xraw[[;;d,;;d]]]];
SeedRandom[11];
Row[Framed[Image[#,ImageSize->d],FrameMargins->None,FrameStyle->AbsoluteThickness[1]]&/@Reverse[Cases[FoldList[upsample,{{1}},Partition[ArcCosh[Exp[2#]]/2&/@Reverse@NestList[wfun,wc+0.01,22],2]],x_/;Length[x]==d]],"→"]
]

Spins are locally aligned, reinforcing the orientation of local representative spins, mediating in a mutual alignment of representative spins, enhancing spin correlation at larger scales.

The two phases are separated by a critical point at w=(Subscript[w, c]), where the phase transition happens.

At the critical point, w stops flowing — realizing a RG fixed point. The critical value (Subscript[w, c]) can be found by solving the fixed point equation

where the RG transformation ℛ was defined in Eq. (). One solution of Eq. () is

(Subscript[w, c])≈0.486,

which is close to the known exact solution of 2D Ising model

(Subscript[w, c])=(1/2)log(1+(Sqrt[2]))≈0.4407.

At the RG fixed point, the probability model p(s) looks the same at different scales — an emergent scale invariance. Under coarse graining:

Block[{s,ptgt,pmdl,L,rec,wfun,wc,bern,upsample,d=64},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
L[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]],"w1"->Log[Cosh[2 #2[[1,2]]]]/2|>&@@If[w<=5,FindMinimum[L[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{0,{wt->2w-0.3462}}],{w,Join[Range[0.01,2.5,0.01],Range[3,10],Range[20,90,10],10^Range[2,5]]}];
wfun=Interpolation[Prepend[{#["w"],#["w1"]}&/@rec,{0,0}]];
wc=w/.FindRoot[wfun[w]==w,{w,0.4}];
bern[h_]:=Sign[Tanh[h]-RandomReal[{-1,1},Dimensions[h]]];
upsample[x_,{w0_,w1_}]:=Block[{zs,xs,xraw},
zs=bern[w0(x+RotateLeft[x,#])&/@{{0,-1},{-1,0},{0,1},{1,0}}];
xs=bern[w1(zs+RotateLeft[zs])];
xraw=Flatten[{{xs[[1]],xs[[2]]},{xs[[4]],xs[[3]]}},{{3,1},{4,2}}];
If[Length[xraw]<=d,xraw,xraw[[;;d,;;d]]]];
SeedRandom[7];
Row[Framed[Image[#,ImageSize->d],FrameMargins->None,FrameStyle->AbsoluteThickness[1]]&/@Reverse[Cases[NestList[upsample[#,{#,#}&@(ArcCosh[Exp[2wc]]/2)]&,{{1}},11],x_/;Length[x]==d]],"→"]
]

The spin configuration looks (statistically) the same as we zoom out (or in) — exhibiting a self-similarity. As a consequence, all thermodynamic properties must be scale-free at the critical point.

RG Flow Diagram

Block[{a=4+2I,col,cmap,icmap,c,stream},
col=<|"phase1"->Green,"phase2"->Blue,"critical"->Red|>;
cmap=#(1+#/a)&;
icmap=a(Sqrt[1+4#/a]-1)/2&;
c[{θ1_,θ2_}]:=ComplexExpand[(-#1^2+#2^2+#1^4/2)&@@ReIm[cmap@(θ1+I θ2)]];
stream=First@StreamPlot[Evaluate[-Grad[c[{θ1,θ2}],{θ1,θ2}]],{θ1,-2,2},{θ2,-1,1},StreamPoints->{{{ReIm@icmap[0.1I],col["critical"]},{ReIm@icmap[-0.2I],col["critical"]},{ReIm@icmap[-0.5],col["phase1"]},{ReIm@icmap[0.5],col["phase2"]},Automatic}},StreamStyle->Directive[GrayLevel[0.7],AbsoluteThickness[1]],StreamScale->{0.08,All,0.025}];
Graphics[{Block[{lt=0.8,pts},pts=ReIm[icmap/@(I Range[-0.98,1.16,0.01])];{{Lighter[col["phase1"],lt],Polygon[Join[pts,{{-2,1},{-2,-1}}]]},{Lighter[col["phase2"],lt],Polygon[Join[pts,{{2,1},{2,-1}}]]}}],stream,{AbsolutePointSize[6],{#1,Point[#2]}&@@@Thread@{col/@{"phase1","critical","phase2"},ReIm[icmap/@{-1,0,1}]}}},PlotRangePadding->0,Frame->True,FrameLabel->{"\!\(θ\_1\)","\!\(θ\_2\)"},AspectRatio->Automatic,ImageSize->{Automatic,160}]]

Flow lines: the trajectories of θ under coarse graining. In the continuum limit, the flow of θ can be described by a differential RG flow equation

where β(θ) describes the flow velocity at each θ point, and log ℓ measures the RG step (or RG “time”) with ℓ being the cutoff length scale (such as the lattice spacing).

Fixed points: the points (Subscript[θ, ★]) where the flow lines converge or diverge, i.e. where the flow velocity vanishes

β((Subscript[θ, ★]))=0.

Solutions of (Subscript[θ, ★]) represent scale-invariant states of the system, where the physical behavior does not change under renormalization.

Stable fixed points (the green, blue points): attracting flow lines from all directions in their vicinity, representing stable phases of the system that are robust against perturbations.

Unstable fixed points (the red point): repelling flow lines in at least one directions, marking the phase transition between adjacent phases.

Critical points: unstable fixed points with one relevant perturbation (the small perturbation under which θ flows away from the fixed point).

Multi-critical points: unstable fixed points with more than one relevant perturbations.

Critical surface (the red flow line): the line (or surface, more generally) on which the RG flow converges to a critical point. It separates the RG flow towards different stable fixed points, and serves as the phase boundary.

Perturbations at RG Fixed Points

Derive Eq. () from Eq. ().

((dδθ)/(dlog ℓ))=((d(θ-(Subscript[θ, ★])))/(dlog ℓ))=((dθ)/(dlog ℓ))=-β(θ).

β(θ)=β((Subscript[θ, ★])+δθ)
=β((Subscript[θ, ★]))+δθ·(Subscript[∇, θ])β((Subscript[θ, ★]))+𝒪(((δθ)^2))
=0+δθ·(Subscript[∇, θ])β((Subscript[θ, ★]))+𝒪(((δθ)^2)).

((dδθ)/(dlog ℓ))=-δθ·(Subscript[∇, θ])β((Subscript[θ, ★]))+𝒪(((δθ)^2)).

δθ=(∑_a)(Subscript[δθ, a]) (Subscript[n, a]),

(Subscript[n, a])·(Subscript[∇, θ])β((Subscript[θ, ★]))=(Subscript[λ, a])(Subscript[n, a]),

Derive Eq. ().

((d(Subscript[δθ, a]))/(dlog ℓ))(Subscript[n, a])=((d((Subscript[δθ, a])(Subscript[n, a])))/(dlog ℓ))
=-((Subscript[δθ, a])(Subscript[n, a]))·(Subscript[∇, θ])β((Subscript[θ, ★]))+𝒪(((δθ)^2))
=-(Subscript[δθ, a])(Subscript[n, a])·(Subscript[∇, θ])β((Subscript[θ, ★]))+𝒪(((δθ)^2))
=-(Subscript[λ, a])(Subscript[δθ, a])(Subscript[n, a])+𝒪(((δθ)^2)).

Relevant perturbation: an eigen perturbation (Subscript[δθ, a]) with (Subscript[λ, a])<0, that increases under RG and becomes important at large scale, serving as a driving parameter of the phase transition, perpendicular to the critical surface.

Irrelevant perturbation: an eigen perturbation (Subscript[δθ, a]) with (Subscript[λ, a])>0,  that diminishes under RG and become insignificant at large scale. Irrelevant perturbations span the tangent space of the critical surface.

Marginal perturbation: an eigen perturbations (Subscript[δθ, a]) with (Subscript[λ, a])=0. There is no RG flow under such perturbation — a delicate balance at the threshold between relevance and irrelevance.

Exact marginal perturbation: all higher order terms in the RG flow equation also vanish, such that d(Subscript[δθ, a])/dlog ℓ=0 exactly.

Marginally relevant/irrelevant perturbation: if the higher order terms do not vanish, they will decide whether the perturbation is relevant (increasing) or irrelevant (diminishing) under RG flow.

Universal Scaling Behaviors

Scaling Hypothesis

The only function that is scale-free is the power-law function. For example

y=A (ℓ^(-λ)),

any change in the length scale ℓ→b ℓ can be absorbed into the coefficient A, leaving the form of the function unchanged.

Manipulate[Plot[1/ℓ^λ,{ℓ,0,L},PlotStyle->Black,FrameLabel->{"\!\(ℓ\)","\!\( y\)"},ImagePadding->{{60, 8}, {40, 5}},PlotLabel->StringTemplate["\!\(y∼ℓ\^`1`\)"][NumberForm[-λ,{3,2}]],ImageSize->{Automatic,140}],{{λ,1},-3,3,ImageSize->Small},{{L,1,"scale"},0.1,10,ImageSize->Small},LabelStyle->"InlineFormula",Paneled->False]

All thermodynamic properties are derivatives of free energy density f=F/V (the free energy per site) with respect to model parameters θ (or Lagrangian multipliers). At the critical point, they scale with the length scale ℓ as

where d is the spatial dimension, and (Subscript[λ, a]) is the RG eigenvalue associated with the eigen perturbation (Subscript[δθ, a]).

The free energy density scaling f∼(ℓ^d) follows from the fact that a single site at scale ℓ is a representative of (ℓ^d) microscopic spins, which has a free energy that is (ℓ^d) times the average free energy of each microscopic spin.

The parameter scaling (Subscript[δθ, a])∼(ℓ^(-(Subscript[λ, a]))) was given by the solution of the linearized RG equation near the critical point, see Eq. ().

To make all thermodynamic properties scale invariant, the free energy density f as a function of (Subscript[δθ, a]) must take a scale-free form, described by the scaling hypothesis

ϕ is a universal scaling function relating the dimensionless parameters (Subscript[δθ, a])/(ℓ^(-(Subscript[λ, a]))) with the dimensionless free energy f/(ℓ^d), such that the scale ℓ is canceled everywhere, and Eq. () is scale-free.

Scaling Laws

The critical point is at (T,h)=((Subscript[T, c]),0).

Eigen perturbations scale as

δT=T-(Subscript[T, c])∼(ℓ^(-(Subscript[λ, T]))),
δh=h∼(ℓ^(-(Subscript[λ, h]))).

with two exponents (Subscript[λ, T]) and (Subscript[λ, h]) to be determined. We may call (Subscript[λ, T]) the temperature exponent and (Subscript[λ, h]) the magnetic exponent.

The scaling form Eq. () implies

Specific heat

(Subscript[c, V])=-T(((∂^2)f)/(∂(T^2)))∼(ℓ^(d+2(Subscript[λ, T]))).

Magnetization density

m=-((∂f)/(∂h))∼(ℓ^(d+(Subscript[λ, h]))).

Spin susceptibility

χ=-(((∂^2)f)/(∂(h^2)))∼(ℓ^(d+2(Subscript[λ, h]))).

Derive Eq. (), Eq. (), Eq. ().

(ϕ^((m,n)))(x,y):=(∂_x^m)(∂_y^n)ϕ(x,y).

(Subscript[c, V])=-T(((∂^2)f)/(∂(T^2)))
=-(Subscript[T, c])(((∂^2)f)/(∂(((δT))^2)))
=-(Subscript[(Subscript[T, c])(((∂^2)/(∂(((δT))^2)))(ℓ^d)ϕ((ℓ^(Subscript[λ, T]))δT,(ℓ^(Subscript[λ, h]))h)), δT,h->0])
=-(Subscript[T, c])(ϕ^((2,0)))(0,0)(ℓ^(d+2(Subscript[λ, T]))),

m=-((∂f)/(∂h))
=-(Subscript[((∂/(∂h))(ℓ^d)ϕ((ℓ^(Subscript[λ, T]))δT,(ℓ^(Subscript[λ, h]))h)), δT,h->0])
=-(ϕ^((0,1)))(0,0)(ℓ^(d+(Subscript[λ, h]))),

χ=-(((∂^2)f)/(∂(h^2)))
=-(Subscript[(((∂^2)/(∂(h^2)))(ℓ^d)ϕ((ℓ^(Subscript[λ, T]))δT,(ℓ^(Subscript[λ, h]))h)), δT,h->0])
=-(ϕ^((0,2)))(0,0)(ℓ^(d+2(Subscript[λ, h]))).

({{Definition, Exponent}, {(Subscript[c, V])∼((|T-(Subscript[T, c])|)^(-α)), α=(d+2(Subscript[λ, T]))/(Subscript[λ, T])}, {(Subscript[m, ★])∼((((Subscript[T, c])-T))^β), β=-(d+(Subscript[λ, h]))/(Subscript[λ, T])}, {χ∼(((T-(Subscript[T, c])))^(-γ)), γ=(d+2(Subscript[λ, h]))/(Subscript[λ, T])}, {(Subscript[m, ★])∼(h^(1/δ)), δ=-(Subscript[λ, h])/(d+(Subscript[λ, h]))}})

Rushbrooke’s relation []

Widom’s relation []

G. S. Rushbrooke, J. Chem. Phys., 39, 842, (1963).

B. Widom, J. Chem. Phys. 41, 1633 (1964).

Correlation Function

Block[{L=16,r={20,15}},
Graphics[{Translate[{Opacity[0.15],AbsoluteThickness[1.2],Rotate[Translate[Line[{#,0}&/@{-L-0.5,L+0.5}],{0,#}]&/@Range[-L,L],# π/2]&/@{0,1}},{L/2,L/2}],{Dashed,Line[{{0,0},r}]},{Point[{{0,0},r}]},Text["\!\(s(0)\)",{0,0},{0,1}],Text["\!\(s(r)\)",r,{0,-1}]},Frame->True,FrameLabel->{"\!\(r\_1\)","\!\(r\_2\)"}]]

At the critical point, the self-similar and scale-free properties necessitate a power-law correlation

thereby defining an exponent η.

Away from the critical point, a characteristic length scale — the correlation length ξ — emerges, and the correlation function decays exponentially as

G(r)∼(1/((|r|)^(d-2+η)))(e^(-|r|/ξ)).

The correlation length is expected to diverge (ξ→∞) as the system approaches the critical point (for Eq. () to reduce to Eq. ()). The divergence is also a power-law,

with another exponent ν.

Hyperscaling Relations

The exponent ν regards how the correlation length ξ scales.

The correlation length ξ describes the typical size of a magnetic domain (a region in which spins are predominantly aligned) in the spin configuration.

Block[{samples,bg,ah,plt},
samples=Block[{s,ptgt,pmdl,loss,rec,wfun,wc,bern,upsample,d=32},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
loss[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]],"w1"->Log[Cosh[2 #2[[1,2]]]]/2|>&@@If[w<=5,FindMinimum[loss[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{0,{wt->2w-0.3462}}],{w,Join[Range[0.01,2.5,0.01],Range[3,10],Range[20,90,10],10^Range[2,5]]}];
wfun=Interpolation[Prepend[{#["w"],#["w1"]}&/@rec,{0,0}]];
wc=w/.FindRoot[wfun[w]==w,{w,0.4}];
bern[h_]:=Sign[Tanh[h]-RandomReal[{-1,1},Dimensions[h]]];
upsample[x_,{w0_,w1_}]:=Block[{zs,xs,xraw},
zs=bern[w0(x+RotateLeft[x,#])&/@{{0,-1},{-1,0},{0,1},{1,0}}];
xs=bern[w1(zs+RotateLeft[zs])];
xraw=Flatten[{{xs[[1]],xs[[2]]},{xs[[4]],xs[[3]]}},{{3,1},{4,2}}];
If[Length[xraw]<=d,xraw,xraw[[;;d,;;d]]]];
SeedRandom[0];
FoldList[upsample,{{1}},Partition[ArcCosh[Exp[2#]]/2&/@Reverse@NestList[wfun,wc-0.1,10],2]]];
ah=Graphics[Line[Offset/@{{-4,3},{0,0},{-4,-3}}]];
plt[sample_,ξ_,ℓ_]:=Block[{L=Length[sample]/2,x={20,15},m=5,n=5},
Graphics[{{Raster[Rescale[sample,{-1,1}],{{-1,-1},{1,1}},ColorFunction->(GrayLevel[0.85+0.15#]&)]},{Opacity[0.1],AbsoluteThickness[1.2],Rotate[Translate[Line[{#,0}&/@({-L,L}/L)],{0,#}]&/@(Range[-L,L]/L),# π/2]&/@{0,1}},Block[{r=0.5},Translate[{Red,{AbsoluteThickness[1],Dashed,Circle[{0,0},r]},{Arrowheads[{{0.01,1,ah},{-0.01,0,ah}}],Arrow[ReIm[{0,r}Exp[I π/3]]]},Text[StringTemplate["\!\(ξ≃`1`\)"][ξ],ReIm[r/2 Exp[I π/3]],{-1.2,0.5}]},{0.,-0.375}]]},BaselinePosition->Scaled[0.39],PlotRange->({#,#}&@{-1,1}),PlotLabel->Style[StringTemplate["\!\(ℓ=`1`\)"][ℓ],"InlineFormula",Black],PlotRangePadding->0,Frame->True,FrameTicks->None,FrameLabel->None,ImageSize->100]];
Row[plt[samples[[-#]],2^(4-#),2^(#-1)]&/@Range[3],Style["→","Text"]]]

It worth noting that ξ is measured in units of the lattice spacing ℓ. So as ℓ increases under coarse graining, ξ will decrease inversely,

ξ∼(ℓ^(-1)).

Using the scaling δT∼(ℓ^(-(Subscript[λ, T]))) in Eq. (), one finds

ξ∼(((δT))^(1/(Subscript[λ, T]))).

Comparing with the definition of ν in Eq. (), ξ∼(((δT))^(-ν)), one concludes

ν=-1/(Subscript[λ, T]).

The exponent η concerns how fast the spin correlation decays with distance (larger η → more rapid decay).

Using the fluctuation-response relation, the variance of magnetization M=∫(d^d)r s(r) and the spin susceptibility χ=-(1/V)(∂_h^2)F are related

var M=-T(((∂^2)F)/(∂(h^2)))=V T χ.

Given that var M:=⟨(M^2)⟩-((⟨M⟩)^2), the variance can be expressed as an integration over the correlation function G(r),

var M=V∫(d^d)r G(r).

Derive Eq. ().

M=∫(d^d)r s(r),

var M=⟨(M^2)⟩-((⟨M⟩)^2)
=⟨(((∫(d^d)r s(r)))^2)⟩-((⟨∫(d^d)r s(r)⟩)^2)
=∫(d^d)r ∫(d^d)(r )'(⟨s(r)s(r')⟩-⟨s(r)⟩⟨s(r')⟩).

G(r-r')=⟨s(r)s(r')⟩-⟨s(r)⟩⟨s(r')⟩.

var M=∫(d^d)r ∫(d^d)(r )'G(r-r')
=∫(d^d)(r )'∫(d^d)r G(r)
=V∫(d^d)r G(r).

Using the form of the correlation function in Eq. (), G(r)∼((|r|)^(-d+2-η))(e^(-|r|/ξ)), the integration scales as

var M=V∫(d^d)r ((|r|)^(-d+2-η))(e^(-|r|/ξ))∝V (ξ^(2-η)).

Show the scaling of the integration in Eq. () by nondimentionalization.

var M=V∫(d^d)r ((|r|)^(-d+2-η))(e^(-|r|/ξ))
=V(∫_0^∞)dr (Subscript[A, d])(r^(d-1))(r^(-d+2-η))(e^(-r/ξ))
=V (Subscript[A, d])(∫_0^∞)dr (r^(1-η))(e^(-r/ξ)),

x=r/ξ,

var M=V (Subscript[A, d])(ξ^(2-η))(∫_0^∞)dx (x^(1-η))(e^(-x)).

Integrate[x^(1-η)Exp[-x],{x,0,∞},Assumptions->η<2]

Comparing Eq. () and Eq. (), one arrives at

χ∼(T^( -1)) (ξ^(2-η)).

At the critical point, T=(Subscript[T, c]) is a constant, and by Eq. (),

χ∼δ(T^( -(d+2(Subscript[λ, h]))/(Subscript[λ, T])))∼(ξ^( -(d+2(Subscript[λ, h])))),

where the scaling δT∼(ξ^(Subscript[λ, T])) in Eq. () was used. Comparing Eq. () and Eq. (),

η=2+d+2(Subscript[λ, h]).

Just by computing the behavior of the correlation function G(r) near the critical point, (Subscript[λ, T]) and (Subscript[λ, h]) can be determined.

Josephson’s relation []

Stell’s relation []

Fisher’s relation []

({{ ,  , α, β, γ, δ, ν, η}, {d=2, n=1, 0, 1/8, 7/4, 15, 1, 1/4}, {d=3, n=0, 0.24, 0.30, 1.16, 4.83, 0.59, 0.028}, {, n=1, 0.11, 0.33, 1.24, 4.81, 0.63, 0.034}, {, n=2, 0.00, 0.35, 1.32, 4.79, 0.67, 0.035}, {, n=3, -0.12, 0.37, 1.39, 4.79, 0.71, 0.036}, {d≥4,  , 0, 1/2, 1, 3, 1/2, 0}})

Each distinct set of critical exponents specifies a universality class of critical systems (or continuous phase transitions). Systems within a universality class share the same critical exponents and universal scaling functions.

Universality classes are affected by spatial dimensions d and symmetry groups (such as O(n)).

For O(n) models, above the upper critical dimension (d≥4), all models are governed by the mean-field universality class, as fluctuations becomes less significant in higher dimensions.

There is a special class of scaling relations, often referred to as hyperscaling relations, in which the spatial dimension d appears explicitly, such as  Eq. () and Eq. (). These relations break down for d>4, where the mean-field exponents prevails.

B. D. Josephson, Phys. Lett. 21, 608 (1966).

G. Stell, Phys. Rev. Lett. 20, 533 (1968).

M. E. Fisher, J. Math. Phys. 5, 944 (1964).

J. Zinn-Justin, Scholarpedia, 5(5):8346 (2010).

Critical Exponent Estimations

Critical Exponent ν

(w^((k+1)))=ℛ((w^((k)))),
ℛ(w):=(1/2)log cosh(2 (OverTilde[w])(w)).

Introduce the perturbation of the coupling constant w near the critical point (Subscript[w, c]),

to linearize the RG flow equation

Derive Eq. () based on Eq. ().

(w^((k+1)))=ℛ((w^((k)))),

(Subscript[w, c])+δ(w^((k+1)))=ℛ((Subscript[w, c])+δ(w^((k))))
=ℛ((Subscript[w, c]))+(∂_w)ℛ((Subscript[w, c]))δ(w^((k)))+…
=(Subscript[w, c])+(∂_w)ℛ((Subscript[w, c]))δ(w^((k)))+…

In each RG step, the perturbation δw will be multiplied by a factor (∂_w)ℛ((Subscript[w, c])) (the slop of ℛ(w) at w=(Subscript[w, c])). The fact that (∂_w)ℛ((Subscript[w, c]))>1 makes δw a relevant perturbation.

Block[{s,ptgt,pmdl,L,rec,wfun,wc},
s=Tuples[{-1,1},4];
ptgt[w_]:=#/Total[#]&@Exp[w #.RotateLeft[#]&/@s];
pmdl[wt_]:=#/Total[#]&@Sum[Exp[wt Total[#]sp&/@s],{sp,{-1,1}}];
L[w_,wt_]:=ptgt[w].Log[ptgt[w]/pmdl[wt]];
rec=Table[<|"w"->w,"L"->#1,"wt"->#2[[1,2]],"w1"->Log[Cosh[2 #2[[1,2]]]]/2|>&@@FindMinimum[L[w,wt],{wt,Sqrt[w],Sqrt[w]+w}],{w,0.01,2.5,0.01}];
wfun=Interpolation[Prepend[{#["w"],#["w1"]}&/@rec,{0,0}]];
wc=w/.FindRoot[wfun[w]==w,{w,0.4}];
Plot[{100(wfun[wc+0.01 x]-wc),x},{x,-1,1},PlotStyle->{Black,Directive[AbsoluteThickness[1],Dashed,Black]},Frame->True,FrameLabel->{"\!\( δw\^(k)\)","\!\(δw\^(k+1)\)"},PlotRange->{{-1,1},{-1.5,1.5}},Epilog->{AbsoluteThickness[1],Arrowheads[0.06],Red,Arrow[Flatten[{{#1,#2},{#2,#2}}&@@@Thread[{Most[#],Rest[#]}]&@NestList[100(wfun[wc+0.01 #]-wc)&,0.1,7],1]],
Blue,Arrow[Flatten[{{#1,#2},{#2,#2}}&@@@Thread[{Most[#],Rest[#]}]&@NestList[100(wfun[wc+0.01 #]-wc)&,-0.1,7],1]],Black,Text[StringTemplate["\!\(∂\_w ℛ(w\_c)≈`1`\)"][NumberForm[#,3]],{0,0},{-0.7,-1},{1,#}]&[wfun'[wc]]},AspectRatio->Automatic,ImageSize->160]]

Following Eq. (), starting with an initial perturbation δ(w^((0))), after k steps of iteration, the perturbation grows to

δ(w^((k)))=((((∂_w)ℛ((Subscript[w, c]))))^k)δ(w^((0))).

Meanwhile, the correlation length ξ shrinks to

(ξ^((k)))=(ξ^((0)))/((((Sqrt[2])))^k),

since the lattice scale ℓ enlarges by (Sqrt[2]) in each RG step, and the correlation length ξ∼(ℓ^(-1)) scales inversely with ℓ.

Eliminate k from Eq. () and Eq. (), one obtains

ξ∝δ(w^(-ν)),

Eliminate k and show Eq. ().

log(((δ(w^((k))))/(δ(w^((0))))))=k log (∂_w)ℛ((Subscript[w, c])).

log(((ξ^((k)))/(ξ^((0)))))=-k log (Sqrt[2]).

log(((ξ^((k)))/(ξ^((0)))))/log(((δ(w^((k))))/(δ(w^((0))))))=-((log (Sqrt[2]))/(log (∂_w)ℛ((Subscript[w, c]))))=-ν,

log(((ξ^((k)))/(ξ^((0)))))=-ν log(((δ(w^((k))))/(δ(w^((0)))))).

((ξ^((k)))/(ξ^((0))))=(((((δ(w^((k))))/(δ(w^((0)))))))^(-ν)).

with the exponent ν given by

Critical Exponent η

G(r)∼(1/((|r|)^(d-2+η))).

the correlation function G has been redefined on the lattice, indexed by two sites i and j

(Subscript[G, i j])=⟨(Subscript[s, i])(Subscript[s, j])⟩-⟨(Subscript[s, i])⟩⟨(Subscript[s, j])⟩
=(∑_s)(Subscript[s, i])(Subscript[s, j])p(s)-((∑_s)(Subscript[s, i])p(s))((∑_s)(Subscript[s, j])p(s)).

(Subscript[r, i j]):=|(Subscript[r, i])-(Subscript[r, j])| denotes the distance between sites i and j on a lattice, in unit of the lattice spacing.

Monte Carlo sampling: unbiased, but suffers from critical slowdown (i.e. local update becomes slow to equilibrate the system at criticality).

Inverse RG sampling: biased (approximate), but permits fast direct sampling at the RG fixed point.

Key idea: starting from a single spin, progressively fine-graining, regenerate large spin configurations in logarithmic steps.

Block[{bern,upsample},
bern[h_]:=Sign[Tanh[h]-RandomReal[{-1,1},Dimensions[h]]];
upsample[x_,w_]:=Block[{zs,xs,xraw},
zs=bern[w(x+RotateLeft[x,#])&/@{{0,-1},{-1,0},{0,1},{1,0}}];
xs=bern[w(zs+RotateLeft[zs])];
Flatten[{{xs[[1]],xs[[2]]},{xs[[4]],xs[[3]]}},{{3,1},{4,2}}]];
SeedRandom[42];
Row[Block[{nmax=7},
Framed[Image[#,ImageSize->Length@#],FrameMargins->None]&/@NestList[upsample[#,0.89]&,{{1}},nmax]],Style["→","Text"]]
]

The fine-graining is realized by sampling from the conditional distribution

p(s|s')=((e^(-(OverTilde[S])(s,s')))/((∑_s)(e^(-(OverTilde[S])(s,s'))))),
(OverTilde[S])(s,s')=-(Subscript[(OverTilde[w]), c])(∑_I)(∑_(i∈I))(Subscript[s, i])(Subscript[s, I]^′).

s - fine-grained spin configuration.

s' - coarse-grained spin configuration.

(Subscript[(OverTilde[w]), c]) - fixed the coupling at the critical value (where the variance of magnetization is maximal).

The sampling is efficient because given s', the components of s are independent, and can be sampled independently (no equilibrium process required).

Block[{bern,upsample,corr},
bern[h_]:=Sign[Tanh[h]-RandomReal[{-1,1},Dimensions[h]]];
upsample[x_,w_]:=Block[{zs,xs,xraw},
zs=bern[w(x+RotateLeft[x,#])&/@{{0,-1},{-1,0},{0,1},{1,0}}];
xs=bern[w(zs+RotateLeft[zs])];
Flatten[{{xs[[1]],xs[[2]]},{xs[[4]],xs[[3]]}},{{3,1},{4,2}}]];corr[x_]:=Block[{d=Length@x},ListCorrelate[x,x,{d/2,d/2}]/d^2];
Block[{n=7,M=1000,xs,G,d,G1D,dat,logr,fit,η},
SeedRandom[0];
xs=N@Table[Nest[upsample[#,0.89]&,{{RandomChoice[{-1,1}]}},n],M];
G=Mean[corr/@xs]-corr@Mean[xs];
d=Length@G;
G1D=KeySort@KeyMap[N,Mean/@GroupBy[{Norm[#1-{d/2,d/2}],#2}&@@@Most@ArrayRules[G],First->Last]];
dat=Table[{r,G1D[r]},{r,Cases[#,r_/;0<r<Max[#]/3]&@Keys[G1D]}];
fit=Fit[Log[dat],{1,logr},logr];
η=-Coefficient[fit,logr];
ListLogLogPlot[{Table[Exp@{logr,fit},{logr,Range[#1,#2,(#2-#1)/2]&@@MinMax@Log@dat[[All,1]]}],dat},Joined->{True,False},PlotStyle->{Red,Black},FrameLabel->{"\!\(r\_\(i j\)\)","\!\(G\_\(i j\)\)"},PlotLabel->StringTemplate["\!\(G\_\(i j\)∼r\_\(i j\)\%`1`\)"][NumberForm[-η,{3,2}]],ImageSize->200]
]]

Block[{bern,upsample,corr,expdata},
bern[h_]:=Sign[Tanh[h]-RandomReal[{-1,1},Dimensions[h]]];
upsample[x_,w_]:=Block[{zs,xs,xraw},
zs=bern[w(x+RotateLeft[x,#])&/@{{0,-1},{-1,0},{0,1},{1,0}}];
xs=bern[w(zs+RotateLeft[zs])];
Flatten[{{xs[[1]],xs[[2]]},{xs[[4]],xs[[3]]}},{{3,1},{4,2}}]];corr[x_]:=Block[{d=Length@x},ListCorrelate[x,x,{d/2,d/2}]/d^2];
SeedRandom[0];
expdata=Table[{n,Table[Block[{M=1000,xs,G,G1D,dat,logr,fit,η},
xs=N@Table[Nest[upsample[#,0.89]&,{{RandomChoice[{-1,1}]}},n],M];
G=Mean[corr/@xs]-corr@Mean[xs];
G1D=KeySort@KeyMap[N,Mean/@GroupBy[{Norm[#1-{2^(n-1),2^(n-1)}],#2}&@@@Most@ArrayRules[G],First->Last]];
dat=Table[{r,G1D[r]},{r,Cases[#,r_/;0<r<Max[#]/3]&@Keys[G1D]}];
fit=Fit[Log[dat],{1,logr},logr];
η=-Coefficient[fit,logr]],10]},{n,4,7}];
Block[{nmin,nmax},{nmin,nmax}=MinMax[expdata[[All,1]]];Graphics[{{Green,{Dashed,Line[{#,1/4}&/@{nmin-1,nmax+1}]},Text["exact",{nmin-1,1/4},{-1,-0.5}]},Translate[Block[{w=0.15,col=Red,ηs},ηs=Quantile[#2,{0,0.25,0.5,0.75,1}];
{col,Line[{0,#}&/@ηs[[{1,5}]]],{EdgeForm[col],Lighter[col,0.8],Rectangle[{-w,ηs[[2]]},{w,ηs[[4]]}]},Translate[Line[{0.8#,0}&/@{-w,w}],{0,#}]&/@ηs[[{1,3,5}]]}],{#1,0}]&@@@expdata},PlotRange->{All,{0.18,0.26}},PlotRangePadding->Scaled[0.05],AspectRatio->0.618,Frame->True,FrameLabel->{"\!\( L\)","\!\(η\)"},FrameTicks->{{Automatic,Automatic},{{#,StringTemplate["\!\(2\^`1`\)"][#]}&/@Range[nmin-1,nmax+1],None}},ImageSize->200]]
]