n=10;
Plot[n!/((n/2+ϵ/2)!(n/2-ϵ/2)!),{ϵ,-n,n},PlotLabel->Style["Multiplicity for N particle spin-1/2 paramagnet",Black,28],PlotRange->{{-n,n},{1,300}}(*,PlotLegends->Placed[{Style["",20], Style["E",20], Style["s2",20], Style["s3",20]},{0.8,0.8}]*), Frame->{{True,False},{True,False}},FrameLabel->{{Style["Ω(E,10)",Black,28],None},{Style["x*Emax/h",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black},{Red, Dashed},{Blue, Dashed},{Green, Dashed}}]
Export["C:\\Users\\Sina\\Google Drive\\Spring 2021 Classes\\Statistical Mechanics\\HW1_p2b_image.png",%,ImageResolution->500];

n=10;
Plot[2/Log[(n+ϵ)/(n-ϵ)],{ϵ,-n,n},PlotLabel->Style["T(E) for N=10 particle spin-1/2 paramagnet",Black,28],PlotRange->{{0,n},{0,50}}(*,PlotLegends->Placed[{Style["",20], Style["E",20], Style["s2",20], Style["s3",20]},{0.8,0.8}]*), Frame->{{True,False},{True,False}},FrameLabel->{{Style["(Subscript[k, B])T(E,10)",Black,28],None},{Style["E/h",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black},{Red, Dashed},{Blue, Dashed},{Green, Dashed}}]
Export["C:\\Users\\Sina\\Google Drive\\Spring 2021 Classes\\Statistical Mechanics\\HW1_p2c_image.png",%,ImageResolution->500];

n=100;
Plot[Evaluate[n Tanh[B/t]/.B->{0.1,0.5,1,5,10}],{t,0,10},PlotLabel->Style["m(T, B) for N=100 particle spin-1/2 paramagnet",Black,28],PlotRange->{{0,10},{0,125}},PlotLegends->Placed[{Style["(Subscript[μ, B])B=0.1",20], Style["(Subscript[μ, B])B=0.5",20], Style["(Subscript[μ, B])B=1.0",20], Style["(Subscript[μ, B])B=5.0",20], Style["(Subscript[μ, B])B=10",20]},{0.8,0.7}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["m/(Subscript[μ, B])",Black,28],None},{Style["(Subscript[k, B])T",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black},{Red, Dashed},{Blue, Dashed},{Green, Dashed},{Black,Dashed}}]
Export["C:\\Users\\Sina\\Google Drive\\Spring 2021 Classes\\Statistical Mechanics\\HW1_p2d1_image.png",%,ImageResolution->500];Plot[Evaluate[n Tanh[B/t]/.t->{0.1,0.5,1,5,10}],{B,0,10},PlotLabel->Style["m(T, B) for N=100 particle spin-1/2 paramagnet",Black,28],PlotRange->{{0,10},{0,125}},PlotLegends->Placed[{Style["(Subscript[k, B])T=0.1",20], Style["(Subscript[k, B])T=0.5",20], Style["(Subscript[k, B])T=1.0",20], Style["(Subscript[k, B])T=5.0",20], Style["(Subscript[k, B])T=10",20]},{0.8,0.3}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["m/(Subscript[μ, B])",Black,28],None},{Style["(Subscript[μ, B])B",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black},{Red, Dashed},{Blue, Dashed},{Green, Dashed},{Black,Dashed}}]
Export["C:\\Users\\Sina\\Google Drive\\Spring 2021 Classes\\Statistical Mechanics\\HW1_p2d2_image.png",%,ImageResolution->500];

Plot[Evaluate[{3 (ℏ ω0)/(Exp[ℏ ω0/x]-1),3x,3ℏ ω0 Exp[-ℏ ω0/x]}/.{ℏ->1, ω0->1}],{x,0,10},PlotLabel->Style["E(T) for N 3d harmonic oscillators",Black,28],PlotRange->{{0,3},{0,10}},PlotLegends->Placed[{Style["E(T)/N",20], Style["3(Subscript[k, B])T",20], Style["3(Subscript[ℏω, 0])(e^(-(Subscript[ℏω, 0])/(Subscript[k, B]) T))",20]},{0.2,0.75}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["E/N",Black,28],None},{Style["(Subscript[k, B])T",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black},{Red, Dashed},{Blue, Dashed},{Green, Dashed},{Black,Dashed}}]
Export["C:\\Users\\Sina\\Google Drive\\Spring 2021 Classes\\Statistical Mechanics\\HW1_p3c_image.png",%,ImageResolution->500];

Plot[Evaluate[{3(ℏ ω0)^2/(kB T^2)*Exp[ℏ ω0/(kB T)]/(Exp[ℏ ω0/(kB T)]-1)^2,3 kB,3(ℏ ω0)^2/(kB T^2)Exp[-ℏ ω0/(kB T)]}/.{ℏ->1, ω0->1,kB->1}],{T,0,10},PlotLabel->Style["(Subscript[C, V])(T) for N 3d harmonic oscillators",Black,28],PlotRange->{{0,2},{0,5}},PlotLegends->Placed[{Style["(Subscript[C, V])(T)/N",20], Style["3kB",20], Style["3((Subscript[ℏω, 0]))^2/((Subscript[k, B])(T^2))\!\(\*SuperscriptBox[\(e\), \(\(-ℏ\)\\ \*SubscriptBox[\(ω\), \(0\)]/\*SubscriptBox[\(k\), \(B\)] T\)]\)]",20]},{0.25,0.75}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["(Subscript[C, V])/N",Black,28],None},{Style["T",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black},{Red, Dashed},{Blue, Dashed},{Green, Dashed},{Black,Dashed}}]
Export["C:\\Users\\Sina\\Google Drive\\Spring 2021 Classes\\Statistical Mechanics\\HW1_p3d_image.png",%,ImageResolution->500];

S[T_,V_,n_]= 3 n kB/2(Log[kB T(1-2/(3n))]+Log[m/(2 Pi ℏ^2)])+n kB(Log[V]-Log[n])+kB Log[n/(kB T(3n/2-1))]+kB Log[3 Δ/2];
SE[En_,V_,n_]=3n kB/2 Log[2En/(3n)*(m/2 Pi ℏ^2)]+n kB Log[V/n]+kB Log[n/En*(3 Δ/2)]

T[En_]=D[SE[En, V, n],En]
Solve[T[En]==1/T,En]

Ω[En_]=Δ/En (m/(2 Pi ℏ^2))^(3 n/2)(3 n V^n/2)En^(3n/2)/(n! (3n/2)!)
ΩT[T_]=FullSimplify[Ω[kB T(3n-2)/2]]

soln=Solve[ΩT[T]==1,T][[1]][[1]][[2]]
FullSimplify[soln*(3n-2)kB]

Tc=N[2/(3n-2)(3 Δ n V^n/(n!(3n/2)!))^(2/(2-3n))(m/(2 Pi ℏ^2))^(-3n/(3n-2))/.{V->1,n->10,kB->1,m->1,ℏ->1, Δ->1}]

Plot[ΩT[T]/.{V->1,n->10,kB->1,m->1,ℏ->1, Δ->1},{T,0,11}]
FindRoot[ΩT[T]-1/.{V->1,n->10,kB->1,m->1,ℏ->1, Δ->1},{T,8.66}]

Plot[Evaluate[{S[T,V,n],2 Pi ℏ^2 n/(m n - 2m/3)(n/V)^(1-n/2)Δ^(3n/2-1)}/.{V->1,n->10,kB->1,m->1,ℏ->1, Δ->1}],{T,0,100}]

N[2 Pi ℏ^2 n/(m n - 2m/3)(n/V)^(1-n/2)Δ^(3n/2-1)/.{V->1,n->10,kB->1,m->1,ℏ->1, Δ->1}]

D[S[T,V,n],