140B Statistical Physics

Part 2. Quantum Gases

Thermal Quantum Gases

Bosons and Fermions

Quantum Many-Body System

Each person represents a particle,

Each role they take (be it a student, professor, staff, or service personnel) is akin to a quantum state.

Each role describes a person’s social functionality, just as each quantum state describes a particle’s physical properties, such as momentum, spin, charge, color, flavor, etc.

The occupation number (Subscript[n, α])=0,1,2,… represents the number of particles occupying the αth quantum state (denoted as (Ket[α])). (like the number of people in each role)

Each quantum state (Ket[α]) also has an associated energy (Subscript[ϵ, α]) -- the cost to keep a particle in the state. (like each role in a university comes with a salary -- the cost to keep a person in their position)

Single-Particle v.s. Many-Body State

Single-particle state: the state or mode of a particle (like the role of a person), denoted as (Ket[α]) (with α=1,2,3,… being the mode index).

Many-body state: the state of a system of particles (like the role of an institute), denoted as (Ket[n]) with

n:={(Subscript[n, α])∈ℕ|α=1,2,3…},

which is an integer vector that encodes occupation numbers on all available single-particle states (Ket[α]).

The total number of particles in the system (like the total number of people in our university) is

The total energy of the system (like the total salary expenditure of our university) is

[Note: to differentiate the two energies E and (Subscript[ϵ, α]), we often call E the many-body energy, and (Subscript[ϵ, α]) the single-particle energy.]

Identical Particles

For distinguishable particles, a two-particle system can be describe by the many-body quantum state (Ket[α])⊗(Ket[β]), where the 1st particle is in single-particle state (Ket[α]) and the 2nd in (Ket[β]).

For identical particles, many-body states like  (Ket[α])⊗(Ket[β]) and  (Ket[β])⊗(Ket[α]) are indistinguishable:

Style[Grid[{{"\!\(\[LeftBracketingBar]α⟩⊗\[LeftBracketingBar]β⟩\)","","\!\(\[LeftBracketingBar]β⟩⊗\[LeftBracketingBar]α⟩\)"},{"the 1st particle in \!\(\[LeftBracketingBar]α⟩\)\nthe 2nd particle in \!\(\[LeftBracketingBar]β⟩\)"," ","the 1st particle in \!\(\[LeftBracketingBar]β⟩\)\nthe 2nd particle in \!\(\[LeftBracketingBar]α⟩\)"},{"↘"," ","↙"},{"there is one particle in \!\(\[LeftBracketingBar]α⟩\), and another particle in \!\(\[LeftBracketingBar]β⟩\)",SpanFromLeft,SpanFromLeft},{"","\!\(\[LeftBracketingBar]1\_α 1\_β⟩\)",""}},Frame->{None,None,{{2,1}->True,{2,3}->True,{4,1}->True,{4,2}->True,{4,3}->True}},FrameStyle->LightGray,Alignment->{Center,Center}],"Text",LineSpacing->{0,14}]

The only legitimate description is (Ket[(Subscript[1, α]) (Subscript[1, β])]), as an equal-weight linear combination of  (Ket[α])⊗(Ket[β]) and  (Ket[β])⊗(Ket[α]). But there are still two possible ways to combine them:

Symmetric combination => bosons

(Subscript[(Ket[(Subscript[1, α]) (Subscript[1, β])]), B])=(1/(Sqrt[2]))((Ket[α])⊗(Ket[β])+(Ket[β])⊗(Ket[α])),

Anti-symmetric combination => fermions

(Subscript[(Ket[(Subscript[1, α]) (Subscript[1, β])]), F])=(1/(Sqrt[2]))((Ket[α])⊗(Ket[β])-(Ket[β])⊗(Ket[α])).

For bosons, (Subscript[(Ket[(Subscript[2, α])]), B]) is still a meaningful quantum state, with the state norm enhanced by (Sqrt[2]) due to a constructive interference:

(1/(Sqrt[2]))((Ket[α])⊗(Ket[α])+(Ket[α])⊗(Ket[α]))=(Sqrt[2])(Ket[α])⊗(Ket[α])=(Sqrt[2])(Subscript[(Ket[(Subscript[2, α])]), B]).

Boson enhancement factor: the probability to transfer a boson to a state that has already been occupied by n bosons will get enhanced by a factor of n+1.

=> Bosons like to stay in the same state (exhibiting an extroverted personality).

For fermions, the quantum state vanishes due to destructive interference.

(1/(Sqrt[2]))((Ket[α])⊗(Ket[α])-(Ket[α])⊗(Ket[α]))=0.

Pauli exclusion principle: two (or more) fermions can not occupy the same state simultaneously.

=> Fermions hate to stay in the same state (exhibiting an introverted personality).

Examples of Bosons and Fermions

Elementary particles (the Standard Model).

Block[{l=0.45,r=0.18,δ=0.14,col,pannel,grid},
col=<|"q"->Purple,"l"->Cyan,"g"->Orange,"h"->Red|>;
pannel[symbol_,name_,color_,mass_,charge_,spin_]:=Block[{dl},dl=StringMatchQ[name,__~~"\n"~~__];{{EdgeForm[color],Opacity[0],Rectangle[{-l,-l},{l,l},RoundingRadius->0.1]},Translate[{{color,Disk[{0,0},r]},{White,Text[Style[symbol,SingleLetterItalics->False,FontFamily->"CMU Sans Serif"],{0,0}]}},{0,If[dl,0.1,0]}],{Black,Text[Style[name,If[dl,11,12],LineSpacing->{0,If[dl,10,11]},FontFamily->"CMU Sans Serif"],{0,-l},{0,If[dl,-0.7,-1]}]},Text[Style[#1,8,FontFamily->"CMU Bright"],{-0.9l,l-δ#2},{-1,0}]&@@@Thread[{{mass,charge,spin},{1,2,3}}]}];
grid[pans_]:=Translate[Translate[#1,{#2,0}]&@@@Thread[{#1,Range[0,Length[#1]-1]}],{0,-#2}]&@@@Thread[{pans,Range[0,Length[pans]-1]}];
Graphics[{{{{Lighter[Blue,0.9],Rectangle[{-0.95,-3.6},{2.6,1.3},RoundingRadius->0.2]},Text[Style["Fermions\n(matter constituent)",LineSpacing->{0,12}],{1,1}]},Text[Style[#1,12,FontFamily->"CMU Concrete"],{#2,l},{0,-0.8}]&@@@Thread@{{"I","II","III"},{0,1,2}},
grid@{{pannel["u","up",col@"q","2.2MeV","2/3","1/2"],pannel["c","charm",col@"q","1.28GeV","2/3","1/2"],pannel["t","top",col@"q","173GeV","2/3","1/2"]},{pannel["d","down",col@"q","4.7MeV","-1/3","1/2"],pannel["s","strange",col@"q","96MeV","-1/3","1/2"],pannel["b","bottom",col@"q","4.18GeV","-1/3","1/2"]},{pannel["e","electron",col@"l","0.511MeV","-1","1/2"],pannel["μ","muon",col@"l","106MeV","-1","1/2"],pannel["τ","tauon",col@"l","1.78GeV","-1","1/2"]},{pannel["\!\(ν\_e\)","electron\nneutrino",col@"l","<1eV","0","1/2"],pannel["\!\(ν\_μ\)","muon\nneutrino",col@"l","<0.17MeV","0","1/2"],pannel["\!\(ν\_τ\)","tauon\nneutrino",col@"l","<18.2MeV","0","1/2"]}},
{col@"q",Text[Style["QUARKS",Bold,FontFamily->"CMU Sans Serif Demi Condensed"],{-l,-1-l},{-1.0,-0.8},{0,1}]},
{col@"l",Text[Style["LEPTONS",Bold,FontFamily->"CMU Sans Serif Demi Condensed"],{-l,-3-l},{-1.0,-0.8},{0,1}]},
Translate[{Gray,Text[Style[#1,8,FontFamily->"CMU Bright"],{-0.05,-δ#2},{1,0}]&@@@Thread[{{"mass","charge","spin"},{1,2,3}}]},{-l,l}]},
Translate[{{{Lighter[Brown,0.9],Rectangle[{-0.6,-3.6},{1.6,1.3},RoundingRadius->0.2]},Text[Style["Bosons\n(force mediator)",LineSpacing->{0,12}],{0.5,1}]},grid@{{pannel["g","gluon",col@"g","0","0","1"],pannel["H","higgs",col@"h","125GeV","0","0"]},{pannel["γ","photon",col@"g","0","0","1"]},{pannel["W","W boson",col@"g","80.4GeV","±1","1"]},{pannel["Z","Z boson",col@"g","91.2GeV","0","1"]}},
{col@"g",Text[Style["GAUGE BOSONS",Bold,FontFamily->"CMU Sans Serif Demi Condensed"],{l,-3-l},{-1,1},{0,1}],Text[Style["VECTOR BOSONS",10,Bold,FontFamily->"CMU Sans Serif Demi Condensed"],{l,-3-l},{-1.02,3},{0,1}]},
{col@"h",Text[Style["SCALAR BOSONS",Bold,FontFamily->"CMU Sans Serif Demi Condensed"],{1+l,-l},{1.05,-0.4},{0,1}]}},{3.3,0}]},PlotRangePadding->0,ImageSize->320]]

Fermions are of half-integer spins.

Bosons are of integer spins.

Inside an atom ( isotope):

Block[{a1=0.14,a2=0.24,a3=0.04,b3=0.2,re={1.2,0.7},rγ={1.2,0.3},rc={1.2,-0.4},rp={0.8,-0.3},col},
col=<|"e"->Cyan,"p"->Blue,"n"->Green,"γ"->Orange|>;
Graphics[{{Text[Style["atom",FontFamily->"CMU Sans Serif"],{0,0.8}]},
{col@"e",Opacity[0.02],Disk[{0,0},#^2]&/@Range[1,0.02,-0.02]},{Block[{f},f=BSplineFunction[{{0.34,0.3},{0.7,0.7},re}];{col@"e",Line[f/@Range[0,0.8,0.05]]}],Translate[{{Cyan,Disk[{0,0},a1]},{White,Text[Style["e",SingleLetterItalics->False,FontFamily->"CMU Sans Serif"],{0,0}]},{col@"e",Text[Style["electron (fermion)",FontFamily->"CMU Sans Serif"],{a1,0},{-1.1,0}]}},re]},
{Black,Point[{0,0}],Block[{f},f=BSplineFunction[{{0,0},{0.3,-0.4},rc}];Line[f/@Range[0.1,0.75,0.05]]]},
Translate[{col@"γ",Line[{a3 Sin[4π #]Cos[π #/2],b3 #}&/@Range[-1,1,0.01]],Text[Style["photon (boson)",FontFamily->"CMU Sans Serif"],{a3,0},{-1.24,-0.5}]},rγ],Translate[{{col@#1,Disk[#2,a1],{White,Text[Style[#1,FontFamily->"CMU Sans Serif"],#2]}}&@@@Thread@{Append[Flatten@ConstantArray[{"n","p"},3],"n"],Append[CirclePoints[{a2,π/6},6],{0,0}]},Text[Style["necleus",FontFamily->"CMU Sans Serif"],{0,-a1-a2},{0,1}],
Block[{f},f=BSplineFunction[{a2 ReIm@Exp[-I π/6],rp}];{{col@"p",Line[f/@Range[0.3,0.5,0.1]]},{col@"n",Line[{1,-1}#&/@f/@Range[0.3,0.5,0.1]]}}],
Translate[{col@"p",Disk[{0,0},a2],{White,Text[Style[#1,FontFamily->"CMU Sans Serif"],#2]&@@@Thread@{{"u","u","d"},CirclePoints[{a1,π/6},3]}},{col@"p",Text[Style["proton (fermion)",FontFamily->"CMU Sans Serif"],{a2,0},{-1.1,0}]}},rp],
Translate[{col@"n",Disk[{0,0},a2],{White,Text[Style[#1,FontFamily->"CMU Sans Serif"],#2]&@@@Thread@{{"d","d","u"},CirclePoints[{a1,π/6},3]}},{col@"n",Text[Style["neutron (fermion)",FontFamily->"CMU Sans Serif"],{a2,0},{-1.1,0}]}},rp{1,-1}]},rc]},ImageSize->320]]

Rules (due to the multiplication rule of fermion sign):

boson+boson=boson,
fermion+boson =fermion,
fermion+fermion=boson.

Total number of fermions in an atom

(Subscript[N, F])=(Subscript[N, p])+(Subscript[N, n])+(Subscript[N, e])(Overscript[=, neutral])2(Subscript[N, p])+(Subscript[N, n]).

A neutral atom is a fermion (boson) iff the number of neutrons (Subscript[N, n]) is odd (even).

Statistical Distributions

What is Ideal Gas?

Its energy is a simple sum of energies of particles in each single-particle mode

E(n)=(∑_α)(Subscript[ϵ, α])(Subscript[n, α]).

The energy function will contain higher-order interaction terms

E(n)=(∑_α)(Subscript[ϵ, α])(Subscript[n, α])+(Underscript[(Underscript[(∑_(α β))(Subscript[ϵ, α β])(Subscript[n, α])(Subscript[n, β]), ⏟]), two-body interaction])+(Underscript[(Underscript[(∑_(α β γ))(Subscript[ϵ, α β γ])(Subscript[n, α])(Subscript[n, β])(Subscript[n, γ]), ⏟]), three-body interaction])+…

The interaction, say (Subscript[ϵ, α β]), is

Attractive, if (Subscript[ϵ, α β])<0: energy is reduced when particles attract each other. (like people in a cooperative environment will reduce operational cost)

Repulsive, if (Subscript[ϵ, α β])>0: energy is increased when particles repel each other. (like people in a confrontational environment will increase operational cost)

If the interaction is weak (Subscript[ϵ, α β]),(Subscript[ϵ, α β γ]),…≪(Subscript[ϵ, α]), the system can remain in the gas phase. If the interaction gets strong enough, the system may transition into liquid or solid phases.

Grand Partition Function

Each microstate (Ket[n]) of a quantum many-body system is labeled by

n:={(Subscript[n, α])|α=1,2,…}.

Following Eq. () and Eq. (), the energy and particle number functions are

E(n)=(∑_α)(Subscript[ϵ, α])(Subscript[n, α]), N(n)=(∑_α)(Subscript[n, α]).

The probability p(n) to observe the system in the microstate (Ket[n]) is

p(n)=(1/ℨ)(e^(-β(E(n)-μ N(n)))),

where grand partition function ℨ is given by

ℨ=(∑_n)(e^(-β(E(n)-μ N(n)))).

The probability to observe (Subscript[n, α]) particles in the single-particle state (Ket[α]) is

With the single-mode grand partition function

ℨ=(∏_α)(Subscript[ℨ, α]),

Derive Eq. ( - ) from Eq. (, ).

ℨ=(∑_n)(e^(-β(E(n)-μ N(n))))
=(∑_n) exp(-β((∑_α)(Subscript[ϵ, α])(Subscript[n, α])-μ(∑_α)(Subscript[n, α])))
=(∑_n)exp(-β(∑_α)((Subscript[ϵ, α])(Subscript[n, α])-μ)(Subscript[n, α]))
=(∑_n)(∏_α)(e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α])))
= (∏_α)(∑_(Subscript[n, α]))(e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α]))).

(∑_n)(∏_α)f((Subscript[n, α]))=(∑_(Subscript[n, 1]))(∑_(Subscript[n, 2]))(∑_(Subscript[n, 3]))… (f((Subscript[n, 1]))f((Subscript[n, 2]))f((Subscript[n, 3]))…)
=(∑_(Subscript[n, 1])) f((Subscript[n, 1]))(∑_(Subscript[n, 2]))f((Subscript[n, 2]))(∑_(Subscript[n, 3]))f((Subscript[n, 3]))…
=(∏_α)(∑_(Subscript[n, α]))f((Subscript[n, α])).

(Subscript[ℨ, α])=(∑_(Subscript[n, α]))(e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α]))),

p(n)=(1/ℨ)(e^(-β(E(n)-μ N(n))))
=(1/ℨ)exp(-β((∑_α)(Subscript[ϵ, α])(Subscript[n, α])-μ(∑_α)(Subscript[n, α])))
=(1/ℨ)exp(-β(∑_α)((Subscript[ϵ, α])(Subscript[n, α])-μ)(Subscript[n, α]))
=(1/ℨ)(∏_α)(e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α])))
=(∏_α)((e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α])))/(Subscript[ℨ, α])).

p((Subscript[n, α]))=(1/(Subscript[ℨ, α]))(e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α]))),

p(n)=(∏_α)p((Subscript[n, α])),

log ℨ=(∑_α)log (Subscript[ℨ, α]).

Bosons

(Subscript[ℨ, α])=(∑_((Subscript[n, α])=0)^∞)(e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α])))=(1/(1-(e^(-β ((Subscript[ϵ, α])-μ))))).

Evaluate the summation in Eq. ().

Sum[Exp[-β ξ n],{n,0,∞}]

Fermions

(Subscript[ℨ, α])=(∑_((Subscript[n, α])=0,1))(e^(-β((Subscript[ϵ, α])-μ)(Subscript[n, α])))=1+(e^(-β ((Subscript[ϵ, α])-μ))).

Bose-Einstein Distribution

p((Subscript[n, α]))=(1-(e^(-β  ((Subscript[ϵ, α])-μ))))(e^(-β  ((Subscript[ϵ, α])-μ)(Subscript[n, α]))).

The distribution is tuned by a joint parameter  β((Subscript[ϵ, α])-μ).

Manipulate[DiscretePlot[Exp[- βϵ n](1-Exp[-βϵ]),{n,0,10},PlotStyle->Black,PlotRange->{0,1},PlotRangePadding->Scaled[.05],PlotLabel->StringTemplate[" \!\(β (ϵ\_α-μ)=`1`\)"][βϵ],FrameLabel->{"\!\(n\_α\)","\!\(p(n\_α)\)"},Prolog->{Translate[{Blue,Line[{0,#}&/@{0,1}],Text[StringTemplate["\!\(⟨n\_α⟩=`1`\)"][NumberForm[1/(Exp[βϵ]-1),{2,2}]],{0,1},{-1.1,1}]},{1/(Exp[βϵ]-1),0}]},ImageSize->180],{{βϵ,1," \!\(β (ϵ\_α-μ)\)"},0.01,5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Evaluate the summation in Eq. ().

Sum[n(1-Exp[-β ξ])Exp[-β ξ n],{n,0,∞}]

Manipulate[Block[{β=1/T},Plot[1/(Exp[β ξ]-1),{ξ,0,5},PlotStyle->Blue,PlotRange->{0,5},Filling->Axis,PlotRangePadding->Scaled[.05],
PlotLabel->StringTemplate[" \!\(T =`1`\)"][T],FrameLabel->{"\!\(ϵ\_α-μ\)","\!\(⟨n\_α⟩\)"},ImageSize->180]],{{T,1," \!\(T=1/β\)"},0.01,5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

The Bose-Einstein distribution is well-defined only for (Subscript[ϵ, α])-μ>0.

Fermi-Dirac Distribution

p((Subscript[n, α]))=((e^(-β ((Subscript[ϵ, α])-μ)(Subscript[n, α])))/(1+(e^(-β ((Subscript[ϵ, α])-μ)))))=({{{, ({{(1/((e^(-β ((Subscript[ϵ, α])-μ)))+1)), (Subscript[n, α])=0,}, {(1/((e^( β ((Subscript[ϵ, α])-μ)))+1)), (Subscript[n, α])=1.}})}})

The distribution is tuned by a joint parameter  β((Subscript[ϵ, α])-μ).

Manipulate[DiscretePlot[Exp[- βϵ n]/(1+Exp[-βϵ]),{n,0,1},PlotStyle->Black,PlotRange->{0,1},PlotRangePadding->Scaled[.05],PlotLabel->StringTemplate[" \!\(β (ϵ\_α-μ)=`1`\)"][βϵ],FrameLabel->{"\!\(n\_α\)","\!\(p(n\_α)\)"},Prolog->{Translate[{Red,Line[{0,#}&/@{0,1}],Text[StringTemplate["\!\(⟨n\_α⟩=`1`\)"][NumberForm[1/(Exp[βϵ]+1),{2,2}]],{0,1},{-1.1 Sign[βϵ],1}]},{1/(Exp[βϵ]+1),0}]},ImageSize->180],{{βϵ,0," \!\(β (ϵ\_α-μ)\)"},-5,5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Manipulate[Block[{β=1/T},Plot[1/(Exp[β ξ]+1),{ξ,-5,5},PlotStyle->Red,PlotRange->{0,1},Filling->Axis,PlotRangePadding->Scaled[.05],
PlotLabel->StringTemplate[" \!\(T =`1`\)"][T],FrameLabel->{"\!\(ϵ\_α-μ\)","\!\(⟨n\_α⟩\)"},ImageSize->180]],{{T,1," \!\(T=1/β\)"},0.01,5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

The Fermi-Dirac distribution is applicable for any (Subscript[ϵ, α])-μ∈ℝ.

Maxwell-Boltzmann Distribution

η=({{{, ({{+1, for bosons,}, {-1, for fermions.}})}})

For bosons, the chemical potential  μ can not exceed the lowest energy level (i.e. ∀α:μ<(Subscript[ϵ, α])).

For fermions, there is not restriction on the chemical potential.

In the dilute limit, the average occupancy of each single-particle state is close to zero (compare to 1), i.e. ⟨(Subscript[n, α])⟩→0 (such that N=(∑_α)⟨(Subscript[n, α])⟩ is also small).

Based on Eq. (), this limit is achieved when (e^( β ((Subscript[ϵ, α])-μ)))→∞, such that the fermion sign η in the denominator is negligible, and the average occupation number

will be close to zero.

Block[{β=1},LogPlot[{1/(Exp[β ξ]-1),1/(Exp[β ξ]+1),Exp[-β ξ]},{ξ,0,7},PlotStyle->{Blue,Red,Directive[Dashed,Black]},GridLines->{None,Automatic},PlotRange->{0.001,10},PlotRangePadding->Scaled[.05],
FrameLabel->{"\!\(β(ϵ\_α-μ)\)","\!\(⟨n\_α⟩\)"},PlotLegends->Placed[Style[#,"Text"]&/@{"Bose-Einstein","Fermi-Dirac","Maxwell-Boltzmann"},Right],AspectRatio->1,ImageSize->180]]

Energy and Particle Number

The average total energy (defined as E:=(∑_n)E(n)p(n)):

The average total particle number (defined as N:=(∑_n)N(n)p(n)):

Continuum Limit

Particles in Free Space

α=(k,σ),

Momentum k=((Subscript[k, 1]),(Subscript[k, 2]),…,(Subscript[k, D])) is a D-dimensional vector for particles in D-dimensional space.

Block[{ah,k=2π{3,4},ψ,A=0.03},
ah=Graphics[Line[Offset/@{{-4,-2},{0,0},{-4,2}}]];
ψ[{x_,y_}]:=Re@Exp[I k.{x,y}];
Graphics[{First@ContourPlot[ψ[{x,y}],{x,0,1},{y,0,1},ContourStyle->None,Contours->Range[-1.1,1.1,0.2],ColorFunction->(Lighter[Blend[{Cyan,Gray,Yellow},#],0.8]&)],
Block[{y=0.75/4},{Red,{AbsoluteThickness[1],Opacity[0.5],Line[{#,y}&/@{0,1}]},Line[Table[{x,y+A ψ[{x,y}]},{x,0,1,0.01}]]}],Block[{x=0.75/3},{Green,{AbsoluteThickness[1],Opacity[0.5],Line[{x,#}&/@{0,1}]},Line[Table[{x-A ψ[{x,y}],y},{y,0,1,0.01}]]}],
{Black,Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}]},{AbsoluteThickness[1],Arrowheads[{{-0.01,0,ah},{0.01,1,ah}}],GeometricTransformation[Translate[{Arrow[{{0,0},{1,0}}],Scale[Text[Style[" \!\( L\) ",Background->White],{0.5,0},{0,0.2}],{1,(-1)^#},{0,0}]},{0,-0.1}],MatrixPower[{{0,1},{1,0}},#]]&/@{0,1}},
Translate[{Arrowheads[0.08],Arrow[{{0,0},0.01k}],Text["k",0.01k,{-1,0}]},{0.5,0.5}],Text["(\!\(D=2\))",{0.5,1},{0,-1}]},ImageSize->120]]

In quantum mechanics, particles are waves, and the momentum (Subscript[k, i]) is set by the wave length (Subscript[λ, i]) via the De Broglie relation in the corresponding direction,

(Subscript[k, i])=((2π)/(Subscript[λ, i]))=((2π (Subscript[m, i]))/L)    ((Subscript[m, i])∈ℤ),

which must be quantized in a box of finite size L.

The volume of the system (the box) will be

V=(L^D).

Spin σ labels the internal degrees of freedom of the particle. The number of spin states is called the spin degeneracy, denoted as (Subscript[g, s]).

Electron (spin-1/2 fermion):

(Subscript[g, s])=2.

Two spin states: σ=±1/2.

Photon (spin-1 gauge boson):

(Subscript[g, s])=D-1.

(D-1) spin (polarization) states: σ=±1 for D=3. The σ=0 photon (longitudinal) photon is forbidden by the gauge invariance.

Phonon (spin-1 scalar boson):

(Subscript[g, s])=D.

D spin (polarization) states: σ=0,±1 for D=3.

Usually (without internal symmetry breaking), the energy (Subscript[ϵ, α])=(Subscript[ϵ, k,σ]) does not depend on the spin degrees of freedom σ, so different spin states are degenerated in energy.

Dispersion Relations

Non-relativistic particles (electrons in metal):

with m - the inertial mass of the particle.

Relativistic particles (photons, acoustic phonons, electrons in graphene):

with c - the limiting speed of the particle, e.g. the speed light/sound or fermi-velocity.

Density of States

E=V∫dϵ ((ϵ g(ϵ))/((e^( β (ϵ-μ)))-η)).

N=V∫dϵ ((g(ϵ))/((e^( β (ϵ-μ)))-η)).

General Principle

Block[{a=0.1,R=0.8,dR=0.05,pts,ah},
pts=Tuples[Range[-1,1,a],2];
ah=Graphics[Line[Offset/@{{-4,-2},{0,0},{-4,2}}]];
Graphics[{{Arrowheads[0.05],AbsoluteThickness[1],Arrow[{#,0}&/@{-1,1}],Arrow[{0,#}&/@{-1,1}],Text["\!\(k\_1\)",{1,0},{0,1}],Text["\!\(k\_2\)",{0,1},{1,0}]},
First@RegionPlot[RegionDifference[Disk[{0,0},R+dR],Disk[{0,0},R-dR]],PlotStyle->LightGray,BoundaryStyle->Directive[Dashed,AbsoluteThickness[1],Black]],
{AbsolutePointSize[4],Point@Cases[pts,p_/;R-dR<Norm[p]<R+dR],Gray,Point@Cases[pts,p_/;Norm[p]<=R-dR]},
Translate[{AbsoluteThickness[1],Gray,Line[{{-1,-1},{-1,1},{1,1},{1,-1},{-1,-1}}a/2]},a{2,2}],
Translate[{Gray,{AbsolutePointSize[8],Point[{0,0}]},{AbsoluteThickness[1],Line[{{-1,-1},{-1,1},{1,1},{1,-1},{-1,-1}}a],{Arrowheads[{{-0.01,0,ah},{0.01,1,ah}}],Black,Arrow[a{#,-1.5}&/@{-1,1}],Arrow[a{1.5,#}&/@{-1,1}],Text[Style["\!\(\(2π\)\/L\)","InlineFormula"],{0,-1.8a},{0,1}],Text[Style["\!\(\(2π\)\/L\)","InlineFormula"],{1.8a,0},{-1,0}]}}},{0.85,0.95}],
Text["(\!\(D=2\) case)",{0,-1},{0,1}]}]]

In this picture, each point corresponds to a state (Ket[k,σ]) in the momentum space with (Subscript[g, s])-fold spin degeneracy. The momentum is discretized according to Eq. ().

The continuum limit refers to the limit L→∞, such that the states are densely distributed, forming a continuum.

Consider a ball in the momentum space of radius k, its volume is

Ω=(∫_(|k|<k))(d^D)k 
=(∫_0^k)(Subscript[A, D])(k^(D-1))dk
=((Subscript[A, D])/D)(k^D),

where (Subscript[A, D]) is the area of a (D-1)-dimensional hypersphere of unit radius (in the D-dimensional space), given by

(Subscript[A, D])=((2(π^(D/2)))/(Γ(D/2))),

where Γ(n)=(∫_0^∞)(x^(n-1))(e^(-x))dx is the Euler Gamma function.

Style[Grid[Transpose@Prepend[Append[Table[{d,2π^(d/2)/Gamma[d/2]},{d,1,10}],{"…","…"}],{"\!\( D\)","\!\(A\_D\)"}],Dividers->{{2->True},{2->True}}],"DisplayFormula"]

The number of states (Subscript[Ν, k]) in the ball of radius k is (approximately)

(Subscript[Ν, k])=(∑_α)1=(Subscript[g, s])(∑_(|k|<k))1=(Subscript[g, s])(Ω/(((((2π)/L)))^D))=(((Subscript[g, s]) (Subscript[A, D])V)/((((2π))^D)D)) (k^D).

On the other hand, we expect the DOS function g(ϵ) to be such that

(Subscript[Ν, k])=(∑_α)1=V(∫_0^(Subscript[ϵ, k]))dϵ g(ϵ).

g(ϵ)=(1/V)(Subscript[(((∂(Subscript[Ν, k]))/(∂(Subscript[ϵ, k])))), (Subscript[ϵ, k])=ϵ])=(1/V)(Subscript[(((∂(Subscript[Ν, k]))/(∂k))((∂k)/(∂(Subscript[ϵ, k])))), (Subscript[ϵ, k])=ϵ])
=(((Subscript[g, s])(Subscript[A, D]))/(((2π))^D)) (Subscript[((k^(D-1))((∂k)/(∂(Subscript[ϵ, k])))), (Subscript[ϵ, k])=ϵ]),

Specific Examples

For non-relativistic particles (with the inertial mass m)

Derive Eq. () given the dispersion relation Eq. ().

(Subscript[ϵ, k])=((k^2)/(2m)),

k=(Sqrt[2m (Subscript[ϵ, k])])
=>((∂k)/(∂(Subscript[ϵ, k])))=(Sqrt[(m/(2(Subscript[ϵ, k])))]).

g(ϵ)=(((Subscript[g, s])(Subscript[A, D]))/(((2π))^D)) (Subscript[((k^(D-1))((∂k)/(∂(Subscript[ϵ, k])))), (Subscript[ϵ, k])=ϵ])
=(((Subscript[g, s])(Subscript[A, D]))/(((2π))^D)) (Subscript[(((((Sqrt[2m (Subscript[ϵ, k])])))^(D-1))(Sqrt[(m/(2(Subscript[ϵ, k])))])), (Subscript[ϵ, k])=ϵ])
=(((Subscript[g, s])(Subscript[A, D]))/(2(((2π))^D))) (((2m))^(D/2))(ϵ^(D/2-1)).

Block[{m=1},Column[Table[Plot[x^(d/2-1)UnitStep[x],{x,-1,3},Epilog->{Text[StringTemplate["\!\(D=`1`\)"][d],Scaled@{1,1},{1.2,1}]},PlotStyle->Black,Filling->0,FillingStyle->Opacity[0.15],PlotRange->{0,2.4},PlotRangePadding->Scaled[.03],Exclusions->None,FrameLabel->{"\!\( ϵ\)","\!\(g(ϵ)\)"},ImagePadding->{{40,2},{If[d<3,2,42],2}},ImageSize->160],{d,3}],Spacings->0]]

For relativistic particles (with the limiting speed c)

Derive Eq. () given the dispersion relation Eq. ().

(Subscript[ϵ, k])=c k,

k=((Subscript[ϵ, k])/c)
=>((∂k)/(∂(Subscript[ϵ, k])))=(1/c).

g(ϵ)=(((Subscript[g, s])(Subscript[A, D]))/(((2π))^D)) (Subscript[((((((Subscript[ϵ, k])/c)))^(D-1))(1/c)), (Subscript[ϵ, k])=ϵ])
=(((Subscript[g, s])(Subscript[A, D]))/(((2π))^D)) (c^(-D))(ϵ^(D-1)).

Block[{m=1},Column[Table[Plot[Abs[x]^(d-1)/d UnitStep[x],{x,-1,3},Epilog->{Text[StringTemplate["\!\(D=`1`\)"][d],Scaled@{1,1},{1.2,1}]},PlotStyle->Black,Filling->0,FillingStyle->Opacity[0.15],PlotRange->{0,2.4},PlotRangePadding->Scaled[.03],Exclusions->None,FrameLabel->{"\!\( ϵ\)","\!\(g(ϵ)\)"},ImagePadding->{{40,2},{If[d<3,2,42],2}},ImageSize->160],{d,3}],Spacings->0]]

the dimension of space D,

the dynamical exponent ζ (such that (Subscript[ϵ, k])∼(k^ζ)):

ζ=1 for relativistic particles with linear dispersion,

ζ=2 for non-relativistic particles with quadratic dispersion.

Grid[{{"System","D","ζ","Δ"},{"Electrons in metal","3","2","3/2"},{"Electrons in graphene","2","1","2"},{"Electrons in field\[Hyphen]effect transistors","2","2","1"},{"Electrons in nanotube","1","2","1/2"},{"Electrons on quantum Hall insulator boundary","1","1","1"},{"Phonons in solid","3","1","3"},{"Atomic gas","3","2","3/2"},{"Cold atoms in 2D optical trap","2","2","1"}},ItemStyle->"Text",Alignment->Left,ItemSize->{{15,Automatic},Automatic},Dividers->{None,2->True}]

g(ϵ)dϵ=d((((ϵ/(Subscript[ϵ, 0]))))^Δ),

(1/(Subscript[ϵ, 0]^Δ))=({{{, ({{((Subscript[g, s])/(Γ(Δ/2+1)))((((1/(2 (π^(1/2))c))))^Δ), ζ=1 (and Δ=D),}, {((Subscript[g, s])/(Γ(Δ+1)))((((m/(2π))))^Δ), ζ=2 (and Δ=D/2).}})}})

Derive Eq. ().

(Subscript[A, D])=((2(π^(D/2)))/(Γ(D/2))).

n Γ(n)=Γ(n+1),

For ζ=1, we have Δ=D, and

g(ϵ)=(((Subscript[g, s])(Subscript[A, D]))/(((2π))^D)) (c^(-D))(ϵ^(D-1))=((Δ (ϵ^(Δ-1)))/(Subscript[ϵ, 0]^Δ)),

so that

(1/(Subscript[ϵ, 0]^Δ))=(((Subscript[g, s])(Subscript[A, Δ]))/(Δ (((2π))^Δ))) (c^(-Δ))
=((Subscript[g, s])/(Δ (((2π))^Δ))) ((2(π^(Δ/2)))/(Γ(Δ/2)))(c^(-Δ))
=((Subscript[g, s])/((2^Δ)(π^Δ)))((2(π^(Δ/2)))/(2(Δ/2)Γ(Δ/2))) (1/(c^Δ))
=((Subscript[g, s])/((2^Δ)(π^(Δ/2))))(1/(Γ(Δ/2+1))) (1/(c^Δ))
=((Subscript[g, s])/(Γ(Δ/2+1))) ((((1/(2 (π^(1/2))c))))^Δ).

For ζ=2, we have Δ=D/2, and

g(ϵ)=(((Subscript[g, s])(Subscript[A, 2Δ]))/(2(((2π))^(2Δ)))) (((2m))^Δ)(ϵ^(Δ-1))=((Δ (ϵ^(Δ-1)))/(Subscript[ϵ, 0]^Δ)),

so that

((1 )/(Subscript[ϵ, 0]^Δ))=(((Subscript[g, s])(Subscript[A, 2Δ]))/(2Δ (((2π))^(2Δ)))) (((2m))^Δ)
=((Subscript[g, s])/(2Δ (((2π))^(2Δ))))((2(π^Δ))/(Γ(Δ))) (((2m))^Δ)
=((Subscript[g, s])/(2 (4^Δ)(π^(2Δ))))((2(π^Δ))/(Δ Γ(Δ))) (((2m))^Δ)
=((Subscript[g, s])/((4^Δ)(π^Δ)))(1/(Γ(Δ+1))) (((2m))^Δ)
=((Subscript[g, s])/(Γ(Δ+1))) (((((2m)/(4π))))^Δ)
=((Subscript[g, s])/(Γ(Δ+1))) ((((m/(2π))))^Δ).

Thermodynamic Properties

Free Energy

𝔉=-T log ℨ.

According to Eq. (), for ideal gases, the free energy arises from the independent contributions of each single-particle mode.

𝔉=-T(∑_α)log (Subscript[ℨ, α])
=η T V(∫_0^∞)dϵ g(ϵ)log(1-η (e^(-β (ϵ-μ)))).

For bosons: z∈(0,1) (as  μ<0).

For fermions: z∈(0,+∞).

η=±1 is the fermion sign (+1 for bosons, -1 for fermions).

s counts the total power of x (including dx) in the numerator.

Derive and verify Eq. (, ).

𝔉=η T V(∫_0^∞)dϵ g(ϵ)log(1-η (e^(-β (ϵ-μ))))
=η T V(∫_0^∞)d((((ϵ/(Subscript[ϵ, 0]))))^Δ)log(1-η (e^(β μ))(e^(-β ϵ))).

𝔉=η T V(∫_0^∞)d(((((β ϵ)/(β (Subscript[ϵ, 0])))))^Δ)log(1-η z (e^(-β ϵ)))
=((η T V)/(((β (Subscript[ϵ, 0])))^Δ))(∫_0^∞)d((((β ϵ))^Δ))log(1-η z (e^(-β ϵ)))
=((η (T^(Δ+1)) V)/(Subscript[ϵ, 0]^Δ))(∫_0^∞)d((x^Δ))log(1-η z (e^(-x))).

(∫_0^∞)d((x^Δ))log(1-η z (e^(-x)))=(Subscript[[(x^Δ)log(1-η z (e^(-x)))], x=0]^(x=∞))-(∫_0^∞)dx((η (x^Δ))/((z^(-1))(e^x)-η)).

𝔉=-(((T^(Δ+1)) V)/(Subscript[ϵ, 0]^Δ))(∫_0^∞)dx((x^Δ)/((z^(-1))(e^x)-η)),

Mathematical Properties of (Subscript[ℐ, η]^s)(z)

Expression: (Subscript[ℐ, η]^s)(z) can be expressed in terms of the Euler Gamma function Γ(s) and the polylogarithm function (Subscript[Li, s])(z)

(Subscript[ℐ, η]^s)(z)=(1/η)Γ(s)(Subscript[Li, s])(η z),

where Γ(s) and (Subscript[Li, s])(z) are defined by

Γ(s)=(∫_0^∞)(x^(s-1))(e^(-x))dx,
(Subscript[Li, s])(z)=(∑_(k=1)^∞)((z^k)/(k^s)).

Verify Eq. () given the definition Eq. ().

Integrate[x^(s-1)/(Exp[x]/z-η),{x,0,∞},Assumptions->s>0]

(Subscript[ℐ, η]^s)(z)=(∫_0^∞)dx((x^(s-1))/((z^(-1))(e^x)-η))
=(∫_0^∞)dx (x^(s-1))((z (e^(-x)))/(1-η z (e^(-x))))
=(∫_0^∞)dx (x^(s-1))z (e^(-x))(∑_(k=0)^∞)(((η z (e^(-x))))^k)
=(1/(η ))(∫_0^∞)dx (x^(s-1))(∑_(k=1)^∞)(((η z (e^(-x))))^k)
=(1/(η ))(∑_(k=1)^∞)(((η z))^k)(∫_0^∞)dx (x^(s-1))(e^(-k x))

Γ(s)=(∫_0^∞)(x^(s-1))(e^(-x))dx,

(∫_0^∞)dx (x^(s-1))(e^(-k x))=(1/(k^s))(∫_0^∞)d(k x) (((k x))^(s-1))(e^(-k x))
(Overscript[=, k x→x])(1/(k^s))(∫_0^∞)dx (x^(s-1))(e^(-x))
=((Γ(s))/(k^s)).

(Subscript[ℐ, η]^s)(z)=(1/(η ))(∑_(k=1)^∞)(((η z))^k)((Γ(s))/(k^s))
=(1/(η ))Γ(s)(∑_(k=1)^∞)((((η z))^k)/(k^s)).

(Subscript[Li, s])(z)=(∑_(k=1)^∞)((z^k)/(k^s)),

(Subscript[ℐ, η]^s)(z)=(1/(η ))Γ(s)(Subscript[Li, s])(η z).

Γ(s)=(s-1)! can be expressed as a factorial. Therefore, it has an important property:

s Γ(s)=Γ(s+1).

Graph: Behavior of (Subscript[ℐ, η]^s)(z) for η=±1.

Block[{ℐ,plt},
ℐ[η_,s_,z_?NumericQ]:=Gamma[s]PolyLog[s,η z]/η;
Block[{r=3},
Plot[Evaluate[Tanh[ℐ[#,2,r ArcTanh[ζ]]/r]&/@{1,-1}],{ζ,0,1},Prolog->{{},{Dashed,AbsoluteThickness[1],Line[{#,#}&/@{0,1}]}},
Epilog->{Translate[{Blue,Text["Bosons",{0,0},{0,-0.6}],Text[Style["(\!\(η=+1\))",12],{0,0},{0,0.6}]},{0.18,0.75}],Translate[{Red,Text["Fermions",{0,0},{0,-0.6}],Text[Style["(\!\(η=-1\))",12],{0,0},{0,0.6}]},{0.75,0.2}]},
PlotRange->{{0,1},{0,1}},PlotRangePadding->0.05,PlotRangeClipping->False,GridLines->({#,#}&@Flatten@{{0},Tanh[Range[1,6]/r],{1}}),FrameLabel->{"\!\( z\)","\!\(ℐ\_η\%s(z)\)"},FrameTicks->({{#,Automatic},{#,Automatic}}&@{{0,0},{Tanh[1/r],1},{1,∞}}),ImagePadding->{{Automatic,10},{Automatic,10}},AspectRatio->1,ImageSize->160]]]

For bosons (η=+1): z∈(0,1) and the integral function (Subscript[ℐ, +]^s)(z) monotonically increases and stops (or diverges) at z=1.

For fermions (η=-1): z∈(0,+∞) and the integral function (Subscript[ℐ, -]^s)(z) monotonically increases towards +∞ as z→+∞.

When z->0, (Subscript[ℐ, η]^s)(z) for both bosons and fermions converges to the same linear behavior (Subscript[ℐ, η]^s)(z)∝z.

Derivative: (Subscript[ℐ, η]^s)(z) behaves like (((log z))^s)/s under differentiation

z((∂(Subscript[ℐ, η]^s)(z))/(∂z))≡((∂(Subscript[ℐ, η]^s)(z))/(∂ log z))=(s-1) (Subscript[ℐ, η]^(s-1))(z).

Verify Eq. ().

Block[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
FullSimplify[z D[ℐ[η,s,z],z]==(s-1) ℐ[η,s-1,z]]]

((∂(Subscript[ℐ, η]^s)(z))/(∂ log z))=z(∂_z) (Subscript[ℐ, η]^s)(z)
=z(∂_z)((1/η)Γ(s)(Subscript[Li, s])(η z))
=(1/η)Γ(s)z(∂_z)(Subscript[Li, s])(η z)
=(1/η)Γ(s)z(∂_z)(∑_(k=1)^∞)((((η z))^k)/(k^s))
=(1/η)Γ(s)η z(∑_(k=1)^∞)((k (((η z))^(k-1)))/(k^s))
=(1/η)Γ(s)(∑_(k=1)^∞)(( (((η z))^k))/(k^(s-1)))
=(1/η)Γ(s)(Subscript[Li, s-1])(η z)
=(1/η)(s-1) Γ(s-1)(Subscript[Li, s-1])(η z)
=(s-1) (Subscript[ℐ, η]^(s-1))(z).

Free Energy Derivatives

Module[{a=0.25},Graph[{"𝔉","V","T","z","μ"},{"V"\[DirectedEdge]"𝔉","T"\[DirectedEdge]"𝔉","z"\[DirectedEdge]"𝔉","T"\[DirectedEdge]"z","μ"\[DirectedEdge]"z"},VertexCoordinates->{{2,-1},{1,-3},{2,-3},{3,-2},{3,-3}},
VertexShapeFunction->Function[{p,v,b},Block[{hl=MemberQ[{},v]},{{EdgeForm[None],If[hl,Gray,White],Disk[p,0.25]},,{If[hl,White,Black],Text[Style["\!\( "<>v<>"\)","InlineFormula"],p,{0,0.2}]}}]],EdgeShapeFunction->Function[{ps,e},Block[{p0,p1,l,d},
p0=First@ps;p1=Last@ps;
d=Normalize[p1-p0];
{Black,Arrowheads[{0.08}],Arrow[{p0+a d,p1-a d}]}]],PlotRange->{{0.7,3.3},{-3.3,-0.7}},ImageSize->100]]

d𝔉=-PdV-SdT-Ndμ,

Pressure

Entropy

Derive Eq. ().

Module[{a=0.25},Graph[{"𝔉","V","T","log z","μ"},{"V"\[DirectedEdge]"𝔉","T"\[DirectedEdge]"𝔉","log z"\[DirectedEdge]"𝔉","T"\[DirectedEdge]"log z","μ"\[DirectedEdge]"log z"},VertexCoordinates->{{2,-1},{1,-3},{2,-3},{3,-2},{3,-3}},
VertexShapeFunction->Function[{p,v,b},Block[{hl=MemberQ[{},v]},{{EdgeForm[None],If[hl,Gray,White],Disk[p,0.25]},,{If[hl,White,Black],Text[Style["\!\( "<>v<>"\)","InlineFormula"],p,{0,0.2}]}}]],EdgeShapeFunction->Function[{ps,e},Block[{p0,p1,l,d},
p0=First@ps;p1=Last@ps;
d=Normalize[p1-p0];
{Black,Arrowheads[{0.08}],Arrow[{p0+a d,p1-a d}]}]],PlotRange->{{0.7,3.5},{-3.3,-0.7}},ImageSize->100]]

log z=(μ/T)=>(Subscript[(((∂log z)/(∂T))), μ])=-(μ/(T^2)).

S=-(Subscript[(((∂𝔉)/(∂T))), V,μ])
=-(Subscript[(((∂𝔉)/(∂T))), V,z])-(Subscript[(((∂𝔉)/(∂log z))), V,T])(Subscript[(((∂log z)/(∂T))), μ])
=-(Subscript[((∂/(∂T))(-T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z))), V,z])-(Subscript[((∂/(∂log z))(-T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z))), V,T])(-(μ/(T^2)))
=V(∂/(∂T))(T  ((((T/(Subscript[ϵ, 0]))))^Δ))(Subscript[ℐ, η]^(Δ+1))(z)-T V ((((T/(Subscript[ϵ, 0]))))^Δ)((μ/(T^2)))(((∂(Subscript[ℐ, η]^(Δ+1))(z))/(∂log z)))
=(Δ+1)V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z)-Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)((μ/T))(Subscript[ℐ, η]^Δ)(z)
=V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ+1)(Subscript[ℐ, η]^(Δ+1))(z)-Δ ((μ/T))(Subscript[ℐ, η]^Δ)(z))
=V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ+1)(Subscript[ℐ, η]^(Δ+1))(z)-Δ (log z) (Subscript[ℐ, η]^Δ)(z)).

Particle number

Derive Eq. ().

μ=T log z.

N=-(Subscript[(((∂𝔉)/(∂μ))), V,T])
=-(1/T)(Subscript[(((∂𝔉)/(∂log z))), V,T])
=-(1/T)(Subscript[((∂/(∂log z))(-T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z))), V,T])
=V ((((T/(Subscript[ϵ, 0]))))^Δ)((∂(Subscript[ℐ, η]^(Δ+1))(z))/(∂log z))
=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ) (Subscript[ℐ, η]^Δ)(z).

Module[{a=0.25},Graph[{"P","S","N","V","T","z","μ"},{"T"\[DirectedEdge]"P","z"\[DirectedEdge]"P","V"\[DirectedEdge]"S","T"\[DirectedEdge]"S","z"\[DirectedEdge]"S","V"\[DirectedEdge]"N","T"\[DirectedEdge]"N","z"\[DirectedEdge]"N","T"\[DirectedEdge]"z","μ"\[DirectedEdge]"z"},VertexCoordinates->{{1,-1},{2,-1},{3,-1},{1,-3},{2,-3},{3,-2},{3,-3}},VertexShapeFunction->Function[{p,v,b},Block[{hl=MemberQ[{},v]},{{EdgeForm[None],If[hl,Gray,White],Disk[p,0.25]},,{If[hl,White,Black],Text[Style["\!\( "<>v<>"\)","InlineFormula"],p,{0,0.2}]}}]],EdgeShapeFunction->Function[{ps,e},Block[{p0,p1,l,d},
p0=First@ps;p1=Last@ps;
d=Normalize[p1-p0];
{Black,Arrowheads[{0.08}],Arrow[{p0+a d,p1-a d}]}]],PlotRange->{{0.7,3.3},{-3.3,-0.7}},ImageSize->100]]

Derive Eq. ().

E=𝔉+T S+μ N
=-T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z)+T V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ+1)(Subscript[ℐ, η]^(Δ+1))(z)-Δ ((μ/T))(Subscript[ℐ, η]^Δ)(z))+μ Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^Δ)(z)
=Δ T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z).

For gas of non-relativistic particles, E=(3/2)P V.

For gas of relativistic particles, E=3 P V.

Equation of State

A well-known example is the (classical) ideal gas law:

P V=N T,

where the Boltzmann constant has been set to (Subscript[k, B])=1. However, it only describe the ideal gas in the classical limit (low density and high temperature).

Comparing Eq. () and Eq. (), we found the equation of state for ideal quantum gas (free bosons and fermions)

P V=(Subscript[k, eff])N T,

where the P V/N T ratio is a function of fugacity z=(e^(μ/T)) only, denoted as

(Subscript[k, eff]):=((P V)/(N T))=(((Subscript[ℐ, η]^(Δ+1))(z))/(Δ(Subscript[ℐ, η]^Δ)(z))).

It can be viewed as an effective Boltzmann constant (we will see that (Subscript[k, eff])=1 for ideal classical gas as z→0).

Dilute (Classical) Limit

Block[{ℐ,plt},
ℐ[η_,s_,z_?NumericQ]:=Gamma[s]PolyLog[s,η z]/η;
Block[{r=3},
Plot[Evaluate[Tanh[ℐ[#,2,r ArcTanh[ζ]]/r]&/@{1,-1}],{ζ,0,1},Prolog->{{{Opacity[0.2],Black,Disk[{0,0},0.1]},{Black,Text["Dilute limit",{0.1,0},{-1.1,0}]}},{Dashed,AbsoluteThickness[1],Line[{#,#}&/@{0,1}]}},
Epilog->{Translate[{Blue,Text["Bosons",{0,0},{0,-0.6}],Text[Style["(\!\(η=+1\))",12],{0,0},{0,0.6}]},{0.18,0.75}],Translate[{Red,Text["Fermions",{0,0},{0,-0.6}],Text[Style["(\!\(η=-1\))",12],{0,0},{0,0.6}]},{0.75,0.2}]},
PlotRange->{{0,1},{0,1}},PlotRangePadding->0.05,PlotRangeClipping->False,GridLines->({#,#}&@Flatten@{{0},Tanh[Range[1,6]/r],{1}}),FrameLabel->{"\!\( z\)","\!\(ℐ\_η\%s(z)\)"},FrameTicks->({{#,Automatic},{#,Automatic}}&@{{0,0},{Tanh[1/r],1},{1,∞}}),ImagePadding->{{Automatic,10},{Automatic,10}},AspectRatio->1,ImageSize->160]]]

Verify Eq. ().

Block[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
Series[ℐ[η,s,z],{z,0,2}]]

(Subscript[ℐ, η]^s)(z)=(1/η)Γ(s)(Subscript[Li, s])(η z),

(Subscript[Li, s])(z)=(∑_(k=1)^∞)((z^k)/(k^s)).

(Subscript[ℐ, η]^s)(z)=(1/η)Γ(s)(∑_(k=1)^∞)((((η z))^k)/(k^s))
=Γ(s)(∑_(k=1)^∞)(((η^(k-1)) (z^k))/(k^s)).

(Subscript[ℐ, η]^s)(z)=Γ(s)((((η^(1-1)) (z^1))/(1^s))+(((η^(2-1)) (z^2))/(2^s))+…)
=Γ(s)(z+((η (z^2))/(2^s))+…).

(Subscript[k, eff])=1-((η z)/(2^(Δ+1)))+...,
N=V ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ+1)z+….

Verify the expansions in Eq. ().

Block[{ℐ,k},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
k[η_,s_,z_]:=ℐ[η,Δ+1,z]/(Δ ℐ[η,Δ,z]);
FullSimplify/@{Series[k[η,Δ,z],{z,0,1}],Series[Δ V (T/ϵ0)^Δ ℐ[η,Δ,z],{z,0,1}]}]

For (Subscript[k, eff]):

(Subscript[k, eff])=(((Subscript[ℐ, η]^(Δ+1))(z))/(Δ(Subscript[ℐ, η]^Δ)(z)))
=((Γ(Δ+1) (z+(η/(2^(Δ+1)))(z^2)+…))/(Δ Γ(Δ) (z+(η/(2^Δ))(z^2)+…)))
=((Γ(Δ+1) z (1+(η/(2^(Δ+1)))z+…))/(Γ(Δ+1)z  (1+(η/(2^Δ))z+…)))
=(1+(η/(2^(Δ+1)))z+…)(((1+(η/(2^Δ))z+…))^(-1))
=1+(η/(2^(Δ+1)))z-(η/(2^Δ))z+…
=1+((1/2)-1)(η/(2^Δ))z+…
=1-(η/(2^(Δ+1)))z+…

For N:

N=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^Δ)(z)
=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ)(z+(η/(2^Δ))(z^2)+…).

The (z^2) term can be ignored, as the expansion to the order of z is already sufficient for obtaining the virial expansion to the leading order.

The deviation from the ideal classical gas behavior (P V=N T) will be significant when (N/V)((((Subscript[ϵ, 0])/T))^Δ)≫1, i.e. when

the temperature is low (T→0),

or the density is high (N/V→∞).

Bose statistics (η=+1) causes an decrease in pressure compare to the classical limit, as if there is an effective attraction between bosons.

DynamicModule[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
DynamicModule[{η=1,Δ=3/2,PTz,sTz,nTz,znT,snT,PnT},
PTz[T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
sTz[T_,z_]:=T^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]);
nTz[T_,z_]:=Δ T^Δ ℐ[η,Δ,z];
znT[n_,T_]:=Block[{z},
Re@Switch[η,1,If[n>=nTz[T,1],1,z/.FindRoot[nTz[T,z]==n,{z,0,1}]],-1,z/.FindRoot[nTz[T,z]==n,{z,n/(T^Δ Gamma[Δ+1]),Exp[n^(1/Δ)/T]}]]];
snT[n_,T_]:=sTz[T,znT[n,T]];
PnT[n_,T_]:=PTz[T,znT[n,T]];
Manipulate[Plot[{PnT[n,T],n T},{T,0.01,2},PlotRange->{0,2},PlotRangePadding->0.1,PlotLabel->StringTemplate["\!\(N/V=`1`\)"][n],
PlotStyle->{Blue,Directive[Black,AbsoluteThickness[1],Dashed]},FrameLabel->{"\!\(T/ϵ\_0\)","\!\(P/ϵ\_0\)"},ImageSize->180],{{n,1,Style["\!\(N/V\)","InlineFormula"]},0.2,2,ImageSize->Small},Paneled->False,SaveDefinitions->True]]
]

Fermi statistics (η=-1) causes an increase in pressure compare to the classical limit, as if there is an effective repulsion between fermions.

DynamicModule[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
DynamicModule[{η=-1,Δ=3/2,PTz,sTz,nTz,znT,snT,PnT},
PTz[T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
sTz[T_,z_]:=T^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]);
nTz[T_,z_]:=Δ T^Δ ℐ[η,Δ,z];
znT[n_,T_]:=Block[{z},
Re@Switch[η,1,If[n>=nTz[T,1],1,z/.FindRoot[nTz[T,z]==n,{z,0,1}]],-1,z/.FindRoot[nTz[T,z]==n,{z,n/(T^Δ Gamma[Δ+1]),Exp[n^(1/Δ)/T]}]]];
snT[n_,T_]:=sTz[T,znT[n,T]];
PnT[n_,T_]:=PTz[T,znT[n,T]];
Manipulate[Plot[{PnT[n,T],n T},{T,0.01,2},PlotRange->{0,2},PlotRangePadding->0.1,PlotLabel->StringTemplate["\!\(N/V=`1`\)"][n],
PlotStyle->{Red,Directive[Black,AbsoluteThickness[1],Dashed]},FrameLabel->{"\!\(T/ϵ\_0\)","\!\(P/ϵ\_0\)"},ImageSize->180],{{n,1,Style["\!\(N/V\)","InlineFormula"]},0.2,2,ImageSize->Small},Paneled->False,SaveDefinitions->True]]
]

Isothermal Process

P(T,z)=T  ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z),
V(T,z)=N (((Δ  ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^Δ)(z)))^(-1)).

For bosons (η=+1), the isothermal curves are like

Block[{ℐ,PTz,VTz},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
PTz[η_,Δ_,T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
VTz[η_,Δ_,T_,z_]:=1/(Δ T^Δ ℐ[η,Δ,z]);
Block[{η=1,Δ=3/2,Ts,cols,legend,isotherms,degbdy},
Ts={0.2,0.4,0.6,0.8,1,1.2,1.4,2,4,6};
cols=Blend[{Green,Yellow,Orange},#]&/@Sqrt@Rescale[Ts,{Min@Ts,Max@Ts}];
legend=Column[{Style["\!\(T/ϵ\_0\)","InlineFormula"],Grid[{Graphics[{#1,Line[{#,0}&/@{0,1}]},AspectRatio->0.1,ImageSize->18],Style[ToString@NumberForm[N@#2,{2,1}],12,"InlineFormula"]}&@@@Thread@{cols,Ts},Spacings->{0.2,0}]}];
isotherms=First@ParametricPlot[Evaluate@Table[{VTz[η,Δ,T,Exp[μ/T]],PTz[η,Δ,T,Exp[μ/T]]},{T,Ts}],{μ,-30,0},PlotStyle->cols]/.{Line[pts_]:>Line[Append[pts,{0,Last@Last@pts}]]};
degbdy=Line[Table[{V,(1/V)^(1+1/Δ)Gamma[Δ+1]Zeta[Δ+1]/(Gamma[Δ+1]Zeta[Δ])^(1+1/Δ)},{V,3,0.1,-0.1}]];
Row[{Graphics[{{Lighter[Blue,0.8],FilledCurve[{degbdy,Line[{{0,3},{0,0},{3,0}}]}]},{AbsoluteThickness[1],Blue,Dashed,degbdy},isotherms},PlotRange->{{0,3},{0,3}},PlotRangePadding->Scaled[0.03],PlotRangeClipping->True,Frame->True,FrameLabel->{"\!\( V/N\)","\!\(P/ϵ\_0\)"},PlotLabel->Style[StringTemplate["\!\(Δ=`1`\)"][Δ],14,Black],BaselinePosition->Scaled[0.55],ImageSize->180],legend}]]]

For fermions (η=-1), the isothermal curves are like

Block[{ℐ,PTz,VTz},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
PTz[η_,Δ_,T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
VTz[η_,Δ_,T_,z_]:=1/(Δ T^Δ ℐ[η,Δ,z]);
Block[{η=-1,Δ=3/2,Ts,cols,legend,isotherms,degbdy},
Ts={0.2,0.4,0.6,0.8,1,1.2,1.4,2,4,6};
cols=Blend[{Green,Yellow,Orange},#]&/@Sqrt@Rescale[Ts,{Min@Ts,Max@Ts}];
legend=Column[{Style["\!\(T/ϵ\_0\)","InlineFormula"],Grid[{Graphics[{#1,Line[{#,0}&/@{0,1}]},AspectRatio->0.1,ImageSize->18],Style[ToString@NumberForm[N@#2,{2,1}],12,"InlineFormula"]}&@@@Thread@{cols,Ts},Spacings->{0.2,0}]}];
isotherms=First@ParametricPlot[Evaluate@Table[{VTz[η,Δ,T,Exp[μ/T]],PTz[η,Δ,T,Exp[μ/T]]},{T,Ts}],{μ,-30,30},PlotStyle->cols]/.{Line[pts_]:>Line[Append[pts,{0,Last@Last@pts}]]};
degbdy=Line[Table[{V,1/(Δ+1)/V^(1+1/Δ)},{V,3,0.1,-0.1}]];
Row[{Graphics[{{Lighter[Red,0.8],FilledCurve[{degbdy,Line[{{0,3},{0,0},{3,0}}]}]},{AbsoluteThickness[1],Red,Dashed,degbdy},isotherms},PlotRange->{{0,3},{0,3}},PlotRangePadding->Scaled[0.03],PlotRangeClipping->True,Frame->True,FrameLabel->{"\!\( V/N\)","\!\(P/ϵ\_0\)"},PlotLabel->Style[StringTemplate["\!\(Δ=`1`\)"][Δ],14,Black],BaselinePosition->Scaled[0.55],ImageSize->180],legend}]]]

Degenerated Bose gas: the curves flatten out at low volumes, suggesting transitions to Bose-Einstein condensates (BEC), where pressure cease to increase as more and more bosons condense to a zero-momentum state that does not contribute to pressure.

Degenerated Fermi gas: even at zero temperature, fermions can sustain a finite pressure, called the degeneracy pressure, due to the Pauli exclusion principle preventing fermions from occupying the same quantum state.

Specific Heat

Isochoric (constant-volume) specific heat

It can be computed by evaluating S as a function of T given N,V fixed.

Module[{a=0.25},Graph[{"P","S","N","V","T","z","μ"},{"T"\[DirectedEdge]"P","z"\[DirectedEdge]"P","V"\[DirectedEdge]"S","T"\[DirectedEdge]"S","z"\[DirectedEdge]"S","V"\[DirectedEdge]"N","T"\[DirectedEdge]"N","z"\[DirectedEdge]"N","T"\[DirectedEdge]"z","μ"\[DirectedEdge]"z"},VertexCoordinates->{{1,-1},{2,-1},{3,-1},{1,-3},{2,-3},{3,-2},{3,-3}},VertexShapeFunction->Function[{p,v,b},Block[{hl=MemberQ[{"N","V"},v]},{{EdgeForm[None],If[hl,Gray,White],Disk[p,0.25]},,{If[hl,White,Black],Text[Style["\!\( "<>v<>"\)","InlineFormula"],p,{0,0.2}]}}]],EdgeShapeFunction->Function[{ps,e},Block[{p0,p1,l,d},
p0=First@ps;p1=Last@ps;
d=Normalize[p1-p0];
{Black,Arrowheads[{0.08}],Arrow[{p0+a d,p1-a d}]}]],PlotRange->{{0.7,3.3},{-3.3,-0.7}},ImageSize->100]]

Given N,V, input T.

Solve for z from N(V,T,z)=N.

Output S=S(V,T,z).

Isobaric (constant-pressure) specific heat

It can be computed by evaluating S as a function of T given N,P fixed.

Module[{a=0.25},Graph[{"P","S","N","V","T","z","μ"},{"T"\[DirectedEdge]"P","z"\[DirectedEdge]"P","V"\[DirectedEdge]"S","T"\[DirectedEdge]"S","z"\[DirectedEdge]"S","V"\[DirectedEdge]"N","T"\[DirectedEdge]"N","z"\[DirectedEdge]"N","T"\[DirectedEdge]"z","μ"\[DirectedEdge]"z"},VertexCoordinates->{{1,-1},{2,-1},{3,-1},{1,-3},{2,-3},{3,-2},{3,-3}},VertexShapeFunction->Function[{p,v,b},Block[{hl=MemberQ[{"N","P"},v]},{{EdgeForm[None],If[hl,Gray,White],Disk[p,0.25]},,{If[hl,White,Black],Text[Style["\!\( "<>v<>"\)","InlineFormula"],p,{0,0.2}]}}]],EdgeShapeFunction->Function[{ps,e},Block[{p0,p1,l,d},
p0=First@ps;p1=Last@ps;
d=Normalize[p1-p0];
{Black,Arrowheads[{0.08}],Arrow[{p0+a d,p1-a d}]}]],PlotRange->{{0.7,3.3},{-3.3,-0.7}},ImageSize->100]]

Given N,P, input T.

Solve for z from P(T,z)=P.

Solve for V from N(V,T,z)=N.

Output S=S(V,T,z).

Manipulate[Block[{η=1,ℐ,PTz,sTz,nTz,znT},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
PTz[T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
sTz[T_,z_]:=T^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]);
nTz[T_,z_]:=Δ T^Δ ℐ[η,Δ,z];
znT[n_,T_,z0_:1]:=Block[{z},
If[η==1&&Δ>1&&n>=nTz[T,1],1,z/.FindRoot[nTz[T,z]==n,{z,z0}]]];
Block[{Tmax=3,Ts,zs,ss,sT,data},
Ts=Tmax E Exp[-1/Range[0.16,1,0.005]];
Off[General::munfl,FindRoot::lstol];
zs=Re@Reverse@Rest@FoldList[znT[n,#2,#1]&,0,Reverse@Ts];
ss=Re[sTz[#1,#2]/n]&@@@Thread@{Ts,zs};
sT=Interpolation[Prepend[Thread@{Ts,ss},{0,0}],InterpolationOrder->1];
data=Prepend[Table[{T,Chop[T sT'[T]]},{T,Mean/@Thread@{Most@Ts,Rest@Ts}}],{0,0}];
ListLinePlot[data,PlotStyle->Blue,PlotRange->{0,Max[Max[data[[All,2]],Δ]]},PlotRangePadding->{Scaled[0.02],Scaled[0.05]},FrameLabel->{"\!\(T/ϵ\_0\)","\!\(c\_V\)"},Epilog->{Black,AbsoluteThickness[1],Dashed,Line[{#,Δ}&/@{0,Tmax}]},PlotLabel->StringTemplate["\!\(N/V=`1`\)"][n],ImageSize->{Automatic,170}]
]],
{{Δ,3/2,Style["\!\(Δ\)","InlineFormula"]},{1/2,1,3/2,2,3}},{{n,1,Style["\!\(N/V\)","InlineFormula"]},0.2,2,ImageSize->Small},Paneled->False,SaveDefinitions->True]

Manipulate[Block[{η=1,ℐ,PTz,sTz,nTz,zPT,znT},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
PTz[T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
sTz[T_,z_]:=T^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]);
nTz[T_,z_]:=Δ T^Δ ℐ[η,Δ,z];
zPT[P_,T_,z0_:1]:=Block[{z},
If[η==1&&Δ>0&&P>=PTz[T,1],0.9999,z/.FindRoot[PTz[T,z]==P,{z,z0}]]];
Block[{Tmax=3,Ts,zs,ss,sT,data},
Ts=Tmax E Exp[-1/Range[0.16,1,0.005]];
Off[General::munfl,FindRoot::lstol,FindRoot::nlnum];
zs=Re@Reverse@Rest@FoldList[zPT[P,#2,#1]&,0,Reverse@Ts];
ss=Re[sTz[#1,#2]/nTz[#1,#2]]&@@@Thread@{Ts,zs};
sT=Interpolation[Prepend[Thread@{Ts,ss},{0,0}],InterpolationOrder->1];
data=Prepend[Table[{T,Chop[T sT'[T]]},{T,Mean/@Thread@{Most@Ts,Rest@Ts}}],{0,0}];
ListLinePlot[data,
PlotStyle->Blue,PlotRange->{0,Max[Min[Max[data[[All,2]]],2(Δ+1)],Δ+1]},PlotRangePadding->{Scaled[0.02],Scaled[0.05]},FrameLabel->{"\!\(T/ϵ\_0\)","\!\(c\_P\)"},Epilog->{Black,AbsoluteThickness[1],Dashed,Line[{#,Δ+1}&/@{0,Tmax}]},PlotLabel->StringTemplate["\!\(P=`1`\)"][P],ImageSize->{Automatic,170}]
]],
{{Δ,3/2,Style["\!\(Δ\)","InlineFormula"]},{1/2,1,3/2,2,3}},{{P,1,Style["\!\( P\)","InlineFormula"]},0.2,2,ImageSize->Small},Paneled->False,SaveDefinitions->True]

Manipulate[Block[{η=-1,ℐ,PTz,sTz,nTz,znT},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
PTz[T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
sTz[T_,z_]:=T^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]);
nTz[T_,z_]:=Δ T^Δ ℐ[η,Δ,z];
znT[n_,T_,z0_:1]:=Block[{z},
If[η==1&&n>=nTz[T,1],1,z/.FindRoot[nTz[T,z]==n,{z,z0}]]];
Block[{Tmax=3,Ts,zs,ss,sT,data},
Ts=Tmax E Exp[-1/Range[0.16,1,0.02]];
Off[General::munfl,FindRoot::lstol];
zs=Re@Reverse@Rest@FoldList[znT[n,#2,#1]&,0,Reverse@Ts];
ss=Re[sTz[#1,#2]/n]&@@@Thread@{Ts,zs};
sT=Interpolation[Prepend[Thread@{Ts,ss},{0,0}],InterpolationOrder->1];
data=Prepend[Table[{T,Chop[T sT'[T]]},{T,Mean/@Thread@{Most@Ts,Rest@Ts}}],{0,0}];
ListLinePlot[data,PlotStyle->Red,PlotRange->{0,Max[Max[data[[All,2]]],Δ]},PlotRangePadding->{Scaled[0.02],Scaled[0.05]},FrameLabel->{"\!\(T/ϵ\_0\)","\!\(c\_V\)"},Epilog->{Black,AbsoluteThickness[1],Dashed,Line[{#,Δ}&/@{0,Tmax}]},PlotLabel->StringTemplate["\!\(N/V=`1`\)"][n],ImageSize->{Automatic,170}]
]],
{{Δ,3/2,Style["\!\(Δ\)","InlineFormula"]},{1/2,1,3/2,2,3}},{{n,1,Style["\!\(N/V\)","InlineFormula"]},0.2,2,ImageSize->Small},Paneled->False,SaveDefinitions->True]

Manipulate[Block[{η=-1,ℐ,PTz,sTz,nTz,zPT},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
PTz[T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
sTz[T_,z_]:=T^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]);
nTz[T_,z_]:=Δ T^Δ ℐ[η,Δ,z];
zPT[P_,T_,z0_:1]:=Block[{z},
If[η==1&&P>=PTz[T,1],1,z/.FindRoot[PTz[T,z]==P,{z,z0}]]];
Block[{Tmax=3,Ts,zs,ss,sT,data},
Ts=Tmax E Exp[-1/Range[0.16,1,0.02]];
Off[General::munfl,FindRoot::lstol,FindRoot::nlnum];
zs=Re@Reverse@Rest@FoldList[zPT[P,#2,#1]&,0,Reverse@Ts];
ss=Re[sTz[#1,#2]/nTz[#1,#2]]&@@@Thread@{Ts,zs};
sT=Interpolation[Prepend[Thread@{Ts,ss},{0,0}],InterpolationOrder->1];
data=Prepend[Table[{T,Chop[T sT'[T]]},{T,Mean/@Thread@{Most@Ts,Rest@Ts}}],{0,0}];
ListLinePlot[data,
PlotStyle->Red,PlotRange->{0,Max[Max[data[[All,2]]],Δ+1]},PlotRangePadding->{Scaled[0.02],Scaled[0.05]},FrameLabel->{"\!\(T/ϵ\_0\)","\!\(c\_P\)"},Epilog->{Black,AbsoluteThickness[1],Dashed,Line[{#,Δ+1}&/@{0,Tmax}]},PlotLabel->StringTemplate["\!\(P=`1`\)"][P],ImageSize->{Automatic,170}]
]],
{{Δ,3/2,Style["\!\(Δ\)","InlineFormula"]},{1/2,1,3/2,2,3}},{{P,1,Style["\!\( P\)","InlineFormula"]},0.2,2,ImageSize->Small},Paneled->False,SaveDefinitions->True]

Dilute (Classical) Limit

Block[{ℐ,plt},
ℐ[η_,s_,z_?NumericQ]:=Gamma[s]PolyLog[s,η z]/η;
Block[{r=3},
Plot[Evaluate[Tanh[ℐ[#,2,r ArcTanh[ζ]]/r]&/@{1,-1}],{ζ,0,1},Prolog->{{{Opacity[0.2],Black,Disk[{0,0},0.1]},{Black,Text["Dilute limit",{0.1,0},{-1.1,0}]}},{Dashed,AbsoluteThickness[1],Line[{#,#}&/@{0,1}]}},
Epilog->{Translate[{Blue,Text["Bosons",{0,0},{0,-0.6}],Text[Style["(\!\(η=+1\))",12],{0,0},{0,0.6}]},{0.18,0.75}],Translate[{Red,Text["Fermions",{0,0},{0,-0.6}],Text[Style["(\!\(η=-1\))",12],{0,0},{0,0.6}]},{0.75,0.2}]},
PlotRange->{{0,1},{0,1}},PlotRangePadding->0.05,PlotRangeClipping->False,GridLines->({#,#}&@Flatten@{{0},Tanh[Range[1,6]/r],{1}}),FrameLabel->{"\!\( z\)","\!\(ℐ\_η\%s(z)\)"},FrameTicks->({{#,Automatic},{#,Automatic}}&@{{0,0},{Tanh[1/r],1},{1,∞}}),ImagePadding->{{Automatic,10},{Automatic,10}},AspectRatio->1,ImageSize->160]]]

P=T ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ+1)z+…,

S=V ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ+1)(1+Δ-log z)z+…,

N=V ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ+1)z+….

Verify the expansions in Eq. (, , ).

Block[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
FullSimplify[Series[#,{z,0,1}]]&/@{T(T/ϵ0)^Δ ℐ[η,Δ+1,z],V(T/ϵ0)^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]),Δ V (T/ϵ0)^Δ ℐ[η,Δ,z]}]

To compute the isochoric (constant volume) specific heat, N and V are held fixed. First solve z from Eq. (),

z=(1/(Γ(Δ+1)))(N/V)(((((Subscript[ϵ, 0])/T)))^Δ),

then substitute Eq. () to Eq. (), to the leading order,

S=N  (Δ log T+log V-log N)+const.

Derive Eq. ().

S=V ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ+1)(1+Δ-log ((1/(Γ(Δ+1)))(N/V)(((((Subscript[ϵ, 0])/T)))^Δ)))(1/(Γ(Δ+1)))(N/V)(((((Subscript[ϵ, 0])/T)))^Δ)
=N (1+Δ-log ((1/(Γ(Δ+1)))(N/V)(((((Subscript[ϵ, 0])/T)))^Δ)))
=N (1+Δ-log N-Δ log (Subscript[ϵ, 0])+log Γ(Δ+1)+log V+Δ log T)
=N(Δ log T+log V-log N)+const.

therefore, (Subscript[c, V])=(Subscript[(∂S/∂log T), N,V])/N is

To compute the isobaric (constant pressure) specific heat, N and P must be fixed. First solve z from Eq. (),

z=(1/(Γ(Δ+1)))(P/T)(((((Subscript[ϵ, 0])/T)))^Δ),

then substitute Eq. () to Eq. (), to the leading order,

S=N  ((Δ+1) log T-log P)+const.

Derive Eq. ().

S=V ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ+1)(1+Δ-log ((1/(Γ(Δ+1)))(P/T)(((((Subscript[ϵ, 0])/T)))^Δ)))(1/(Γ(Δ+1)))(P/T)(((((Subscript[ϵ, 0])/T)))^Δ)
=((P V)/T) (1+Δ-log ((1/(Γ(Δ+1)))(P/T)(((((Subscript[ϵ, 0])/T)))^Δ)))
=((P V)/T) (1+Δ-log P-Δ log (Subscript[ϵ, 0])+log Γ(Δ+1)+(Δ+1)log T)
=N ((Δ+1)log T-log P)+const.

therefore, (Subscript[c, P])=(Subscript[(∂S/∂log T), N,P])/N is

(Subscript[c, V])=Δ+(η/((2^(Δ+1))Γ(Δ-1)))(N/V)(((((Subscript[ϵ, 0])/T)))^Δ)+…,
(Subscript[c, P])=Δ+1+((η (Δ+1))/((2^(Δ+1))Γ(Δ)))(P/T)(((((Subscript[ϵ, 0])/T)))^Δ)+….

Derive Eq. ().

Block[{ℐ,Pz,Sz,Νz,zΝ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
{Pz,Sz,Νz}=FullSimplify[Series[#,{z,0,2}]]&/@{T(T/ϵ0)^Δ ℐ[η,Δ+1,z],V(T/ϵ0)^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]),Δ V (T/ϵ0)^Δ ℐ[η,Δ,z]};
zΝ=InverseSeries[Νz,Ν];
FullSimplify[(T/Ν)D[Sz/.{z->zΝ},T]]]

Block[{ℐ,Pz,Sz,Νz,zP},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
{Pz,Sz,Νz}=FullSimplify[Series[#,{z,0,2}]]&/@{T(T/ϵ0)^Δ ℐ[η,Δ+1,z],V(T/ϵ0)^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]),Δ V (T/ϵ0)^Δ ℐ[η,Δ,z]};
zP=InverseSeries[Pz,P];
FullSimplify[T D[Sz/Νz/.{z->zP},T]]]

Adiabatic Process

Prove the statement Eq. (, ).

P=T  ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z),
S=V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ+1)(Subscript[ℐ, η]^(Δ+1))(z)-((Δ μ)/T)(Subscript[ℐ, η]^Δ)(z)),
N=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^Δ)(z),

(S/N)=(1+(1/Δ))(((Subscript[ℐ, η]^(Δ+1))(z))/((Subscript[ℐ, η]^Δ)(z)))- log z,

P (V^(1+1/Δ))=(T  ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z))((((N/(Δ ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^Δ)(z)))))^((Δ+1)/Δ))
=(Subscript[ℐ, η]^(Δ+1))(z)((((N/(Δ (Subscript[ℐ, η]^Δ)(z)))))^((Δ+1)/Δ))(Subscript[ϵ, 0])=const.,

Summary

Grid[{{"", "","Bose gas","Fermi gas","Classical limit"},
{"Energy","\!\(E=Δ P V\)","✓","✓","✓"},{"Eq. of state","\!\(P V = k\_eff N T\)","\!\(k\_eff<1\)","\!\(k\_eff>1\)","\!\(k\_eff=1\)"},
{"Specific heat","\!\(c\_V\)","\!\((∼T\^Δ)\_\(T→0\)\)","\!\((∼ T)\_\(T→0\)\)","\!\(Δ\)"},
{SpanFromAbove,"\!\(c\_P\)","\!\(0\_\(T→0\)\)","\!\((∼ T)\_\(T→0\)\)","\!\(Δ+1\)"},
{"Isochoric","\!\(V=const.\)","✓","✓","✓"},
{"Isobaric","\!\(P=const.\)","✓","✓","✓"},{"Isothermal","\!\(P V=const.\)",Style["×",Red],Style["×",Red],"✓"},
{"Adiabatic","\!\(P V\^γ=const.\)\n(\!\(γ=1+1/Δ\))","✓","✓","✓"}},Dividers->{None,{2->True,4->True,6->True}},ItemStyle->{Automatic,{"Text",{"InlineFormula"}}}]

In reality, gases are consist of quantum particles (either bosons or fermions) at the microscopic level. => The classical gas is not real.

However, in the dilute limit, both Bose and Fermi gases converge to a universal behavior described by the classical gas. Given z∝(N/V)((((Subscript[ϵ, 0])/T))^Δ) [see Eq. ()], the classical limit (z→0) is characterized by

(N/V)(((((Subscript[ϵ, 0])/T)))^Δ)≪1.

Degenerated Quantum Gases

Degenerated Bose Gas

Unsustainable Thermal Bosons

(N/V)=(∫_0^∞)dϵ ((g(ϵ))/((e^( β (ϵ-μ)))-1)).

g(ϵ)=Δ (ϵ^(Δ-1))/(Subscript[ϵ, 0]^Δ) is the (per volume) density of states.

As the temperature T=1/β lowers, the chemical potential  μ (or the fugacity z=(e^(μ/T))) will adjust to maintain the number density N/V (say N/V=1).

DynamicModule[{ℐ,nTz,legend},
ℐ[Δ_,z_]:=Gamma[Δ+1]PolyLog[Δ+1,z];
nTz[Δ_,T_,z_]:=Re[Δ T^Δ ℐ[Δ-1,z]];
legend=Grid[{Graphics[{#1,Line[{#,0}&/@{-1,0}]},AspectRatio->0.1,ImageSize->18],"\!\("<>#2<>"\)"}&@@@{{Green,"g(ϵ)"},{Blue,"n\_B(ϵ)"}},Spacings->{0.2,0}];
Manipulate[Block[{n=1,z0=0.9999,z,g},
z=Which[Δ==1&&n/T>5,1-SetPrecision[Exp[-n/T],20],
Δ>1&&n>=nTz[Δ,T,1],z0,True,z/.FindRoot[nTz[Δ,T,z]==n,{z,z0}]];
g[ϵ_]:=Δ Abs[ϵ]^(Δ-1)UnitStep[ϵ];
Row[{Plot[{g[ϵ],1/(Exp[ϵ/T]/z-1),g[ϵ]/(Exp[ϵ/T]/z-1)},{ϵ,-1,5},Epilog->{Inset[legend,Scaled[{0.95,0.95}],{1,1}]},PlotStyle->{Green,Blue,Orange},Filling->{3->Axis},Exclusions->None,PlotRange->{{-0.4,3},{0,2}},FrameLabel->{"\!\(ϵ/ϵ\_0\)","Integrand"},PlotLabel->StringTemplate["\!\(T/ϵ\_0=`1`\)"][T],ImageSize->200],
Graphics[Block[{w=0.4},{Translate[{EdgeForm[Red],Lighter[Red,0.5],Rectangle[{-w,0},{w,z}]},{1,0}],Translate[{{EdgeForm[Cyan],Lighter[Cyan,0.5],Rectangle[{-w,#},{w,n}]},{EdgeForm[Orange],Lighter[Orange,0.5],Rectangle[{-w,0},{w,#}]}}&@nTz[Δ,T,z],{2,0}]}](*{{z},{#,1-#}&@nTz[Δ,T,z]}*),PlotRange->{{0.5,2.5},{0,1}},PlotRangePadding->{Scaled[.02],Scaled[.02]},Frame->True,FrameTicks->{{Automatic,Automatic},{{{1,"\!\( z\)"},{2,"\!\(N/V\)"}},None}},AspectRatio->2,ImageSize->{Automatic,120}]}]
],
{{Δ,3/2,Style["\!\( Δ\)","InlineFormula"]},{1/2,1,3/2,2,3}},{{T,1,Style["\!\( T\)","InlineFormula"]},0.01,1.5},
Paneled->False,SaveDefinitions->True]]

It is possible that the integral can not provide the required N/V, even if the chemical potential is tuned to  μ=0. To see this, we investigate the divergent behavior of the integrand near ϵ->0,

(N/V)∝(∫_0^∞)dϵ ((ϵ^(Δ-1))/((e^( ϵ/T))-1))(Overscript[~, ϵ→0]) T(∫_0^(...))((dϵ)/(ϵ^(2-Δ))).

If Δ>1, the integral will converge (even if  μ=0), which can not provide the required N/V as T lowers to a certain point.

If Δ≤1, the integral can diverge (if  μ=0), which enables the adjustment of  μ to reach the required N/V no matter how low T is.

The integral Eq. () only accounts for the thermal bosons (bosons occupying the excited states).

Due to the vanishing DOS as ϵ->0, the bosons occupying the lowest-energy ground state (ϵ=0) was not taken into account.

(Subscript[N, th]): the thermal boson number, defined by the integral Eq. (), whose result is given by Eq. (),

(Subscript[N, th])=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^Δ)(z).

(Subscript[N, 0]): the condensed boson number (in the ground state), given by the Bose-Einstein distribution (Subscript[n, B])(ϵ)=((((z^(-1))(e^(β ϵ))-1))^(-1)) at ϵ=0:

(Subscript[N, 0])=(1/((z^(-1))-1)),

meaning that the fugacity will be set by z=(Subscript[N, 0])/((Subscript[N, 0])+1).

Bose-Einstein Condensation

By saying that the ground state occupation number (Subscript[N, 0]) is macroscopic, we mean

(lim_(V->∞))((Subscript[N, 0])/V)≠0.

This is unusual, as for other excited states, the occupation number does not scale with the system size V (i.e. ⟨(Subscript[n, α])⟩/V→0 as V→∞).

The fraction of thermal v.s. condensed (ground-state) bosons:

DynamicModule[{ℐ,nTz,legend},
ℐ[Δ_,z_]:=Gamma[Δ+1]PolyLog[Δ+1,z];
nTz[Δ_,T_,z_]:=Re[Δ T^Δ ℐ[Δ-1,z]];
legend=Grid[{Graphics[{#1,Line[{#,0}&/@{-1,0}]},AspectRatio->0.1,ImageSize->18],"\!\("<>#2<>"\)"}&@@@{{Red,"N\_th/N"},{Blue,"N\_0/N"}},Spacings->{0.2,0.1}];
Manipulate[
Block[{n=1,z,Ts,zs,ns,Tc},
z[T_,z0_:0.9999]:=Which[Δ==1&&n/T>5,1-SetPrecision[Exp[-n/T],20],
Δ>1&&n>=nTz[Δ,T,1],0.99999999,True,z/.FindRoot[nTz[Δ,T,z]==n,{z,z0}]];
Ts=Range[0.01,1.6,0.01];
zs=Reverse@Rest@FoldList[z[#2,#1]&,0.9999,Reverse@Ts];
ns=nTz[Δ,#1,#2]&@@@Thread[{Ts,zs}];
Tc=If[Δ>1,(n/(Gamma[Δ+1]Zeta[Δ]))^(1/Δ),0];
ListLinePlot[{Thread[{Ts,ns}],Thread[{Ts,1-ns}]},PlotStyle->{Red,Blue},GridLines->{{Tc},None},Epilog->{Inset[legend,Scaled[{0.95,0.9}],{1,1}]},FrameLabel->{"\!\(T/ϵ\_0\)","Fraction"}]
],
{{Δ,3/2,Style["\!\( Δ\)","InlineFormula"]},{1/2,1,3/2,2,3}},
Paneled->False,SaveDefinitions->True]]

The condensate fraction in the BEC phase

Derive Eq. ().

(Subscript[N, th])=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^Δ)(1).

N=Δ V (((((Subscript[T, c])/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^Δ)(1).

((Subscript[N, th])/N)=((Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^Δ)(1))/(Δ V (((((Subscript[T, c])/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^Δ)(1)))=((((T/(Subscript[T, c]))))^Δ).

((Subscript[N, 0])/N)=((N-(Subscript[N, th]))/N)
=1-((Subscript[N, th])/N)
=1-((((T/(Subscript[T, c]))))^Δ).

The BEC transition temperature (Subscript[T, c]) is set by the temperature at which z→1,

where ζ(s):=(Subscript[Li, s])(1) is the Riemann zeta function, defined by

ζ(s)=(∑_(k=1)^∞)(1/(k^s)).

Derive Eq. ().

N=(Subscript[N, th])=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^Δ)(1)
=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)Γ(Δ)(Subscript[Li, Δ])(1)
=Γ(Δ+1)ζ(Δ) V ((((T/(Subscript[ϵ, 0]))))^Δ).

(Subscript[T, c])=(((Subscript[ϵ, 0])((N/(Γ(Δ+1)ζ(Δ) V))))^(1/Δ)).

For non-relativistic particles,

(Subscript[T, c])=((2π)/m)((((1/(ζ(D/2) ))(N/V)))^(2/D)).

To increase the BEC transition temperature, one can

either reduce the boson mass m,

or increase the boson density N/V.

Dimension dependence of (Subscript[T, c]) (following Eq. ())

Manipulate[Block[{Tc},
Tc[d_]:=If[d<=2,0,(n/Zeta[d/2])^(2/d)];Plot[{Tc[d],0,1.8},{d,0,6},PlotStyle->{LightGray,None,None},Filling->{2->{{1},Lighter[Blue,0.9]},3->{{1},Lighter[Red,0.9]}},Epilog->{{Black,Table[Point[{d,Tc[d]}],{d,0,6}]},{Blue,Text["BEC",{4.5,Tc[4]/2.5}]},{Red,Text["Normal",({4,#+0.5}/2+{-#,0.3}/Sqrt[1+#^2])&@Tc[3]]}},PlotLabel->StringTemplate["\!\(N/V\)=`1`"][n],PlotRange->{{0,6},{0,1.8}},PlotRangePadding->{Scaled[.025],Scaled[.04]},FrameLabel->{"\!\( D\)","\!\(T\_c/(2π/m)\)"},ImageSize->200]],
{{n,1,"\!\(N/V\)"},0.2,5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula",SaveDefinitions->True]

BEC transition only happens in 3D and higher dimensions (for non-relativistic bosons).

In general, BEC transition only happens in bosonic systems with Δ=D/ζ>1. According to Eq. (), (Subscript[T, c])→0 as Δ->(SubPlus[1]) due to the divergence of the Riemann zeta function ζ(Δ→1)→∞.

Ultracold Atoms

Examples are found in bosonic atoms, such as Rb, Na, Li, of a number N∼(10^4)→(10^7) below the temperature (Subscript[T, c])∼(10^(-7))K. [2001 Nobel Prize]

Release of Atoms: Initially, the atoms are confined in an optical trap. The trap (trapping laser) is tuned off, allowing atoms to expand freely.

Free Expansion: The atoms expand for a set duration. During this time, their positions evolve according to their initial momenta: atoms with larger momenta will travel further.

Imaging: Following the expansion period, the spatial distribution of the atoms is captured using imaging techniques. This image reflects the momentum distribution of the atoms.

Column[{Show[,ImageSize->240],Show[,ImageSize->160]}]

K. B. Davis, M. -O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Bose-Einstein Condensation in a Gas of Sodium Atoms, Phys. Rev. Lett. 75, 3969 (1995)

Aggregating within each momentum shell 4π (k^2)dk gives the distribution for the momentum magnitude k≡|k|:

(Subscript[N, k])=((4π (k^2))/(exp((1/T)(((k^2)/(2m))-μ))-1)).

However, this only accounts for the thermal boson. One should also include the potential contribution from (Subscript[N, 0]) bosons condensed in the ground state.

(Subscript[N, k])=((4π (k^2))/(exp((1/T)(((k^2)/(2m))-μ))-1))+(Subscript[N, 0])δ(k).

DynamicModule[{ℐ,nTz},
ℐ[Δ_,z_]:=Gamma[Δ+1]PolyLog[Δ+1,z];
nTz[Δ_,T_,z_]:=Re[Δ T^Δ ℐ[Δ-1,z]];
Manipulate[Block[{Δ=3/2,m,n=1,z0=0.99999,z,f0,f1,Tc,T},
m=(2π) Gamma[Δ+1]^(1/Δ);
Tc=(n/(Gamma[Δ+1]Zeta[Δ]))^(1/Δ);
T=Tc x;
z=If[n>=nTz[Δ,T,1],1,z/.FindRoot[nTz[Δ,T,z]==n,{z,z0}]];
f1[k_]:=4π k^2/(Exp[((k^2)/(2m))/T]/z-1)/(2π)^3;
f0[k_]:=(n-nTz[Δ,T,z])Block[{σ=0.1},Exp[-k^2/(2σ^2)]/(Sqrt[2π]σ)];
Plot[f0[k]+f1[k],{k,0,8},Filling->Axis,PlotRange->{0,1},PlotRangePadding->{Scaled[0.02],Scaled[0.03]},FrameLabel->{"\!\( k\)","\!\(N\_k\)"},PlotLabel->StringTemplate["\!\(T/T\_c=`1`\)"][T/Tc],ImageSize->200]
],{{x,0.9,"\!\( T/T\_c\)"},0.01,2,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula",SaveDefinitions->True]
]

Condensation Transition

Near the BEC transition, the fugacity is close to one. So we expand (Subscript[ℐ, +]^Δ)(z) around z->1.

Block[{ℐ,plt},
ℐ[η_,s_,z_?NumericQ]:=Gamma[s]PolyLog[s,η z]/η;
Block[{r=3},
Plot[Evaluate[Tanh[ℐ[#,2,r ArcTanh[ζ]]/r]&/@{1,-1}],{ζ,0,1},Prolog->{{{Opacity[0.3],Blue,Disk[Tanh[{1,ℐ[1,2,1]}/r],0.1]},{Blue,Text["BEC",Tanh[{1,ℐ[1,2,1]}/r]+0.1 ReIm[Exp[I π/4]],{-1.1,0.3},{1,1}]}},{Dashed,AbsoluteThickness[1],Line[{#,#}&/@{0,1}]}},
Epilog->{Translate[{Blue,Text["Bosons",{0,0},{0,-0.6}],Text[Style["(\!\(η=+1\))",12],{0,0},{0,0.6}]},{0.18,0.75}],Translate[{Red,Text["Fermions",{0,0},{0,-0.6}],Text[Style["(\!\(η=-1\))",12],{0,0},{0,0.6}]},{0.75,0.2}]},
PlotRange->{{0,1},{0,1}},PlotRangePadding->0.05,PlotRangeClipping->False,GridLines->({#,#}&@Flatten@{{0},Tanh[Range[1,6]/r],{1}}),FrameLabel->{"\!\( z\)","\!\(ℐ\_η\%s(z)\)"},FrameTicks->({{#,Automatic},{#,Automatic}}&@{{0,0},{Tanh[1/r],1},{1,∞}}),ImagePadding->{{Automatic,10},{Automatic,10}},AspectRatio->1,ImageSize->160]]]

Verify Eq. ().

General case:

Block[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
Assuming[{ϵ>0,s>2},(First[#]+FullSimplify@Series[Rest[#],{ϵ,0,-1}])&@Series[ℐ[1,s,1-ϵ],{ϵ,0,1}]]/.{ϵ->HoldForm[1-z]}]

Special cases:

Block[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
Assuming[{0<z<1},FullSimplify@Series[ℐ[1,1,z],{z,1,1}]]]

Block[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
Assuming[{0<z<1},FullSimplify@Series[ℐ[1,2,z],{z,1,1}]]]

The expansion Eq. () enables us to compute the thermal boson number (Subscript[N, th]) [see Eq. ()] and the energy E [see Eq. ()]:

(Subscript[N, th])=((3(Sqrt[π])ζ(3/2))/4)V ((((T/(Subscript[ϵ, 0]))))^(3/2))(1-((2(Sqrt[π]))/(ζ(3/2)))(Sqrt[1-z])+…),
E=((9(Sqrt[π])ζ(5/2))/8)T V ((((T/(Subscript[ϵ, 0]))))^(3/2))(1-((ζ(3/2))/(ζ(5/2)))(1-z)+…).

Derive Eq. ().

(Subscript[ℐ, +]^(3/2))(z) around z->1

(Subscript[ℐ, +]^(3/2))(z)=Γ(3/2)(ζ(3/2)-ζ(1/2)(1-z)+…)+((π (((1-z))^(1/2)))/(sin(3π/2)))+…
=((Sqrt[π])/2)(ζ(3/2)-ζ(1/2)(1-z)+…)-π (((1-z))^(1/2))+…,

given that Γ(3/2)=(Sqrt[π])/2. Note that (1-z) is a higher order term compared to (((1-z))^(1/2)) and thus can be dropped. So, to the leading order, we have

(Subscript[ℐ, +]^(3/2))(z)=((Sqrt[π])/2)ζ(3/2)-π (Sqrt[1-z])+…
=((Sqrt[π])/2)ζ(3/2)(1-((2(Sqrt[π]))/(ζ(3/2))) (Sqrt[1-z])+…).

Therefore,

(Subscript[N, th])=(3/2)V ((((T/(Subscript[ϵ, 0]))))^(3/2))(Subscript[ℐ, +]^(3/2))(z)
=((3(Sqrt[π])ζ(3/2))/4)V ((((T/(Subscript[ϵ, 0]))))^(3/2))(1-((2(Sqrt[π]))/(ζ(3/2)))(Sqrt[1-z])+…).

(Subscript[ℐ, +]^(5/2))(z) around z->1

(Subscript[ℐ, +]^(5/2))(z)=Γ(5/2)(ζ(5/2)-ζ(3/2)(1-z)+…)+((π (((1-z))^(3/2)))/(sin(5π/2)))+…
=((3(Sqrt[π]))/4)(ζ(5/2)-ζ(3/2)(1-z)+…)+π (((1-z))^(3/2))+…,

given that Γ(5/3)=3(Sqrt[π])/4. Note that (((1-z))^(3/2)) is a higher order term compared to (1-z) and thus can be dropped. So, to the leading order, we have

(Subscript[ℐ, +]^(5/2))(z)=((3(Sqrt[π]))/4)(ζ(5/2)-ζ(3/2)(1-z)+…)
=((3(Sqrt[π])ζ(5/2))/4)(1-((ζ(3/2))/(ζ(5/2)))(1-z)+…).

Therefore,

E=(3/2)T V ((((T/(Subscript[ϵ, 0]))))^(3/2))(Subscript[ℐ, +]^(5/2))(z)
=((9(Sqrt[π])ζ(5/2))/8)T V ((((T/(Subscript[ϵ, 0]))))^(3/2))(1-((ζ(3/2))/(ζ(5/2)))(1-z)+…).

These expressions in Eq. () are not very friendly. Let us simplify them by introducing:

The BEC transition temperature (Subscript[T, c]) (adapt from Eq. () or Eq. ()),

(Subscript[T, c])=(((Subscript[ϵ, 0])((4/(3(Sqrt[π])ζ(3/2)))(N/V)))^(2/3)).

The energy (Subscript[E, c]) of the system at T=(Subscript[T, c]),

(Subscript[E, c])=((9(Sqrt[π])ζ(5/2))/8)(Subscript[T, c]) V (((((Subscript[T, c])/(Subscript[ϵ, 0]))))^(3/2))=(3/(2 B))N (Subscript[T, c]).

Two numerical constants A and B,

A=((2(Sqrt[π]))/(ζ(3/2)))≈1.35697,
B=((ζ(3/2))/(ζ(5/2)))≈1.9473724663169563`.

Then Eq. () can be organized into a relatively neat form:

((Subscript[N, th])/N)=((((T/(Subscript[T, c]))))^(3/2))(1-A(Sqrt[1-z])),

E=(Subscript[E, c]) ((((T/(Subscript[T, c]))))^(5/2))(1-B (1-z)).

The behavior of the fugacity z near the transition can be analyzed as follows:

When T>(Subscript[T, c]), (Subscript[N, th])=N (thermal excitations can host all bosons), then z can be determined by solving Eq. (),

1=((((T/(Subscript[T, c]))))^(3/2))(1-A(Sqrt[1-z])),

whose solution is

z=1-(1/(A^2))(((1-(((((Subscript[T, c])/T)))^(3/2))))^2).

When T<(Subscript[T, c]), it is no longer possible for (Subscript[N, th]) to match N. The remaining (Subscript[N, 0])=N-(Subscript[N, th]) bosons will condense. In this case, z=(Subscript[N, 0])/((Subscript[N, 0])+1)→1 will become very close to 1 (given that (Subscript[N, 0]) is macroscopically large).

Block[{A=2Sqrt[π]/Zeta[3/2]},
Plot[1-Max[0,1-(1/T)^(3/2)]^2/A^2,{T,0.8,1.2},PlotStyle->Black,GridLines->{{1},{1}},FrameLabel->{"\!\(T/T\_c\)","\!\( z\)"},FrameTicks->{{Table[{z,If[Mod[z/0.01,1]==0.,NumberForm[z,{3,2}],""],If[Mod[z/0.01,1]==0.,{0.02,0},{0.015,0}]},{z,0.9,1,0.002}],Automatic},{Automatic,Automatic}},ImageSize->180]]

Note that the piecewise nature of z (as a function of T) is the fundamental reason that leads to singularities in various thermodynamic functions at T=(Subscript[T, c]), as will be seen later.

Plugging the solution Eq. () to Eq. (), we obtain the energy E as a function of temperature T (with fixed N and V, packed in (Subscript[T, c]) and (Subscript[E, c]))

Block[{A=2Sqrt[π]/Zeta[3/2],B=Zeta[3/2]/Zeta[5/2],Ε},
Ε[T_]:=T^(5/2)(1-B/A^2 Max[0,1-(1/T)^(3/2)]^2);
Plot[Ε[T],{T,0,1.6},PlotStyle->Black,GridLines->{{1},{1}},FrameLabel->{"\!\(T/T\_c\)","\!\( E/E\_c\)"},ImageSize->180]]

1st order derivative: (isochoric) heat capacity -- the amount of heat (per particle) required to change the temperature.

(Subscript[C, V])=(Subscript[(((∂E)/(∂T))), N,V]).

Show that the heat capacity (Subscript[C, V]) in Eq. () is related to the specific heat (Subscript[c, V]) in Eq. () by (Subscript[C, V])=(Subscript[c, V])N.

(Subscript[C, V])=(Subscript[(((∂E)/(∂T))), N,V])
=(Subscript[(((∂E)/(∂S))), N,V])(Subscript[(((∂S)/(∂T))), N,V])
=(Subscript[T(((∂S)/(∂T))), N,V]).

(Subscript[c, V])=(T/N)(Subscript[(((∂S)/(∂T))), N,V]),

Block[{A=2Sqrt[π]/Zeta[3/2],B=Zeta[3/2]/Zeta[5/2],Ε},
Ε[T_]:=T^(5/2)(1-B/A^2 Max[0,1-(1/T)^(3/2)]^2);
Plot[Ε'[T],{T,0,1.6},PlotStyle->Black,GridLines->{{1},{5/2}},FrameLabel->{"\!\(T/T\_c\)","\!\( C\_V=∂\_T E\)"},ImageSize->180]]

The heat capacity exhibits a cusp at T=(Subscript[T, c]), which is rooted in the similar cusp behavior in (∂z/∂T)∼max(0,1-((((Subscript[T, c])/T))^(3/2))).

2nd order derivative: starts to exhibit a discontinuity at T=(Subscript[T, c]).

Block[{A=2Sqrt[π]/Zeta[3/2],B=Zeta[3/2]/Zeta[5/2],Ε},
Ε[T_]:=T^(5/2)(1-B/A^2 Max[0,1-(1/T)^(3/2)]^2);
Plot[Ε''[T],{T,0,1.6},Exclusions->{1},PlotStyle->Black,GridLines->{{1},None},FrameLabel->{"\!\(T/T\_c\)","\!\( ∂\_T\%2E\)"},ImageSize->180]]

1st order phase transitions:

Boiling of water

Melting of ice

Structural transitions in some materials

2nd order (or continuous) phase transitions:

BEC transitions

Superconducting transitions

Magnetic transitions in some materials

Liquid crystal transitions

3rd order phase transitions:
(They are less common or more theoretical.)

Certain topological quantum phase transitions (e.g. mass generation transition of 2D Dirac fermions)

Infinite order phase transitions:

Berezinskii-Kosterlitz-Thouless (BKT) transitions (e.g. superfluid transition in 2D) [2016 Nobel Prize]

Superfluid Helium

Graphics[{{Blue,Line[{2.2,0.1}+# ReIm[Exp[0.54I π]]&/@{-0.06,2.8}]},{Green,BSplineCurve[{{0,0},{2.2,0},{4.2,0.1},{4.9,0.18},{5.15,0.28}}],Point[{5.15,0.28}]},{Red,BSplineCurve[{{0,2.5},{0.9,2.5},{1.3,2.6},{1.6,2.7},{2,3.0},{2.1,3.5}}]},Translate[{Text["He\[Hyphen]II",{0,0},{0,-0.6}],Text[Style["(superfluid)",12],{0,0},{0,0.6}]},{1.02,0.7}],Translate[{Text["He\[Hyphen]I",{0,0},{0,-0.6}],Text[Style["(normal fluid)",12],{0,0},{0,0.6}]},{4,1.2}],Text["He solid",{0.9,3.2}],{Blue,Text["λ\[Hyphen]line",{2,1.6},{0,-0.7},ReIm[Exp[-0.46I π]]]},{Green,Text["vaporizing",{4,0.4},{0,0},{1,0.15}]},{Red,Text["melting",{2.15,3.3},{-1,0}]}},GridLines->{{},{}},PlotRange->{{0,6},{0,3.5}},PlotRangePadding->{Scaled[0.03],Scaled[0.04]},Frame->True,FrameLabel->{"\!\( T\) [K]","\!\( P\) [MPa]"},AspectRatio->0.7,ImageSize->220]

A classical liquid (normal fluid) phase, called He-I: finite viscosity, limited heat conductivity, staying in its container normally.

A quantum liquid (superfluid) phase, called He-II: zero viscosity, nearly infinite thermal conductivity, can spontaneously climb walls and escape its container.

Block[{Tc=2.2,CV,plt},
CV[T_]:=-0.8 Log[Abs[T-Tc]/(1+Abs[T-Tc])](1-0.5UnitStep[T-Tc])(T^3+0.02T^5);
plt[{Tmin_,Tmax_},ref_,scl_,side_]:=Plot[CV[ref Tc+scl T],{T,(Tmin-ref Tc)/scl,(Tmax-ref Tc)/scl},Exclusions->None,PlotRange->{0,70},Axes->None,PlotRangePadding->{Scaled[0.05],3},PlotStyle->Black,GridLines->{{Switch[ref,0,Tc,1,0]},None},FrameLabel->{Style[Switch[ref,0,"\!\( T\)",1,"\!\(T-T\_c\)"]<>Switch[scl,1," [K]",0.1,"\n[\!\(10\^-1\)K]",0.01,"\n[\!\(10\^-2\)K]"],LineSpacing->{0,1}],"\!\(C\_V\) [J/mol·K]"},ImagePadding->{{38side+2, 2}, {65, 2}},AspectRatio->1.5,ImageSize->{Automatic,180}];
Row@{plt[{0,4},0,1,1],plt[Tc+0.2{-1,1},1,0.1,0],plt[Tc+0.02{-1,1},1,0.01,0]}]

BEC refers to the condensate of non-interacting bosons (gas),

Superfluid is the condensate of (repulsively) interacting bosons (liquid).

Black-Body Radiation

Radiations Carried by Bosons

Electromagnetic radiation (light) is carried by photons (spin-1).

Acoustic radiation (sound) is carried by phonons (spin-0).

Spin wave radiation is carried by magnons (spin-1).

Gravitational radiation is carried by gravitons (spin-2).

*Dark matter radiation could hypothetically be carried by axions (spin-0).

The physical laws of black-body radiation applies to all different types of radiations, although we will mainly focus on photons gas in the following discussions.

Grid[{{"System","T","Dominant radiation"},{"Cosmic microwave background (CMB)","2.7 K","Microwave"},
{"Cold interstellar dusts","∼10K","Terahertz"},
{"Night hemisphere of the Earth","300 K","Infrared"},
{"Electric stove/oven","500K","Infrared"},
{"Incandescent light","2500K","Infrared to Visible"},{"The Sun","6000K","Visible"},{"Black hole accretion disk","∼\!\(10\^6\)K","X\[Hyphen]ray"}},Dividers->{None,2->True},ItemStyle->"Text"]

Zero Chemical Potential for Radiations

This is because a heat bath could absorb one high-energy boson and emit two lower-energy bosons (causing boson number to +1), or vice versa (causing boson number to -1).

For examples, atoms can serve as heat bath to absorb and emit photons.

Row[{Block[{photon,atom},
photon[col_,θ_,r_]:={AbsoluteThickness[1],col,Translate[Rotate[Line@Table[{x,0.2Cos[π x]Sin[6π x]},{x,-0.5,0.5,0.01}],θ,{}],r]};
SeedRandom[1];
Graphics[{Line[4.8{{-1,-1},{-1,1},{1,1},{1,-1},{-1,-1}}],Block[{n=10},Disk[#1,0.2]&/@RandomReal[{-4.5,4.5},{n,2}]],Block[{n=20},photon@@@Thread@{RandomChoice[{Red,Orange,Purple},n],RandomReal[{-π,π},n],RandomReal[{-4,4},{n,2}]}]},PlotRange->{{-5,5},{-5,5}},ImageSize->{Automatic,100}]],
DynamicModule[{wiggle,arrow,bubble,ah,bg,dir=1},
ah[r_,θ_,w_:3,h_:5]:=Translate[Rotate[Line[Offset/@{{-h,-w},{0,0},{-h,w}}],θ,{0,0}],r];
wiggle[{p1_,p2_}]:=Block[{t,n},t=p2-p1;n=Normalize[{{0,1},{-1,0}}.t];
{Line@Table[p1+x t+0.1n Sin[12π x] ,{x,0,1,0.01}],Switch[dir,1,ah[p2,ArcTan@@t],-1,ah[p1,π+ArcTan@@t]]}];
arrow[{p1_,p2_}]:=Block[{t},t=p2-p1;
{Line[{p1,p2}],Switch[dir,1,ah[p2,ArcTan@@t],-1,ah[p1,π+ArcTan@@t]]}];
EventHandler[Dynamic[Graphics[{
{AbsoluteThickness[1],Gray,arrow[{-0.5,#}&/@{0,8}],arrow[{0.5,#}&/@{8,2}],arrow[{0.5,#}&/@{2,0}]},Translate[Line[{#,0}&/@{-1,1}],{0,#}]&/@{0,2,8},Translate[{Purple,wiggle[{#,0}&/@{-5,-0.5}],Text["\!\(ℏ ω\_1\)",{-5,0},{-1,-1}]},{0,4}],
Translate[{Orange,wiggle[{#,0}&/@{0.5,5}],Text["\!\(ℏ ω\_2\)",{5,0},{1,-1}]},{0,5}],Translate[{Red,wiggle[{#,0}&/@{0.5,5}],Text["\!\(ℏ ω\_3\)",{5,0},{1,-1}]},{0,1}]},ImageSize->{Automatic,100}],TrackedSymbols:>{dir}],{"MouseClicked":>(dir=-dir)}]]},Spacer[5]]

In such processes,

while energy remains conserved => ⟨E⟩=E,

the particle number does not => ⟨N⟩=??.

If the particle number is not a conserved quantity, we should not require a definitive expectation value for the number of particles in the system. Therefore, the grand canonical ensemble is not a good starting point for our discussion.

Our problem: how to infer the chemical potential  μ with no knowledge about the particle number N.

According to Eq. (), the entropy S is related to chemical potential μ (or the fugacity z=(e^(μ/T))) by

S=V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ+1)(Subscript[ℐ, +]^Δ)(z)-Δ (Subscript[ℐ, +]^(Δ-1))(z) log z).

Given V and T fixed, S is a monotonically increasing function of z∈(0,1) for any choice of Δ>0.

Manipulate[
Plot[Block[{ℐ},ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;(Δ+1)ℐ[1,Δ+1,z]-Δ Log[z]ℐ[1,Δ,z]],{z,0,1},PlotRange->All,PlotStyle->Black,FrameLabel->{"\!\( z\)","\!\( S\)"},ImagePadding->{{45,5},{40,5}},ImageSize->{Automatic,140}],{{Δ,3/2},{1/2,1,3/2,2,3}},LabelStyle->"InlineFormula",Paneled->False,SaveDefinitions->True]

Prove that S is a monotonically increasing function of z∈(0,1).

((∂S)/(∂z))=V ((((T/(Subscript[ϵ, 0]))))^Δ)(∂_z)((Δ+1)(Subscript[ℐ, +]^(Δ+1))(z)-Δ (Subscript[ℐ, +]^Δ)(z)log z)
(Overscript[=, [1]])V ((((T/(Subscript[ϵ, 0]))))^Δ)((((Δ+1)Δ)/z) (Subscript[ℐ, +]^Δ)(z)-(Δ/z) (Subscript[ℐ, +]^Δ)(z)-((Δ (Δ-1))/z) (Subscript[ℐ, +]^(Δ-1))(z)log z)
=V ((((T/(Subscript[ϵ, 0]))))^Δ)(Δ/z)(Δ (Subscript[ℐ, +]^Δ)(z)- (Δ-1) (Subscript[ℐ, +]^(Δ-1))(z)log z)
=V ((((T/(Subscript[ϵ, 0]))))^Δ)(Δ/z)(Δ Γ(Δ)(Subscript[Li, Δ])(z)- (Δ-1) Γ(Δ-1)(Subscript[Li, Δ-1])(z)log z)
(Overscript[=, [2]])V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ Γ(Δ))/z)(Δ (Subscript[Li, Δ])(z)-(Subscript[Li, Δ-1])(z)log z)
(Overscript[=, [3]])V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ Γ(Δ))/z)(Δ (∑_(k=1)^∞)((z^k)/(k^Δ))-(∑_(k=1)^∞)((z^k)/(k^(Δ-1)))log z)
=V ((((T/(Subscript[ϵ, 0]))))^Δ)Δ Γ(Δ)(∑_(k=1)^∞)((Δ-k log z)/(k^Δ))(z^(k-1)).

The chemical potential is always zero for radiation-carrying bosons.

In fact, there is no need to introduce the chemical potential  μ as a Lagrangian multiplier in the first place, if the statistical ensemble is not constrained by the particle number N.

Spectral Energy Density

(E/V)=∫ρ(ω)dω.

In quantum mechanics, due to the wave-particle duality, the radiation frequency ω is related to the boson single-particle energy ϵ by

(We may set ℏ=1).

Knowing that radiations are carried by bosons (η=+1) of zero chemical potential ( μ=0), Eq. () implies that the spectral energy density for the black-body radiation should take the following form

ρ(ω)=((ω g(ω))/((e^( β ω))-1))
=((Δ )/(Subscript[ϵ, 0]^Δ))((ω^Δ)/((e^( β ω))-1)).

Derive Eq. () and Eq. ().

E=V∫dϵ ((ϵ g(ϵ))/((e^( β ϵ))-1))
=V∫dω ((ω g(ω))/((e^( β ω))-1)).

ρ(ω)=((ω g(ω))/((e^( β ω))-1))
=(ω/((e^( β ω))-1))((Δ (ω^(Δ-1)))/(Subscript[ϵ, 0]^Δ))
=((Δ )/(Subscript[ϵ, 0]^Δ))((ω^Δ)/((e^( β ω))-1)).

(1/(Subscript[ϵ, 0]^Δ))=((Subscript[g, s])/(Γ(D/2+1)))((((1/(2 (π^(1/2))c))))^D),

ρ(ω)=(((Subscript[g, s]) D)/(Γ(D/2+1)))((((1/(2 (π^(1/2))c))))^D) ((ω^D)/((e^( β ω))-1))
=(((Subscript[g, s]) D)/((2^D)(π^(D/2))Γ(D/2+1)))((((ω/c))^D)/((e^( β ω))-1)).

ρ(ω)=(1/((e^( β ω))-1))({{{, ({{((2ω)/(π c)), D=1,}, {((ω^2)/(π (c^2))), D=2,}, {((ω^3)/(((π )^2)(c^3))), D=3.}})}})

Table[2d/Gamma[d/2+1](ω/(2Sqrt[π] c))^d,{d,1,3}]

Manipulate[Plot[ω^3/(Exp[ω/T]-1),{ω,0,10},PlotRange->All,PlotStyle->Black,GridLines->{{2.821 T},None},FrameLabel->{"\!\(ω\)","\!\(ρ(ω)\)"},ImagePadding->{{60,5},{40,5}},ImageSize->{Automatic,140}],{{T,1},0.1,2,ImageSize->Small},LabelStyle->"InlineFormula",Paneled->False,SaveDefinitions->True]

Derive Eq. ().

(∂_ω)ρ(ω)=0.

(∂_ω)ρ(ω)=-(1/((π^2)(c^3)))(((ω^2)(3+(e^(β ω))(β ω-3)))/((((e^(β ω))-1))^2)),

D[ω^3/(Exp[β ω]-1),ω]//Together//Simplify

3+(e^x)(x-3)=0,

FindRoot[(3-3 (E^x)+(E^x) x)==0,{x,2.1}]

β ω=x≈2.821.

Thermodynamic Properties

Pressure

P=κ (T^(Δ+1)).

Entropy

S=(Δ+1) κ V (T^Δ).

Particle number

N=((ζ(Δ))/(ζ(Δ+1)))κ V (T^Δ).

Energy

E=Δ κ V (T^(Δ+1)).

Free energy (F=E-T S=𝔉+μ N)

F=𝔉=-κ V (T^(Δ+1)),

κ=(1/(Subscript[ϵ, 0]^Δ))Γ(Δ+1)ζ(Δ+1).

Derive Eq. (-).

P=T  ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^(Δ+1))(1),
S=(Δ+1)V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^(Δ+1))(1),
N=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^Δ)(1),
E=Δ T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^(Δ+1))(1),
F=𝔉=-T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, +]^(Δ+1))(1).

κ:=(1/(Subscript[ϵ, 0]^Δ))(Subscript[ℐ, +]^(Δ+1))(1)=(1/(Subscript[ϵ, 0]^Δ))Γ(Δ+1)ζ(Δ+1),

P=κ (T^(Δ+1)) ,
S=(Δ+1) κ V (T^Δ),
N=Δ κ V (T^Δ)(((Subscript[ℐ, +]^Δ)(1))/((Subscript[ℐ, +]^(Δ+1))(1))),
E=Δ κ V (T^(Δ+1)),
F=𝔉=-κ V (T^(Δ+1)).

N=Δ κ V (T^Δ)(((Subscript[ℐ, +]^Δ)(1))/((Subscript[ℐ, +]^(Δ+1))(1)))
=Δ κ V (T^Δ)((Γ(Δ)ζ(Δ))/(Γ(Δ+1)ζ(Δ+1)))
=((ζ(Δ))/(ζ(Δ+1)))κ V (T^Δ).

Equation of state:

which deviated from ideal classical gas behavior (P V=N T) by a factor ζ(Δ+1)/ζ(Δ). For relativistic bosons in three dimensional space, we have Δ=3, and ζ(4)/ζ(3)≈0.9004, a bit smaller than 1.

In terms of energy E,

we can express

Thermodynamics processes:

Grid[{{"Process","Relation","Notes"},{"Isobaric","\!\(P=const.\)\n\!\(T=const.\)","as \!\( P\) and \!\( T\) are locked together by \!\(P = κ T\^\(Δ+1\)\)"},
{"Isothermal",SpanFromAbove,SpanFromAbove},
{"Isochoric","\!\(V=const.\)",""},
{"Adiabatic","\!\(P V\^γ=const.\)","with \!\(γ=1+1/Δ\)"}},Alignment->{Left,Center},ItemStyle->"InlineFormula",Dividers->{{False,{True},False},{False,{True},False}}]

Specific Heat

Isochoric (constant-volume) specific heat

It exhibits a (Subscript[c, V])∼(T^Δ) behavior, a important feature of radiation gas.

Isobaric (constant-pressure) specific heat can not be defined as making P constant prevents us to change T.

Degenerated Fermi Gas

Fermi Energy

⟨(Subscript[n, α])⟩=(1/((e^(β((Subscript[ϵ, α])-μ)))+1))(Overscript[→, T->0])({{{, ({{1, for (Subscript[ϵ, α])<μ,}, {0, for (Subscript[ϵ, α])>μ.}})}})

Manipulate[Block[{β=1/T},Plot[1/(Exp[β ξ]+1),{ξ,-5,5},PlotStyle->Red,PlotRange->{0,1},Filling->Axis,PlotRangePadding->Scaled[.05],
PlotLabel->StringTemplate[" \!\(T =`1`\)"][T],FrameLabel->{"\!\(ϵ\_α-μ\)","\!\(⟨n\_α⟩\)"},ImageSize->180]],{{T,0.01," \!\(T=1/β\)"},0.01,5,ImageSize->Small},Paneled->False,LabelStyle->"InlineFormula"]

Each fermion added to the system will just fill up the available single-particle states from low energy to high energy.

The energy of the last filled state is called the Fermi energy, defined by the chemical potential at T=0,

N=V∫dϵ ((g(ϵ))/((e^( β (ϵ-μ)))+1))
(Overscript[=, T→0])V(∫_0^(Subscript[ϵ, F]))dϵ g(ϵ).

N=V (((((Subscript[ϵ, F])/(Subscript[ϵ, 0]))))^Δ),

While (Subscript[ϵ, F]) is the highest energy that a fermion can take, the average energy E/N carried by a fermion is a fraction of (Subscript[ϵ, F]), given by

Derive Eq. ().

E=V∫dϵ ((ϵ g(ϵ))/((e^( β (ϵ-μ)))+1))
(Overscript[=, T→0])V(∫_0^(Subscript[ϵ, F]))dϵ ϵ g(ϵ)

E=V(∫_0^(Subscript[ϵ, F]))dϵ ϵ ×((Δ (ϵ^(Δ-1)))/(Subscript[ϵ, 0]^Δ))
=((Δ V)/(Subscript[ϵ, 0]^Δ))(∫_0^(Subscript[ϵ, F]))dϵ (ϵ^Δ)
=((Δ V)/(Subscript[ϵ, 0]^Δ))((Subscript[ϵ, F]^(Δ+1))/(Δ+1))
=(Δ/(Δ+1))((V(((Subscript[ϵ, F])/(Subscript[ϵ, 0]))))^Δ)(Subscript[ϵ, F])
=(Δ/(Δ+1))N (Subscript[ϵ, F])

It goes as E=(3/5)N (Subscript[ϵ, F]) for Δ=(3/2) (for non-relativistic fermions in three dimensional space).

Given E=Δ P V [recall Eq. ()], the pressure P of the degenerated Fermi gas (at T=0) is

Amazingly, P does not vanish at zero temperature for Fermi gas (unlike the ideal classical gas P V=N T).

The finite pressure of Fermi gas at low temperature is called the degeneracy pressure.

It is a consequence of Pauli exclusion principle and is important in understanding the stability of white dwarf stars and neutron stars in astrophysics.

(Subscript[T, F])=((Subscript[ϵ, F])/(Subscript[k, B])).

Under the setting of (Subscript[k, B])=1, temperature and energy are essentially the same unit, with

1.60218×(10^(-19)) J=1 eV=11604.5 K≃(10^4) K.

(Subscript[T, F]) can vary by orders of magnitudes in different systems:

Grid[{{"System","\!\(ϵ\_F\)","\!\(T\_F\)","\!\( T\)"},
{"Ultracold \!\(\*TemplateBox[<|"boxes" -> FormBox[
RowBox[{
SuperscriptBox["\[Null]", "6"], "", 
StyleBox["Li", FontSlant -> "Plain"]}], TraditionalForm], "errors" -> {}, "input" -> "{}^6\\!\\mathrm{Li}", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\) atoms","∼ 0.1 neV","∼ \!\(10\^-6\)K","∼ \!\(10\^-7\)K"},{"Superfluid \!\(\*TemplateBox[<|"boxes" -> FormBox[
RowBox[{
SuperscriptBox["\[Null]", "3"], "", 
StyleBox["He", FontSlant -> "Plain"]}], TraditionalForm], "errors" -> {}, "input" -> "{}^{3}\\!\\mathrm{He}", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\) atoms","∼ 0.1 meV","∼ 1 K","∼ \!\(10\^-3\)K"},{"Electrons in metal","∼ 1 eV","∼ \!\(10\^4\) K","300 K"},{"Electrons in white dwarf stars", "∼ 0.1  MeV","∼ \!\(10\^9\)K","∼ \!\(10\^\(6∼8\)\)K "},{"Nucleons in atomic nucleus","∼ 10 MeV","∼ \!\(10\^11\)K","--"},{"Neutrons in neutron starts","∼ 0.1 GeV","∼ \!\(10\^12\)K","∼ \!\(10\^\(6∼11\)\)K"}},Dividers->{None,{2->True}},ItemStyle->{Automatic,{"InlineFormula",{"Text"}}}]

Compact Stars

Block[{colfun,cls,seq,xs},
colfun=Blend[{Lighter[LightBlue,0.4],Lighter[Cyan,0.8],White,Yellow,Orange,Red},Interpolation[{{32,0},{17,0.15},{8.4,0.3},{6.4,0.4},{5.3,0.5},{4,0.85},{2.7,1}},InterpolationOrder->1][#]]&;
cls=<|"O"->32,"B"->17,"A"->8.4,"F"->6.4,"G"->5.3,"K"->4,"M"->2.7|>;
xs=<|"I"->{0.1,0.4,0.7,0.9},"III"->{0,0.3,0.6,1},"V"->Range[0,1,0.05],"W"->{0,0.06,0.12,0.2,0.3,0.45,0.65,0.8,1}|>;
seq=({-Log[#1],#2}&@@@#&)/@<|"I"->{{20,4.5},{10,4},{6,3.8}},"III"->{{5,1.8},{2.8,2.4}},"V"->{{25,3.2},{7.5,1},{6,0},{5,-1},{3,-4}},"W"->{{30,-2},{6,-4}}|>;
Graphics[{White,Table[Block[{f},f=BSplineFunction[seq[c]];
Table[{colfun[Exp[-f[x][[1]]]],AbsolutePointSize[6.2+f[x][[2]]],Point[f[x]]},{x,xs[c]}]],{c,Keys[seq]}],{GrayLevel[0.9],Text[Style[#,12,FontFamily->"CMU Bright",LineSpacing->{0,10}],{-Log[#2[[1]]],#2[[2]]},##3]&@@@{{"White dwarfs",{16,-2.6},{0,-1},{1,-0.4}},{"Main sequence",{8,0},{0,1},{1,-1.1}},{"Giants",{3.8,1.8},{0,1},{1,0.4}},{"Supergiants",{9,3.7},{0,1},{1,-0.2}}}}},Background->GrayLevel[0.15],PlotRange->{{-Log[35],-Log[2]},{-5,5}},Frame->True,FrameLabel->{{"Luminosity (\!\(Sun=1\))","Absolute magnitude"},{"Temperature [kK]","Spectral color"}},FrameStyle->GrayLevel[0.85],FrameTicksStyle->GrayLevel[0.85],FrameTicks->{{{#,If[#==0,"1",StringTemplate["\!\(10\^`1`\)"][#]]}&/@Range[-4,4,2],{#+4.83-5,Block[{M=5(1-#/2)},If[M>0,"+",""]<>ToString[M]]}&/@Range[-4,4,2]},{{-Log[#],#}&/@{30,20,10,7,5,4,3,2},
Table[{-Log[cls[lb]],Style[lb,colfun[cls[lb]]],{0,0}},{lb,Keys@cls}]}},AspectRatio->1.2,ImageSize->240]]

Main sequence stars: the majority of stars, including the Sun, in stable phase of stellar evolution, fusing H into He in their cores.

White dwarf stars: remnants of low to medium mass stars, like our Sun, no nuclear fusion anymore and are supported against gravitational collapse by electron degeneracy pressure.

Neutron stars: remnants of supernova from more massive stars,  gravity is strong enough to push electrons into nucleus,

p+(SuperMinus[e])→n+(Subscript[ν, e]),

leaving neutrons alone to fight against gravitational collapse by neutron degeneracy pressure. (Too faint to be visible on the HR diagram.)

Black holes: extremely dense remnants formed from the gravitational collapse of very massive stars post-supernova. Nothing can sustained against gravity, even light can not escape.

E=(Subscript[E, gas])+(Subscript[E, grav])

The energy of the degenerated electron gas

(Subscript[E, gas])=(Δ/(Δ+1))N (Subscript[ϵ, F])∝((M^(1+1/Δ))/(R^(3/Δ))).

Verify Eq. ().

(Subscript[E, gas])=(Δ/(Δ+1))N (Subscript[ϵ, F])
=(Δ/(Δ+1))N (Subscript[ϵ, 0]) ((((N/V)))^(1/Δ))
=((Δ (Subscript[ϵ, 0]))/(Δ+1)) ((N^(1+1/Δ))/(V^(1/Δ)))

The number N of electrons in the star equals to the number protons in the star (due to electric neutrality), which is then proportional to the mass of the star, so N∝M.

The volume V of the star is determined by the radius R via V=((4π)/3)(R^3), so V∝(R^3).

(Subscript[E, gas])∝((N^(1+1/Δ))/(V^(1/Δ)))∝((M^(1+1/Δ))/(R^(3/Δ)))

The gravitational energy

(Subscript[E, grav])∝-((M^2)/R).

Derive Eq. ().

E=a((M^(1+1/Δ))/(R^(3/Δ)))-b((M^2)/R),

((∂E)/(∂R))=-((3a)/Δ)((M^(1+1/Δ))/(R^(1+3/Δ)))+b((M^2)/(R^2))=0,

((3a)/Δ)((M^(1+1/Δ))/(R^(1+3/Δ)))=b((M^2)/(R^2))
=>((M^(1/Δ-1))/(R^(3/Δ-1)))=((Δ b)/(3 a))=const.

R∼(M^((1/Δ-1)/(3/Δ-1)))=(M^((1-Δ)/(3-Δ))).

For non-relativistic electrons: Δ=3/2, Eq. () predicts

R∼(M^(-1/3)),

meaning that the star will shrink as it gets more massive.

Because as M increases, gravity compresses the star more, forcing electrons to get closer, generating a higher degeneracy pressure (Fermi energy) to balance the stronger gravity, thus resulting in a smaller star radius R.

However, this trend does not sustain for ever. As (Subscript[ϵ, F]) grows to the point that it becomes comparable to the electron rest mass energy (Subscript[m, e])(c^2), relativity effects sets in.

For relativistic electrons: Δ=3, Eq. () predicts

R∼(M^∞),

a strange expression ... but it just means M∼(R^0)=const., i.e. M will top out when R reduces to the point that electrons become relativistic.

In the relativistic regime, the Fermi gas energy ((Subscript[E, gas])∼(R^(-1))) does not increase as rapidly as it does in the non-relativistic regime ((Subscript[E, gas])∼(R^(-2))). Consequently, there comes a point where the pressure can no longer counterbalance the gravity for a more massive star.

(Subscript[M, limit])≈1.44(Subscript[M, ⊙]).

a type Ia supernova (the “standard candle” in astronomy),

or a more compact object, like a neutron star or a black hole.

Block[{curve},
curve[M0_,γ_,{m0_,m1_}]:=Block[{α,Edeg,Egra,RM},
α=(8 M0/(5Sqrt[5]))^(2/3)/γ;
Edeg[M_,R_]:=Block[{x,y,B},y=γ M^(1/3)/R;x=Sqrt[1+y^2]-1;B=3x(x+1)/(4x+5);α M B];
Egra[M_,R_]:=-(3/5)M^2/R;
RM[M_]:=If[M>=M0,0,If[M>0.99M0,4γ/Sqrt[13]M0^(1/3)((M0/M)^(2/3)-1)^(1/2),Block[{lnR,sol},sol=FindMinimum[Edeg[M,Exp[lnR]]+Egra[M,Exp[lnR]],{lnR,Log[4γ/Sqrt[13]M0^(1/3)]}];Exp[lnR]/.Last@sol]]];Off[FindMinimum::lstol];
Table[Block[{M=M0 (Tanh[m]+1)/2},{Log[10,RM[M]],Log[10,M]}],{m,m0,m1,0.1}]];
Graphics[{{Lighter[Red,0.5],AbsoluteThickness[1],Arrow@BSplineCurve[{{0,0},{0.2,0},{0.8,0},{0,-0.3},{-1.8,-0.3}}]},Translate[Block[{r=0.08,a=0.03},{Red,Disk[{0,0},r],Rotate[Translate[Triangle[a{{0,1},{3,0},{0,-1}}],{1.2r,0}],2π #/12,{0,0}]&/@Range[12],Text[Style["the Sun",12],{0,r+2a},{0.35,-0.6}],Text[Style["death",12],{-0.5,-0.3},{0,0.8}]}],{0,0}],Translate[{{Blue,Line@#,Dashed,Line[{Last@#,Last@#-{1.2,0}}]}&@curve[1.44,1,{-2.35,3}],{Blue,Text[Style["\!\(1.44M\_⊙\)",12],{-1.3,0.45}],Text[Style["\!\(R∼M\^\(-1/3\)\)",10],{0.4,-1.2},{0,-1},{1,-2.7}],Text[Style["White\ndwarf",LineSpacing->{0,12}],{-0.7,-0.5}]}},{-1.88,0}],Translate[{{Green,Line@#}&@curve[5.76,1,{-1.6,0.3}],{Green,Text[Style["Neutron\nstar",LineSpacing->{0,12}],{0.6,-0.4},{0,1},{1,-2.5}]}},{-4.84,0}],{Black,Line[{#,#+5.372}&/@{-5.4,-3.45}],Text["Black hole",{#,#+5.37}&@-4.4,{0,-0.5},{1,1}]}},PlotRange->{{-5.5,0.5},{-2,2}},GridLines->{{0},{0}},Frame->True,FrameLabel->{"\!\(R/R\_⊙\)","\!\(M/M\_⊙\)"},FrameTicks->{{{#,If[#==0,"1",StringTemplate["\!\(10\^`1`\)"][#]]}&/@Range[-1,2],{#,""}&/@Range[-1,2]},{{#,If[#==0,"1",StringTemplate["\!\(10\^`1`\)"][#]]}&/@Range[-5,1],{#,""}&/@Range[-5,1]}},AspectRatio->Automatic,ImageSize->240]]

first expand into a red giant (swallowing the Earth),

then blow away its outer layers,

leaving behind a white dwarf (cooling over time).

Low-Temperature Expansion

z=(e^(μ/T))≃(e^((Subscript[ϵ, F])/T))→+∞ .

Block[{ℐ,plt},
ℐ[η_,s_,z_?NumericQ]:=Gamma[s]PolyLog[s,η z]/η;
Block[{r=3},
Plot[Evaluate[Tanh[ℐ[#,2,r ArcTanh[ζ]]/r]&/@{1,-1}],{ζ,0,1},Prolog->{{{Opacity[0.3],Red,Disk[{1,1},0.1]},{Red,Text["Low\[Hyphen]\!\( T\)",{0.9,1},{1.1,0.6}]}},{Dashed,AbsoluteThickness[1],Line[{#,#}&/@{0,1}]}},
Epilog->{Translate[{Blue,Text["Bosons",{0,0},{0,-0.6}],Text[Style["(\!\(η=+1\))",12],{0,0},{0,0.6}]},{0.18,0.75}],Translate[{Red,Text["Fermions",{0,0},{0,-0.6}],Text[Style["(\!\(η=-1\))",12],{0,0},{0,0.6}]},{0.75,0.2}]},
PlotRange->{{0,1},{0,1}},PlotRangePadding->0.05,PlotRangeClipping->False,GridLines->({#,#}&@Flatten@{{0},Tanh[Range[1,6]/r],{1}}),FrameLabel->{"\!\( z\)","\!\(ℐ\_η\%s(z)\)"},FrameTicks->({{#,Automatic},{#,Automatic}}&@{{0,0},{Tanh[1/r],1},{1,∞}}),ImagePadding->{{Automatic,10},{Automatic,10}},AspectRatio->1,ImageSize->160]]]

Verify Eq. ().

(Subscript[ℐ, -]^s)(z):=(∫_0^∞)dx((x^(s-1))/((z^(-1))(e^x)+1))
=(∫_0^∞)dx((x^(s-1))/((e^(x-log z))+1)).

Θ(x)=({{{, ({{1, x>0,}, {0, x<0.}})}})

(1/((e^(x-log z))+1))=Θ(log z-x)+Φ(x-log z),

the leading step function part: Θ(log z-x),

and the residual part: Φ(x-log z), with the function Φ defined by

Φ(y):=(1/((e^y)+1))-Θ(-y).

Plot[1/(Exp[y]+1)-HeavisideTheta[-y],{y,-10,10},PlotStyle->Black,PlotRange->All,FrameLabel->{"\!\( y\)","\!\(Φ(y)\)"},ImageSize->180]

It is an odd function, i.e. Φ(-y)=-Φ(y).

It converges to 0 exponentially fast as |y| increases from 0.

(Subscript[ℐ, -]^s)(z)=(∫_0^∞)dx((x^(s-1))/((e^(x-log z))+1))
=(∫_0^∞)dx (x^(s-1))Θ(log z-x)+(∫_0^∞)dx (x^(s-1))Φ(x-log z)

The first integral simply reduces to a integration of x up to log z

(∫_0^∞)dx (x^(s-1))Θ(log z-x)
=(∫_0^(log z))dx (x^(s-1))
=((((log z))^s)/s).

The second integral can be analyzed by a change of variable. Define y=x-log z (such that x=y+log z),

(∫_0^∞)dx (x^(s-1))Φ(x-log z)
=(∫_(-log z)^∞)dy (((y+log z))^(s-1))Φ(y).

In the z->+∞ limit, we have -log z→-∞, so the integration lower bound -log z can be replaced by -∞. This is a controlled approximation as Φ(y) decays to zero exponentially with increasing |y|, so the integral is not sensitive to the lower bound far from 0. Therefore,

=(∫_(-log z)^∞)dy (((y+log z))^(s-1))Φ(y)
≃(∫_(-∞)^∞)dy (((y+log z))^(s-1))Φ(y)
=(∫_(-∞)^∞)dy (∑_(n=0)^∞)(({{s-1}, {n}}))(((y^n)(log z))^(s-1-n))Φ(y)
=(∑_(n=0)^∞)(({{s-1}, {n}}))(((log z))^(s-1-n))(∫_(-∞)^∞)dy (y^n)Φ(y),

where (({{a}, {b}}))=((Γ(a+1))/(Γ(b+1)Γ(a-b+1))) denotes the binomial function. Since Φ(y) is an odd function of y, the integral (∫_(-∞)^∞)dy (y^n)ε(y) will be non-vanishing only if n is odd. So we can set n=2m-1 and only focus on the odd order terms in the expansion.

=(∑_(n=0)^∞)(({{s-1}, {n}}))(((log z))^(s-1-n))(∫_(-∞)^∞)dy (y^n)Φ(y)
=(∑_(m=1)^∞)(({{s-1}, {2m-1}}))(((log z))^(s-2m))(∫_(-∞)^∞)dy (y^(2m-1))Φ(y)
=2(∑_(m=1)^∞)(({{s-1}, {2m-1}}))(((log z))^(s-2m))(∫_0^∞)dy (y^(2m-1))Φ(y)
=2(∑_(m=1)^∞)(({{s-1}, {2m-1}}))(((log z))^(s-2m))(∫_0^∞)dy (y^(2m-1))((1/((e^y)+1))-Θ(-y))
=2(∑_(m=1)^∞)(({{s-1}, {2m-1}}))(((log z))^(s-2m))(∫_0^∞)dy ((y^(2m-1))/((e^y)+1))
=2(∑_(m=1)^∞)(({{s-1}, {2m-1}}))(((log z))^(s-2m))(Subscript[ℐ, -]^(2m))(1).

We have use the definition (Subscript[ℐ, -]^s)(1)=(∫_0^∞)dy((y^(s-1))/((e^y)+1)) in the last step.

(Subscript[ℐ, -]^s)(z)=((((log z))^s)/s)+2(∑_(m=1)^∞)(({{s-1}, {2m-1}}))(((log z))^(s-2m))(Subscript[ℐ, -]^(2m))(1)
=(((log z))^s)((1/s)+2(∑_(m=1)^∞)(({{s-1}, {2m-1}}))(((Subscript[ℐ, -]^(2m))(1))/(((log z))^(2m)))).

(Subscript[ℐ, -]^2)(1)=((π^2)/12),(Subscript[ℐ, -]^4)(1)=((7(π^4))/120),…,

Block[{ℐ},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
ℐ[-1,#,1]&/@{2,4}]

(Subscript[ℐ, -]^s)(z)
=(((log z))^s)((1/s)+2(({{s-1}, {1}}))(((Subscript[ℐ, -]^2)(1))/(((log z))^2))+2(({{s-1}, {3}}))(((Subscript[ℐ, -]^4)(1))/(((log z))^4))+…)
=(((log z))^s)((1/s)+(((π^2)(s-1))/(6(((log z))^2)))+((7(π^4)(s-1)(s-2)(s-3))/(360(((log z))^4)))+…)
=((((log z))^s)/s)(1+(((π^2)s (s-1))/(6(((log z))^2)))+((7(π^4)s (s-1)(s-2)(s-3))/(360(((log z))^4)))+…).

The expansion Eq. () enables us to compute the fermion number N [see Eq. ()] and the energy E [see Eq. ()] at low temperature,

N=V ((((μ/(Subscript[ϵ, 0]))))^Δ)(1+(((π^2)Δ (Δ-1))/(6 )) ((((T/μ)))^2)+…),

E=(Δ/(Δ+1)) μ V ((((μ/(Subscript[ϵ, 0]))))^Δ)(1+(((π^2)Δ (Δ+1))/(6 ))((((T/μ)))^2)+…).

Derive Eq. (), Eq. ().

N=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, -]^Δ)(z)
=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)((((log z))^Δ)/Δ)(1+(((π^2)Δ (Δ-1))/(6 (((log z))^2))) +…)
=V ((((T/(Subscript[ϵ, 0]))))^Δ)((((μ/T)))^Δ)(1+(((π^2)Δ (Δ-1))/(6 )) ((((T/μ)))^2)+…)
=V ((((μ/(Subscript[ϵ, 0]))))^Δ)(1+(((π^2)Δ (Δ-1))/(6 )) ((((T/μ)))^2)+…).

E=Δ T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, -]^(Δ+1))(z)
=Δ T V ((((T/(Subscript[ϵ, 0]))))^Δ)((((log z))^(Δ+1))/(Δ+1))(1+(((π^2) (Δ+1)Δ)/(6 (((log z))^2)))+…)
=(Δ/(Δ+1)) T V ((((T/(Subscript[ϵ, 0]))))^Δ)((((μ/T)))^(Δ+1))(1+(((π^2)Δ (Δ+1))/(6 ))((((T/μ)))^2)+…)
=(Δ/(Δ+1)) μ V ((((μ/(Subscript[ϵ, 0]))))^Δ)(1+(((π^2)Δ (Δ+1))/(6 ))((((T/μ)))^2)+…).

For fixed fermion number N, we can solve Eq. () for the chemical potential  μ,

μ=(Subscript[ϵ, F])(1-(((π^2)  (Δ-1) )/6)((((T/(Subscript[ϵ, F]))))^2)+…),

where (Subscript[ϵ, F])=(Subscript[ϵ, 0]) (((N/V))^(1/Δ)) is the Fermi energy.

Derive Eq. ().

(N/V)=(((((Subscript[ϵ, F])/(Subscript[ϵ, 0]))))^Δ).

(((((Subscript[ϵ, F])/(Subscript[ϵ, 0]))))^Δ)=(N/V)=((((μ/(Subscript[ϵ, 0]))))^Δ)(1+(((π^2)Δ (Δ-1))/(6 )) ((((T/μ)))^2)+…),

μ=(Subscript[ϵ, F]) (((1+(((π^2)Δ (Δ-1))/(6 )) ((((T/μ)))^2)+…))^(-1/Δ))
(Overscript[≃, [1]])(Subscript[ϵ, F]) (((1+(((π^2)Δ (Δ-1))/(6 )) ((((T/(Subscript[ϵ, F]))))^2)+…))^(-1/Δ))
(Overscript[≃, [2]])(Subscript[ϵ, F]) (1-(1/Δ)(((π^2)Δ (Δ-1))/(6 )) ((((T/(Subscript[ϵ, F]))))^2)+…)
=(Subscript[ϵ, F]) (1-(((π^2) (Δ-1))/(6 )) ((((T/(Subscript[ϵ, F]))))^2)+…).

Block[{n=2,SeriesFactor,ℐ,Nμ,Eμ,μT,μ},
SeriesFactor[series_]:=Times@@(HoldForm/@{First@series[[3]],SeriesData[series[[1]],series[[2]],Factor[series[[3]]/First@series[[3]]],series[[4]],series[[5]],series[[6]]]});
ℐ[Δ_]:=(μ/T)^(Δ+1)(1/(Δ+1)+SeriesData[T,0,Flatten@Table[{0,FunctionExpand[2Gamma[Δ+1]/Gamma[Δ-2m+2](-PolyLog[2m,-1])]/μ^(2m)},{m,1,n}],1,2n+1,1]);
Nμ=PowerExpand[Δ Ν (T/ϵF)^Δ ℐ[Δ-1]];
Eμ=PowerExpand[Δ Ν  T(T/ϵF)^Δ ℐ[Δ]];
μT=ϵF+Sum[a[m]T^(2m),{m,1,n}]+O[T]^(2n+1);
μ=μT/.Factor[First@Solve[(Nμ/.μ->μT)==Ν,a/@Range[n]]];
SeriesFactor[μ]]

We can see  μ=(Subscript[ϵ, F]) at T=0, as expected.

As T increases,  μ will deviates from (Subscript[ϵ, F]) slowly, in quadratic order of T/(Subscript[ϵ, F]).

Substitute the solution of  μ in Eq. () to the expression of E in Eq. (), we obtain E as a function of T at fixed N and V.

Derive Eq. ().

E=(Δ/(Δ+1)) μ V ((((μ/(Subscript[ϵ, 0]))))^Δ)(1+(((π^2)Δ (Δ+1))/(6 ))((((T/μ)))^2)+…)
=(Δ/(Δ+1)) (Subscript[ϵ, F]) V (((((Subscript[ϵ, F])/(Subscript[ϵ, 0]))))^Δ)(((1-(((π^2)  (Δ-1) )/6)((((T/(Subscript[ϵ, F]))))^2)+…))^(Δ+1))(1+(((π^2)Δ (Δ+1))/(6 ))((((T/(Subscript[ϵ, F]))))^2)+…)
=(Δ/(Δ+1)) N (Subscript[ϵ, F]) (1-(((π^2)  (Δ+1) (Δ-1) )/6)((((T/(Subscript[ϵ, F]))))^2)+…)(1+(((π^2)Δ (Δ+1))/(6 ))((((T/(Subscript[ϵ, F]))))^2)+…)
=(Δ/(Δ+1)) N (Subscript[ϵ, F]) (1-(((π^2)  (Δ+1) (Δ-1) )/6)((((T/(Subscript[ϵ, F]))))^2)+(((π^2)Δ (Δ+1))/(6 ))((((T/(Subscript[ϵ, F]))))^2)+…)
=(Δ/(Δ+1)) N (Subscript[ϵ, F]) (1+(((π^2)  (Δ+1)  )/6)((((T/(Subscript[ϵ, F]))))^2)+…)

Block[{n=2,SeriesFactor,ℐ,Nμ,Eμ,μT,μ},
SeriesFactor[series_]:=Times@@(HoldForm/@{First@series[[3]],SeriesData[series[[1]],series[[2]],Factor[series[[3]]/First@series[[3]]],series[[4]],series[[5]],series[[6]]]});
ℐ[Δ_]:=(μ/T)^(Δ+1)(1/(Δ+1)+SeriesData[T,0,Flatten@Table[{0,FunctionExpand[2Gamma[Δ+1]/Gamma[Δ-2m+2](-PolyLog[2m,-1])]/μ^(2m)},{m,1,n}],1,2n+1,1]);
Nμ=PowerExpand[Δ Ν (T/ϵF)^Δ ℐ[Δ-1]];
Eμ=PowerExpand[Δ Ν  T(T/ϵF)^Δ ℐ[Δ]];
μT=ϵF+Sum[a[m]T^(2m),{m,1,n}]+O[T]^(2n+1);
μ=μT/.Factor[First@Solve[(Nμ/.μ->μT)==Ν,a/@Range[n]]];
SeriesFactor[Eμ]]

Specific Heat

(Subscript[c, V])=(1/N)(Subscript[(((∂E)/(∂T))), N,V]),

Derive Eq. ().

(Subscript[c, V])=(1/N)(Subscript[(((∂E)/(∂T))), N,V])
=(Δ/(Δ+1))  (Subscript[ϵ, F]) ((((π^2)  (Δ+1)  )/6)(∂_T)((((T/(Subscript[ϵ, F]))))^2)+…)
=(((π^2)Δ )/6)(Subscript[ϵ, F])((2T)/(Subscript[ϵ, F]^2))+…
=(((π^2)Δ )/3)((T/(Subscript[ϵ, F])))+….

Block[{n=2,SeriesFactor,ℐ,Nμ,Eμ,cV,μT,μ},
SeriesFactor[series_]:=Times@@(HoldForm/@{First@series[[3]],SeriesData[series[[1]],series[[2]],Factor[series[[3]]/First@series[[3]]],series[[4]],series[[5]],series[[6]]]});
ℐ[Δ_]:=(μ/T)^(Δ+1)(1/(Δ+1)+SeriesData[T,0,Flatten@Table[{0,FunctionExpand[2Gamma[Δ+1]/Gamma[Δ-2m+2](-PolyLog[2m,-1])]/μ^(2m)},{m,1,n}],1,2n+1,1]);
Nμ=PowerExpand[Δ Ν (T/ϵF)^Δ ℐ[Δ-1]];
Eμ=PowerExpand[Δ Ν  T(T/ϵF)^Δ ℐ[Δ]];
μT=ϵF+Sum[a[m]T^(2m),{m,1,n}]+O[T]^(2n+1);
μ=μT/.Factor[First@Solve[(Nμ/.μ->μT)==Ν,a/@Range[n]]];
cV=D[Eμ,T]/Ν;
SeriesFactor[cV]]

The specific heat of a Fermi gas is linear in T at low temperature.

Recall Eq. () that (Subscript[c, V])=Δ is the result for a classical gas. Compare to that, the specific heat of a degenerated Fermi gas is smaller by a factor of T/(Subscript[ϵ, F]) (or T/(Subscript[T, F])).

The specific heat essentially counts the number of modes (per particle) that can be thermally excited to participate in energy exchange.

In classical gas, every particle can participate in energy exchange with the heat bath.

In degenerated Fermi gas, only a T/(Subscript[ϵ, F]) fraction of the fermions near the Fermi surface are excited and can contribute to the energy exchange. The remaining fermions deep in the Fermi sea are effective frozen.

Manipulate[Block[{η=-1,ℐ,PTz,sTz,nTz,znT},
ℐ[η_,s_,z_]:=Gamma[s]PolyLog[s,η z]/η;
PTz[T_,z_]:=T^(Δ+1)ℐ[η,Δ+1,z];
sTz[T_,z_]:=T^Δ((Δ+1)ℐ[η,Δ+1,z]-Δ Log[z]ℐ[η,Δ,z]);
nTz[T_,z_]:=Δ T^Δ ℐ[η,Δ,z];
znT[n_,T_,z0_:1]:=Block[{z},
If[η==1&&n>=nTz[T,1],1,z/.FindRoot[nTz[T,z]==n,{z,z0}]]];
Block[{Tmax=3,Ts,zs,ss,sT,data},
Ts=Tmax E Exp[-1/Range[0.16,1,0.02]];
Off[General::munfl,FindRoot::lstol];
zs=Re@Reverse@Rest@FoldList[znT[n,#2,#1]&,0,Reverse@Ts];
ss=Re[sTz[#1,#2]/n]&@@@Thread@{Ts,zs};
sT=Interpolation[Prepend[Thread@{Ts,ss},{0,0}],InterpolationOrder->1];
data=Prepend[Table[{T,Chop[T sT'[T]]},{T,Mean/@Thread@{Most@Ts,Rest@Ts}}],{0,0}];
ListLinePlot[data,PlotStyle->Black,PlotRange->{0,Max[Max[data[[All,2]]],Δ]},PlotRangePadding->{Scaled[0.02],Scaled[0.05]},FrameLabel->{"\!\(T/ϵ\_0\)","\!\(c\_V\)"},Prolog->{{Gray,AbsoluteThickness[1],Dashed,Line[{#,Δ}&/@{0,Tmax}],Text["\!\(Δ\)",{0,Δ},{-1.2,0.2}]},Translate[{Gray,AbsoluteThickness[1],Dashed,HalfLine[{0,#}&/@{0,Δ}],Text["\!\(ϵ\_F/ϵ\_0\)",{0,0},{-1.2,-0.6}]},{n^(1/Δ),0}],{Gray,AbsoluteThickness[1],Dashed,Line[{3#/(π^2 Δ)n^(1/Δ),#}&/@{0,Δ}]}},PlotLabel->StringTemplate["\!\(N/V=`1`\)"][n],ImageSize->200]
]],
{{Δ,3/2,Style["\!\(Δ\)","InlineFormula"]},{1/2,1,3/2,2,3}},{{n,1,Style["\!\(N/V\)","InlineFormula"]},0.01,2,ImageSize->Small},Paneled->False,SaveDefinitions->True]

Susceptibilities

(Subscript[χ, N])=(1/N)(Subscript[(((∂N)/(∂μ))), V,T]).

Derive Eq. ().

N=V ((((μ/(Subscript[ϵ, 0]))))^Δ)+…,

(Subscript[χ, N])=(1/N)(Subscript[(((∂N)/(∂μ))), V,T])
=(Subscript[(((∂log N)/(∂μ))), V,T])
=(∂/(∂μ))(Δ log μ+const.+…)
=(Δ/μ)+….

(Subscript[χ, N])=(Δ/(Subscript[ϵ, F]))+….

Block[{n=2,SeriesFactor,ℐ,Nμ,Eμ,χN,μT,μ},
SeriesFactor[series_]:=Times@@(HoldForm/@{First@series[[3]],SeriesData[series[[1]],series[[2]],Factor[series[[3]]/First@series[[3]]],series[[4]],series[[5]],series[[6]]]});
ℐ[Δ_]:=(μ/T)^(Δ+1)(1/(Δ+1)+SeriesData[T,0,Flatten@Table[{0,FunctionExpand[2Gamma[Δ+1]/Gamma[Δ-2m+2](-PolyLog[2m,-1])]/μ^(2m)},{m,1,n}],1,2n+1,1]);
Nμ=PowerExpand[Δ Ν (T/ϵF)^Δ ℐ[Δ-1]];
Eμ=PowerExpand[Δ Ν  T(T/ϵF)^Δ ℐ[Δ]];
χN=D[Nμ,μ]/Ν;
μT=ϵF+Sum[a[m]T^(2m),{m,1,n}]+O[T]^(2n+1);
μ=μT/.Factor[First@Solve[(Nμ/.μ->μT)==Ν,a/@Range[n]]];
SeriesFactor[χN]
]

Charge susceptibility: the rate that the charge density Q/N changes with the electric potential ℰ.

(Subscript[χ, Q])=(1/N)(Subscript[(((∂Q)/(∂ℰ))), V,T]).

The total charge of the electron gas is

Q=-e N.

Varying in electric potential dℰ will cause a change in the electron chemical potential by

dμ=-e dℰ.

Therefore,  (Subscript[χ, Q]) is related to  (Subscript[χ, N]) by

Spin susceptibility: the rate that the magnetization density M/N changes with the magnetic field B.

(Subscript[χ, M])=(1/N)(Subscript[(((∂M)/(∂B))), V,T]).

The total magnetization of the electron gas is

M=(Subscript[μ, B])((Subscript[N, ↑])-(Subscript[N, ↓])),

where (Subscript[N, ↑]) and (Subscript[N, ↓]) are respectively the number of up-spin and down-spin electrons in a magnetic field.

Varying the magnetic field dB will cause a change in the chemical potential for up-spin and down-spin electrons,

d(Subscript[μ, ↑])=(Subscript[μ, B])dB,
d(Subscript[μ, ↓])=-(Subscript[μ, B])dB.

Therefore,  (Subscript[χ, M]) is related to  (Subscript[χ, N]) by

This is also known as the Pauli susceptibility, which (partially) explains the paramagnetism of many metals.

(Subscript[c, V])=(((π^2) Δ  )/3)((T/(Subscript[ϵ, F])))+…,
  (Subscript[χ, Q])=((Δ (e^2))/(Subscript[ϵ, F]))+…,
  (Subscript[χ, M])=((Δ (Subscript[μ, B]^2))/(Subscript[ϵ, F]))+…,

For non-interacting electrons, the Wilson ratio is expected to be unity. This prediction holds for most metals, where electron-electron interactions are well screened and thus negligible.

However, deviations from unity in the Wilson ratio are indicative of significant electron-electron interactions within the materials. -- A defining feature of strongly correlated electron systems, including many unconventional high-temperature superconductors.

Summary

Partition Function is All You Need

General Idea

Given the energy function E(x) for some random variable x.

Under maximal entropy estimation, the equilibrium probability distribution of x is given by the exponential form, with β=1/T,

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

The partition function Z is not just a normalization constant

its also encodes all the thermodynamic properties of the statistical ensemble.

Or equivalently, in terms of the free energy F

Derivative of log Z with respect to β, energy expectation value is obtained,

((∂ log Z)/(∂β))=-⟨E(x)⟩.

Derive Eq. ().

((∂ log Z)/(∂β))=(1/Z)((∂Z)/(∂β))
=(1/Z)(∂/(∂β))(∑_x)(e^(-β E(x)))
=(1/Z)(∑_x)(∂/(∂β))((e^(-β E(x))))
=-(1/Z)(∑_x)E(x)(e^(-β E(x)))
=-(∑_x)E(x)p(x)
=-⟨E(x)⟩.

Suppose the energy function (Subscript[E, λ])(x) further depends on a Lagrangian multiplier λ for a function f(x) of x

derivative of log Z with respect to λ gives the expectation value of f(x)

((∂ log Z)/(∂λ))=-β⟨f(x)⟩.

Derive Eq. ().

((∂ log Z)/(∂λ))=(1/Z)((∂Z)/(∂λ))
=(1/Z)(∂/(∂λ))(∑_x)(e^(-β (Subscript[E, λ])(x)))
=(1/Z)(∑_x)(∂/(∂λ))((e^(-β (Subscript[E, λ])(x))))
=-(β/Z)(∑_x)(((∂(Subscript[E, λ])(x))/(∂λ)))(e^(-β (Subscript[E, λ])(x)))
=-β(∑_x)f(x)p(x)
=-β⟨f(x)⟩.

Or more directly as

Cumulant Expansion

To be more systematic, the first few orders of derivatives of log Z are given by (denote ∂/∂λ as (∂_λ))

(∂_λ) log Z=(((∂_λ)Z)/Z),
(∂_λ^2) log Z=(((∂_λ^2)Z)/Z)-((((((∂_λ)Z)/Z)))^2),
(∂_λ^3) log Z=(((∂_λ^3)Z)/Z)-3(((∂_λ^2)Z)/Z)(((∂_λ)Z)/Z)+2((((((∂_λ)Z)/Z)))^2),
….

Verify Eq. ().

NestList[Dt,Log[Z],3]

The partition function derivatives are related to different moments of f(x),

(((∂_λ^n)Z)/Z)=⟨(((-β f(x)))^n)⟩,

for n=1,2,3,….

Prove Eq. ().

(((∂_λ^n)Z)/Z)=(1/Z)(∑_x)(∂_λ^n)(e^(-β (Subscript[E, λ])(x)))
=(1/Z)(∑_x)(∂_λ^n)(e^(-β (E(x)+λ f(x))))
=(1/Z)(∑_x)(∂_λ^n)((e^(-β E(x)))(e^( λ (-β f(x)))))
=(1/Z)(∑_x)(e^(-β E(x)))(∂_λ^n)(e^( λ (-β f(x))))

(∂_λ^n)(e^( λ (-β f(x))))=(((-β f(x)))^n)(e^( λ (-β f(x)))),

(((∂_λ^n)Z)/Z)=(1/Z)(∑_x)(((e^(-β E(x)))(-β f(x)))^n)(e^( λ (-β f(x))))
=(1/Z)(∑_x)(((-β f(x)))^n)(e^(-β (Subscript[E, λ])(x)))
=(∑_x)(((-β f(x)))^n)p(x)
=⟨(((-β f(x)))^n)⟩.

Substitute Eq. () into Eq. (), the derivatives of log Z gives us different cumulants of f(x),

They corresponds to the mean, variance, skewness of f(x).

The particle number fluctuation is quantified by the variance of N(n),

var N:=⟨((N(n))^2)⟩-((⟨N(n)⟩)^2)
=(((-β))^(-2))(∂_μ^2) log Z
=-T(((∂^2)F)/(∂(μ^2)))

The particle number susceptibility is defined as Eq. ()

(Subscript[χ, N])=(1/N)(((∂N)/(∂μ)))=-(1/N)(((∂^2)F)/(∂(μ^2))),

given that N=-(∂_μ)F.

Mid-Term Sample Problem

Both particle and anti-particle occupations cost the same amount of energy (Subscript[ϵ, α]) on the (Ket[α]) state, while the empty occupation does not cost energy.

Anti-particle counts as a “negative” particle in terms of particle number, as it can annihilate with the particle.

E(n)=(∑_α)(Subscript[ϵ, α])(Subscript[n, α]^2),
N(n)=(∑_α)(Subscript[n, α]).

g(ϵ)=((Δ (ϵ^(Δ-1)))/(Subscript[ϵ, 0]^Δ)),

Nondimentionalization

Rule 1: the argument of a transcendental function must be dimensionless.

Examples of transcendental functions: exp, log, trigonometric functions (cos, sin, tan, …, cosh, sinh, tanh, …).

This is because things like exp(5Joule) does not make sense.

Rule 2: perform a change of the integration variable to make it dimensionless.

I=(∫_0^∞)dϵ (ϵ^s)f(β ϵ),

I=(1/(β^(s+1)))(∫_0^∞)dx (x^s)f(x).

I∝(1/(β^(s+1))).

Quantum Gas

Grand Canonical Ensemble

The microstate of a quantum gas is described by a collection of occupation numbers

n:={(Subscript[n, α])|α=1,2,…}.

Energy and particle number for a specific collection n

E(n)=(∑_α)(Subscript[ϵ, α])(Subscript[n, α]), N(n)=(∑_α)(Subscript[n, α]).

Probability distribution

p(n)=(1/ℨ)(e^(-β(E(n)-μ N(n)))),

and the partition function

ℨ=(∑_n)(e^(-β(E(n)-μ N(n)))).

Bosons and Fermions

Introduce the fermion sign

η=({{{, ({{+1, bosons,}, {-1, fermions.}})}})

Partition function ℨ=(∏_α)(Subscript[ℨ, α]), with

log (Subscript[ℨ, α])=-η log(1-η (e^(-β ((Subscript[ϵ, α])-μ)))).

Average occupation number

⟨(Subscript[n, α])⟩=(1/((e^( β ((Subscript[ϵ, α])-μ)))-η)).

Density of State

Definition

∀f: (∑_α)f((Subscript[ϵ, α]))=V∫dϵ g(ϵ)f(ϵ).

General form

g(ϵ)=((Δ (ϵ^(Δ-1)))/(Subscript[ϵ, 0]^Δ)),

with

Δ=(D/ζ),

D - dimension of space,

ζ - dynamical exponent ((Subscript[ϵ, k])∼(k^ζ)).

Thermodynamics

Introduce fugacity

z=(e^(β μ))=(e^(μ/T)).

Free energy 𝔉=-T log ℨ,

𝔉=-T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z),

Pressure

P=-(Subscript[(((∂𝔉)/(∂V))), T,μ])=T  ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z).

Entropy

S=-(Subscript[(((∂𝔉)/(∂T))), V,μ])=V ((((T/(Subscript[ϵ, 0]))))^Δ)((Δ+1)(Subscript[ℐ, η]^(Δ+1))(z)-((Δ μ)/T)(Subscript[ℐ, η]^Δ)(z)).

Particle number

N=-(Subscript[(((∂𝔉)/(∂μ))), V,T])=Δ V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^Δ)(z).

Energy

E=𝔉+T S+μ N
=Δ T V ((((T/(Subscript[ϵ, 0]))))^Δ)(Subscript[ℐ, η]^(Δ+1))(z)
=Δ P V.

Equation of state

P V=(((Subscript[ℐ, η]^(Δ+1))(z))/(Δ(Subscript[ℐ, η]^Δ)(z)))N T.

Higher order properties

Specific heats

(Subscript[c, V])=(T/N)(Subscript[(((∂S)/(∂T))), N,V]),
(Subscript[c, P])=(T/N)(Subscript[(((∂S)/(∂T))), N,P]).

Susceptibilities

(Subscript[χ, N])=(1/N)(Subscript[(((∂N)/(∂μ))), V,T]),
 (Subscript[χ, Q])=(1/N)(Subscript[(((∂Q)/(∂ℰ))), V,T])=(e^2)(Subscript[χ, N]),
 (Subscript[χ, M])=(1/N)(Subscript[(((∂M)/(∂B))), V,T])=(Subscript[μ, B]^2)(Subscript[χ, N]).

Special Cases

z->0: Dilute (classical) limit

Equation of state

P V = N T

Energy

E=Δ P V=Δ N T.

Specific heats

(Subscript[c, V])=Δ,
(Subscript[c, P])=Δ+1.

z->1, η=+1: Bose-Einstein condensation (BEC)

Condensed fractions

((Subscript[N, 0])/N)=1-((((T/(Subscript[T, c]))))^Δ).

Transition temperature

(Subscript[T, c])=({{{, ({{(Subscript[ϵ, 0])((((1/(Γ(Δ+1)ζ(Δ) ))(N/V)))^(1/Δ)), Δ>1}, {0, Δ<=1}})}}).

z=1, η=+1: Black-body radiation or in the BEC phase

Equation of state

P V=((ζ(Δ+1))/(ζ(Δ)))N T.

Energy

E=Δ κ V (T^(Δ+1)).

Specific heats

(Subscript[c, V])=(1+(1/Δ)) κ (T^Δ).

z→∞, η=-1: Degenerated fermi gas

Fermi energy

Definition (Subscript[ϵ, F])=μ(T=0),

Formula

(Subscript[ϵ, F])=(Subscript[ϵ, 0]) ((((N/V)))^(1/Δ)).

Equation of state (T≪(Subscript[ϵ, F]))

P V=(1/(Δ+1)) N (Subscript[ϵ, F]).

Energy (T≪(Subscript[ϵ, F]))

E=(Δ/(Δ+1)) N (Subscript[ϵ, F]).

Specific heats (T≪(Subscript[ϵ, F]))

(Subscript[c, V])=(((π^2) Δ  )/3)((T/(Subscript[ϵ, F]))).

Susceptibilities (T≪(Subscript[ϵ, F]))

(Subscript[χ, N])=(Δ/(Subscript[ϵ, F])),
(Subscript[χ, Q])=((Δ (e^2))/(Subscript[ϵ, F])),
  (Subscript[χ, M])=((Δ (Subscript[μ, B]^2))/(Subscript[ϵ, F])).