12 - Grand canonical potential

Determining μ((Subscript[n, 0]),T) for the 1D tight-binding chain

Clear[BosonicIntegral,FermionicIntegral];
BosonicIntegral[β_?NumberQ,μ_?NumberQ]:=(1/(2π))NIntegrate[(1/(Exp[β(-2Cos[k]-μ)]-1)),{k,-π,π},Method->{Automatic,"SymbolicProcessing"->False}];
FermionicIntegral[β_?NumberQ,μ_?NumberQ]:=(2/(2π))NIntegrate[(1/(Exp[β(-2Cos[k]-μ)]+1)),{k,-π,π},Method->{Automatic,"SymbolicProcessing"->False}];

FindRoot[FermionicIntegral[1/0.5,μ]==0.5,{μ,1.0}]

Clear[tabF];
n0range = {0.5,0.9,1.1,1.5};
Do[tabF[n0] = Table[{T,μ/.FindRoot[FermionicIntegral[1/T,μ]==n0,{μ,1.0}]},{T,0.01,1.,0.02}],{n0,n0range}]

Table[ListLinePlot[tabF[n0],FrameLabel->{"T/J","μ/J"},PlotLabel->"n0 = "<>ToString@n0],{n0,n0range}]

Clear[tabF2];
Trange = {0.05,0.1,0.2,0.5};
Do[tabF2[T] = Table[{n0,μ/.FindRoot[FermionicIntegral[1/T,μ]==n0,{μ,1.0}]},{n0,0.01,1.99,0.02}],{T,Trange}]

Table[ListLinePlot[tabF2[T],PlotRange->{{0,2},All},FrameLabel->{"(Subscript[n, 0])","μ/J"},PlotLabel->"T = "<>ToString@T],{T,Trange}]

Clear[tabB];
n0range = {0.5,1.0,1.5};
Do[tabB[n0] = Quiet@Table[{T,μ/.FindRoot[BosonicIntegral[1/T,μ]==n0,{μ,-2.01}]},{T,0.1,5,0.05}],{n0,n0range}]

Table[ListLinePlot[tabB[n0],PlotRange->{{0,All},{All,-1.9}},GridLinesStyle->Directive[Gray,Dashed],GridLines->{None,{-2}},FrameLabel->{"T/J","μ/J"},PlotLabel->"n0 = "<>ToString@n0],{n0,n0range}]

Clear[tabB2];
Trange = {0.5,1.0};
Do[tabB2[T] = Quiet@Table[{n0,μ/.FindRoot[BosonicIntegral[1/T,μ]==n0,{μ,-2.01}]},{n0,0.01,1.99,0.02}],{T,Trange}]

Table[ListLinePlot[tabB2[T],PlotRange->{{0,2},All},FrameLabel->{"(Subscript[n, 0])","μ/J"},PlotLabel->"T = "<>ToString@T],{T,Trange}]

Bose-Einstein condensation

Define the “fugacity”  z = (e^(β μ)). 
In terms of z we get

Integrate[((Sqrt[x])/((1/z)Exp[x]-1)),{x,0,∞}]

Mathematica has this integral implemented as the PolyLog function.
Now let’s try to solve the numerical equation

(1/2) (Sqrt[π]) PolyLog[(3/2),z] = ((Sqrt[π])/2)(N/V)(λ^3)

PolyLog[(3/2),z] = (N/V)(λ^3)

int[z_]:=((Sqrt[π])/2) PolyLog[(3/2),z];
Plot[int[z],{z,0,1}]

FindRoot[int[z] == 0.01,{z,0}]

The chemical potential is μ = (1/β)ln(z)

The right-hand side is proportional to (λ^(3 ))∝ (T^(-3/2)).
Let us choose units such that every other constant is one. 
We then solve

tab = Table[{T,T Log[z]/.FindRoot[int[z]==(T^(-3/2)),{z,0}]},{T,4,0.3,-0.05}]

ListLinePlot[tab,
Frame->True,
FrameLabel->{"T","μ"}]

Integrate[((Sqrt[x])/(Exp[x]-1)),{x,0,∞}]

N@%

Solve[((Sqrt[π])/2)n(((((2π (ℏ^2))/(m Tc))))^(3/2))==2.31516,Tc]

(1/2)(((((4(π^2))/( Gamma[3/2] Zeta[3/2]))))^(2/3))//N

(1/(2π (((2.612))^(2/3))))

1/%

x