En[kBT_]=-J Tanh[J/kBT]
Cv[kBT_] = (J/kBT)^2/Cosh[J/kBT]^(2)

Numerator[(((J^2) ((Sech[(J/kBT)])^2))/(kBT^2))]

Plot[En[T]/.J->1,{T,0,10},PlotLabel->Style["Total Energy and Heat Capacity (per spin)",Black,28],PlotRange->{{0,3},{-1.5,0.5}},PlotLegends->Placed[{},{0.9,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["(Subscript[k, B])T",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Darker[Green]},{Blue}},Epilog->{Gray,Dashed,Line[{{1,-10},{1,5}}]}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p1(1)_image.png",%,ImageResolution->500];

Plot[Cv[T]/.J->1,{T,0,10},PlotLabel->Style["Heat Capacity per Spin",Black,28],PlotRange->{{0,3},{0,0.5}},PlotLegends->Placed[{},{0.9,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["(Subscript[C, V])/(Subscript[Nk, B])",Black,28],None},{Style["(Subscript[k, B])T",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{Blue},Epilog->{Gray,Dashed,Line[{{1,-10},{1,10}}]}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p1(2)_image.png",%,ImageResolution->500];

UnitConvert[Quantity[1,"ElementaryCharge"^2/"BohrRadius"/"ElectricConstant"]/(4 Pi),"ElectronVolts"]

UnitConvert[Quantity[1,"MagneticConstant""BohrMagneton"^2/"BohrRadius"^3]/(4 Pi),"ElectronVolts"]

N[UnitConvert[Quantity[300,"BoltzmannConstant"*"Kelvins"],"ElectronVolts"]]

γ[α_,κ_,B_]= FullSimplify[Coth[α Beff]-Coth[α Beff/2]/2/.{Beff->(κ m+B),J0->κ μ^2}]

γ[α,κ,B]

Plot[Evaluate[{m,γ[0.25,1,0],γ[0.5,1,0],γ[1,1,0],γ[2,1,0],γ[5,1,0],γ[10,1,0]}],{m,-5,5},PlotLabel->Style["MFT Magnetization at B=0",Black,28],PlotRange->{{-3,3},{-0.75,0.75}},PlotLegends->Placed[{Style["m",20],Style["(1/2)tanh[((ακm)/2)], ακ=(1/4)",20],Style["(1/2)tanh[((ακm)/2)], ακ=(1/2)",20],Style["(1/2)tanh[((ακm)/2)], ακ=1",20],Style["(1/2)tanh[((ακm)/2)], ακ=2",20],Style["(1/2)tanh[((ακm)/2)], ακ=5",20],Style["(1/2)tanh[((ακm)/2)], ακ=10",20]},{1.1,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {1100,400},AspectRatio->Full,PlotStyle->{{Black,Dashed},{Pink},{Green},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2c(1)_image.png",%,ImageResolution->500];

Series[Tanh[x],{x,0,4}]

B=1;
Plot[Evaluate[{m,γ[0.25,1,B],γ[0.5,1,B],γ[1,1,B],γ[2,1,B],γ[5,1,B],γ[10,1,B]}],{m,-5,5},PlotLabel->Style["MFT Magnetization at B=1 (κ=1)",Black,28],PlotRange->{{-3,3},{-0.75,0.75}},PlotLegends->Placed[{Style["m",20],Style["(1/2)tanh[((αB + ακm)/2)], α=(1/4)",20],Style["(1/2)tanh[((αB + ακm)/2)], α=(1/2)",20],Style["(1/2)tanh[((αB + ακm)/2)], α=1",20],Style["(1/2)tanh[((αB + ακm)/2)], α=2",20],Style["(1/2)tanh[((αB + ακm)/2)], α=5",20],Style["(1/2)tanh[((αB + ακm)/2)], α=10",20]},{1.1,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {1100,400},AspectRatio->Full,PlotStyle->{{Black,Dashed},{Pink},{Green},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2c(2)_image.png",%,ImageResolution->500];

B=-1;
Plot[Evaluate[{m,γ[0.25,1,B],γ[0.5,1,B],γ[1,1,B],γ[2,1,B],γ[5,1,B],γ[10,1,B]}],{m,-5,5},PlotLabel->Style["MFT Magnetization at B = –1 (κ=1)",Black,28],PlotRange->{{-3,3},{-0.75,0.75}},PlotLegends->Placed[{Style["m",20],Style["(1/2)tanh[((αB + ακm)/2)], α=(1/4)",20],Style["(1/2)tanh[((αB + ακm)/2)], α=(1/2)",20],Style["(1/2)tanh[((αB + ακm)/2)], α=1",20],Style["(1/2)tanh[((αB + ακm)/2)], α=2",20],Style["(1/2)tanh[((αB + ακm)/2)], α=5",20],Style["(1/2)tanh[((αB + ακm)/2)], α=10",20]},{1.1,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {1100,400},AspectRatio->Full,PlotStyle->{{Black,Dashed},{Pink},{Green},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2c(3)_image.png",%,ImageResolution->500];

Plot[Evaluate[{m,γ[0.5,1,-0.2],γ[0.5,1,-0.1],γ[0.5,1,0],γ[0.5,1,0.1],γ[0.5,1,0.2]}],{m,-5,5},PlotLabel->Style["MFT Magnetization at α = 0.5 (κ=1)",Black,28],PlotRange->{{-0.5,0.5},{-0.2,0.2}},PlotLegends->Placed[{Style["",20],Style["B=-0.2",20],Style["B=-0.1",20],Style["B=+±0.0",20],Style["B=+0.1",20],Style["B=+0.2",20]},{1.1,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black,Dashed},{Pink},{Green},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2c(4)_image.png",%,ImageResolution->500];Plot[Evaluate[{m,γ[5,1,-0.2],γ[5,1,-0.1],γ[5,1,0],γ[5,1,0.1],γ[5,1,0.2]}],{m,-5,5},PlotLabel->Style["MFT Magnetization at α = 5 (κ=1)",Black,28],PlotRange->{{-1,1},{-0.75,0.75}},PlotLegends->Placed[{Style["",20],Style["B=-0.2",20],Style["B=-0.1",20],Style["B=+±0.0",20],Style["B=+0.1",20],Style["B=+0.2",20]},{1.1,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {800,400},AspectRatio->Full,PlotStyle->{{Black,Dashed},{Pink},{Green},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2c(5)_image.png",%,ImageResolution->500];

Plot3D[Evaluate[{m,γ[0.25,1,b],γ[0.5,1,b],γ[1,1,b],γ[2,1,b],γ[5,1,b],γ[10,1,b]}],{m,-5,5},{b,-3,3},PlotLabel->Style["MFT Magnetization(κ=1)",Black,28],PlotRange->{{-3,3},{-0.75,0.75}},PlotLegends->Placed[{Style["m",20],Style["(1/2)tanh[((αB + ακm)/2)], α=(1/4)",20],Style["(1/2)tanh[((αB + ακm)/2)], α=(1/2)",20],Style["(1/2)tanh[((αB + ακm)/2)], α=1",20],Style["(1/2)tanh[((αB + ακm)/2)], α=2",20],Style["(1/2)tanh[((αB + ακm)/2)], α=5",20],Style["(1/2)tanh[((αB + ακm)/2)], α=10",20]},{1.1,0.5}],ImageSize-> {1100,400},AspectRatio->Full,PlotStyle->Automatic]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2c(4)_image.png",%,ImageResolution->500];

Plot3D[Evaluate[{m,γ[0.25,1,g],γ[0.5,1,g],γ[1,1,g],γ[2,1,g],γ[5,1,g],γ[10,1,g]}],{m,-5,5},{g,-3,3},PlotRange->{-1,1}]

FullSimplify[Sinh[μ Beff β]/Sinh[μ Beff β/2]]

D[Cosh[μ β f[x]/2],x]

f[m_,B_,κ_,kBT_]=κ m^2/2-1/kBT Log[2 Cosh[(κ m + B)/(2 kBT)]];
Plot3D[{f[m,B,1,0.5],f[m,B,1,1],f[m,B,1,5]},{m,-3,3},{B,-1,1},PlotLegends->"Expressions",AxesLabel->Automatic]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2d(1)_image.png",%,ImageResolution->500];
Plot3D[{f[m,B,1,0.5]},{m,-3,3},{B,-1,1},PlotLegends->"Expressions",AxesLabel->Automatic]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2d(2)_image.png",%,ImageResolution->500];
Plot3D[{f[m,B,1,1]},{m,-3,3},{B,-1,1},PlotLegends->"Expressions",AxesLabel->Automatic]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2d(3)_image.png",%,ImageResolution->500];
Plot3D[{f[m,B,1,5]},{m,-3,3},{B,-1,1},PlotLegends->"Expressions",AxesLabel->Automatic]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2d(4)_image.png",%,ImageResolution->500];

Plot[Evaluate[{f[m,-0.5,1,0.5],f[m,0,1,0.5],f[m,0.5,1,0.5]}],{m,-20,20},PlotLabel->Style["Free Energy at (Subscript[k, B])T=0.5 (κ=1)",Black,28],PlotRange->{{-6,6},{-4,3}},PlotLegends->Placed[{Style["B = -0.5",20],Style["B =  0.0",20],Style["B = +0.5",20],Style["B=+±0.0",20],Style["B=+0.1",20],Style["B=+0.2",20]},{0.65,0.75}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {600,400},AspectRatio->Full,PlotStyle->{{Blue},{Black,Dashed},{Red},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2d(5)_image.png",%,ImageResolution->500];Plot[Evaluate[{f[m,-0.5,1,1],f[m,0,1,1],f[m,0.5,1,1]}],{m,-20,20},PlotLabel->Style["Free Energy at (Subscript[k, B])T=1 (κ=1)",Black,28],PlotRange->{{-6,6},{-1,3}},PlotLegends->Placed[{Style["B = -0.5",20],Style["B =  0.0",20],Style["B = +0.5",20],Style["B=+±0.0",20],Style["B=+0.1",20],Style["B=+0.2",20]},{0.85,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {600,400},AspectRatio->Full,PlotStyle->{{Blue},{Black,Dashed},{Red},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2d(6)_image.png",%,ImageResolution->500];Plot[Evaluate[{f[m,-0.5,1,5],f[m,0,1,5],f[m,0.5,1,5]}],{m,-20,20},PlotLabel->Style["Free Energy at (Subscript[k, B])T=5 (κ=1)",Black,28],PlotRange->{{-6,6},{-1,3}},PlotLegends->Placed[{Style["B = -0.5",20],Style["B =  0.0",20],Style["B = +0.5",20],Style["B=+±0.0",20],Style["B=+0.1",20],Style["B=+0.2",20]},{0.85,0.5}], Frame->{{True,False},{True,False}},FrameLabel->{{Style["",Black,28],None},{Style["",Black,28], None}}, FrameTicks->All,FrameTicksStyle->Directive[Gray,14],ImageSize-> {600,400},AspectRatio->Full,PlotStyle->{{Blue},{Black,Dashed},{Red},{Cyan},{Blue},{Purple},{Black}}]
Export["C:\\Users\\Sina\\Google Drive\\Classes\\Spring 2021\\Statistical Mechanics\\HW6_p2d(7)_image.png",%,ImageResolution->500];

Series[Tanh[y(1+x)],{x,0,4}]

Series[Log[2 Cosh[α/2 (B +κ m)]],{m,0,4}]

FullSimplify[m-μ/2*Series[Tanh[μ β/2 (κ m +B)],{m,0,3}]/.{κ->J0/μ^2}]

Collect[FullSimplify[m-μ/2(x-x^3/3)/.x->μ β(κ m +B)/2],m]/.κ->(J0/μ^2)

Series[Log[2Cosh[x]],{x,0,4}]

Collect[FullSimplify[κ/2 m^2-1/β(Log[2]+x^2/2-x^4/12)/.x->μ β(B+κ m)/2]/.κ->J0/μ^2,m]

4^3/48

Series[Sqrt[1-x],{x,0,3}]

48/4^3

f[m_,B_,κ_,kBT_]=κ m^2/2-1/kBT Log[2 Cosh[(κ m + B)/(2 kBT)]];
FullSimplify[Series[f[m,0,κ,kBT],{m,0,4}]/κ]/.κ->J0/μ^2
f2=m^2(1-J0 β)/2+b/4 m^4
D[f2,m]
FullSimplify[Solve[%==0,m]/.{J0->4 kB Tc,β->1/(kB T)}]

Solve[(Collect[m-μ/2(μ β (J0 m/μ^2)/2-(μ β (J0 m/μ^2)/2)^3/3),m]/.{J0->(4 kB Tc),β->1/(kB T)})==0,m]

FullSimplify[Solve[x(-μ/2+m^2/(2μ))+4/3 m^3/μ^2==0,m]/.x->μ B β/2]

Series[Cosh[x],{x,0,4}]

Collect[Expand[(x^2+δx^2-2 x δx Cos[θ])^2],δx]

f=tx mx^2/2+ty my^2/2 + ux mx^4/4+ uy my^4/4 +uxy mx^2 my^2/2
FullSimplify[Solve[{D[f,mx]==0,D[f,my]==0},{mx,my},Reals],Assumptions->{tx<0,ty>0,ux>0,uy>0,uxy>0,uxy^2>ux uy}]
FullSimplify[f/.%]

tx=5;ty=10;
ux=1;uy=2;
uxy=10;
Plot3D[f,{mx,-3,3},{my,-3,3},AxesLabel->Automatic]

tx=-5;ty=10;
ux=1;uy=2;
uxy=10;
Plot3D[f,{mx,-4,4},{my,-3,3},AxesLabel->Automatic,PlotRange->{-20,100}]

tx=5;ty=-10;
ux=1;uy=2;
uxy=10;
Plot3D[f,{mx,-3,3},{my,-3,3},AxesLabel->Automatic]

tx=-5;ty=-10;
ux=1;uy=2;
uxy=10;
Plot3D[f,{mx,-3,3},{my,-3,3},AxesLabel->Automatic,PlotRange->{-20,75}]

{f/.{mx-> -Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],my-> -Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]},-Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],-Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]}

Graphics3D[{{Red,PointSize[0.1],Point[{f/.{mx-> Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],my-> Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]},Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]}]},{Red,PointSize[0.1],Point[{f/.{mx-> -Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],my-> Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]},-Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]}]},{Red,PointSize[0.1],Point[{f/.{mx-> Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],my-> -Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]},Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],-Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]}]},{Red,PointSize[0.1],Point[{f/.{mx-> -Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],my-> -Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]},-Sqrt[(uxy ty-uy tx)/(ux uy - uxy^2)],-Sqrt[(uxy tx-ux ty)/(ux uy - uxy^2)]}]}}]