Maths

ImF[α_,Δ_,H_][t_,h_]:=((Abs[t])^(2-α))H[h ((Abs[t])^(-Δ))]

ℋ[d_,A_,B_][X_]:=With[{ζ=If[d==3,-7/3,-(d-3)d/2]},-A HeavisideTheta[-X](((-B X))^ζ)Exp[-1/(((-B X))^(d-1))]]

F4d=Integrate[(1/π)((ImF[α,Δ,ℋ[4,A,B]][t,Xp ((Abs[t])^Δ)])/(Xp-X)),{Xp,-∞,∞},Assumptions->{B>0,X>0}]

FullSimplify[F4d/.X->Y/B,Assumptions->{Y>0}]/((Abs[t])^(2-α))/.Y->B X

F3d=Integrate[(1/π)((ImF[α,Δ,ℋ[3,A,B]][t,Xp ((Abs[t])^Δ)])/(Xp-X)),{Xp,-∞,∞},Assumptions->{B>0,X>0}]

FullSimplify[F3d/.X->Y/B,Assumptions->{Y>0}]/((Abs[t])^(2-α))/.Y->B X

χ2d=Integrate[(2/π)((ImF[α,Δ,ℋ[2,A,B]][t,Xp ((Abs[t])^Δ)])/(((Xp -X))^3)),{Xp,-∞,∞},Assumptions->{B>0,X>0}]

Limit[χ2d/((Abs[t])^(2-α)),X->0,Assumptions->{B>0}]

ℳ2d=C[0]((Abs[t])^(31/8-α))+Integrate[χ2d/.{X->h ((Abs[t])^(-Δ))},h]/.h->X ((Abs[t])^Δ)//Simplify

solc0=Solve[m0==Limit[ℳ2d/((Abs[t])^(31/8-α)),X->0,Assumptions->{B>0}],C[0]]

F2d=Integrate[ℳ2d/.solc0[[1]]/.X->h ((Abs[t])^(-Δ)),h]/.h->X ((Abs[t])^Δ)//Simplify

Limit[F2d/((Abs[t])^(23/4-α)),X->0,Assumptions->{B>0}]

Ssurf[d_]:=((2 (π^(d/2)))/(Gamma[d/2]))

Svol[d_]:=((π^(d/2))/(Gamma[d/2+1]))

ΔF=Ssurf[d] (R^(d-1)) Σ-Svol[d] (R^d)M H;

solR=Solve[D[ΔF,R]==0,R]//FullSimplify

ΔFcrit=FullSimplify[ΔF/.solR[[1]]]

b1=T FullSimplify[((FullSimplify[ΔFcrit (X^(d-1))/.{M->(t^β)ℳ[0],Σ->(t^μ)𝒮[0],H->h T}/.μ->-ν+γ+2β/.γ->2-α-2β/.ν->(2-α)/d/.α->2-β(1+δ)/.h->X (t^(β δ)),Assumptions->{t>0,d>1,𝒮[0]>0,ℳ[0]>0,T>0,h>0,β>0,δ>0,X>0}])^(-1/(d-1)))/.((x_)^(-1+d)):>(((1/x))^(1-d)),Assumptions->{t>0,d>1,𝒮[0]>0,ℳ[0]>0,T>0,h>0,β>0,δ>0,X>0}]

b2=FullSimplify[b1/.{𝒮[0]->4,ℳ[0]->ℳ0,d->2,T->Tc}]

Fitting

Needs["ErrorBarPlots`"];

Tc:=2/Log[1+(Sqrt[2])]

t[Tc_][T_]:=Abs[((T-Tc)/T)]

poly[{As__},λ_][X_]:=Total@({As}[[#]](X^(#))&/@Range[Length[{As}]])/(1+(((λ X))^(Length[{As}]+1)))

polyint[{As__},λ_][X_]:=With[{n=Length[{As}]},Total@({As}[[#]](X^(1+#)) Hypergeometric2F1[1,((1+#)/(1+n)),((2+#+n)/(1+n)),-(((X λ))^(1+n))]/(1+#)&/@Range[n])]

poly[{As___},λ_][X_]:=Total@({As}[[#]](X^(#))/(1+(((λ X))^(#+1)))&/@Range[Length[{As}]])

polyint[{As___},λ_][X_]:=Total@({As}[[#]]((λ^(-(#+1)))/(#+1))Log[1+(((λ X))^(#+1))]&/@Range[Length[{As}]])

f[A_,B_,{As___},λ_][X_]:=((A (B^2))/(π (((B X))^3)))(B X(B X-1)-Exp[1/(B X)]ExpIntegralEi[-1/(B X)])+poly[{As},λ][B X]

ℳ[ℳ0_,A_,B_,{As___},λ_][X_]:=((A B)/π)(1-(1/(B X))(B X-1)Exp[1/(B X)]ExpIntegralEi[-1/(B X)])+ℳ0+polyint[{As},λ][B X]/B

ℳ0:=(((2((Sqrt[2])-1)))^(1/4))((((4+3(Sqrt[2]))ArcSinh[1]))^(1/8))

B:=(Tc^2)ℳ0/((4^2) π)

χ0:=0.0255369719/Tc

Δ:=15/8

γ:=7/4

β:=1/8

data=Import["~/doc/research/sethna/doc/essential-ising_doc/data/data_square-4096.dat"];

assData=Thread[Rule[{"L","T","H","E","δE","M","δM","C","δC","χ","δχ"},#]]&/@data;

htΔvχtγ={{"H"/"T"(((t[Tc]["T"]))^(-Δ)),"χ" (((t[Tc]["T"]))^γ)},ErrorBar["δχ"(((t[Tc]["T"]))^γ)]}/.#&/@assData;

loghtΔvlogχtγ={Log10@{"H"/"T"(((t[Tc]["T"]))^(-Δ)),"χ" (((t[Tc]["T"]))^γ)},ErrorBar["δχ"/"χ"]}/.#&/@assData;

htΔvMtβ={{"H"/"T"(((t[Tc]["T"]))^(-Δ)),"M" (((t[Tc]["T"]))^(-β))},ErrorBar["δM"(((t[Tc]["T"]))^(-β))]}/.#&/@assData;

loghtΔvMtβ={{Log10["H"/"T"(((t[Tc]["T"]))^(-Δ))],"M" (((t[Tc]["T"]))^(-β))},ErrorBar["δM"(((t[Tc]["T"]))^(-β))]}/.#&/@assData;

loghtΔvlogMtβ={Log10@{"H"/"T"(((t[Tc]["T"]))^(-Δ)),"M" (((t[Tc]["T"]))^(-β))},ErrorBar["δM"/"M"]}/.#&/@assData;

χfitData={"H"/"T" ((t[Tc]["T"])^(-Δ)),"χ" ((t[Tc]["T"])^γ)}/.#&/@assData;

χfitWeights=1/((("δχ"((t[Tc]["T"])^γ)))^2)/.#&/@assData;

MfitData={"H"/"T" ((t[Tc]["T"])^(-Δ)),"M" ((t[Tc]["T"])^(-β))}/.#&/@assData;

MfitWeights=1/((("δM"((t[Tc]["T"])^(-β))))^2)/.#&/@assData;

fitData=Join[Prepend[#,1]&/@χfitData,Prepend[#,2]&/@MfitData];

fitWeights=Join[χfitWeights,MfitWeights];

fitFunc[i_][A_,B_,ℳ0_,{As___},λ_][X_]:=KroneckerDelta[i,1]f[A,B,{As},λ][X]+KroneckerDelta[i,2]ℳ[ℳ0,Tc A,B,Tc{As},λ][X]

fitFunc[i_][A_,B_,ℳ0_,c0_,{As___},λ_][X_]:=KroneckerDelta[i,1]((1-c0/(2 A (B^2)/π))f[A,B,{As},λ][X]+c0/(1+λ B X))+KroneckerDelta[i,2]((1-c0/(2 A (B^2)/π))ℳ[ℳ0,Tc A,B,Tc{As},λ][X]+c0/(2 A (B^2)/π)ℳ0+Tc c0 Log[1+B X λ]/(B λ))

nn:=3

fit=NonlinearModelFit[Prepend[#,1]&/@χfitData,fitFunc[i][π χ0/(2 (B^2)),B,ℳ0,A/@Range[nn],λ][x],Join[{λ},Thread[{A/@Range[nn],ConstantArray[0,nn]}]],{i,x},PrecisionGoal->70,MaxIterations->1000,Method->Automatic,Weights->χfitWeights];

fit["ParameterTable"]

pχlog=Show[ErrorListPlot[loghtΔvlogχtγ,PlotRange->All],Plot[{Log10@f[π χ0/(2(B^2)),B,{},0][(10^X)],Log10@fit[1,10^X]},{X,-3,4},PlotStyle->{Black,{Black,Dashed}}]]

pχ=Show[ErrorListPlot[htΔvχtγ,PlotRange->All],Plot[{f[π χ0/(2(B^2)),B,{},0][X],fit[1,X]},{X,0,2000},PlotPoints->500,PlotRange->{0,0.02},PlotStyle->{Black,{Black,Dashed}}]]

pMlog=Show[ErrorListPlot[loghtΔvMtβ,PlotRange->All],Plot[{ℳ[ℳ0,Tc π χ0/(2(B^2)),B,{},0][(10^X)],fit[2,10^X]},{X,-6,40},PlotRange->All,PlotStyle->{Black,{Black,Dashed}}]]

pM=Show[ErrorListPlot[{htΔvMtβ},PlotRange->All],Plot[{ℳ[ℳ0,Tc π χ0/(2(B^2)),B,{},0][X],fit[2,X]},{X,0,1200},PlotRange->All,PlotStyle->{Black,{Black,Dashed}}]]