In[]:=
d=1;​​T=0.0001d;​​ϵf0=-0.102001d;​​p=0.875;​​q0=0.5;​​v=0.5d;​​ds[ϵ_]:=
3
4Sqrt[2]d
Sqrt
(ϵ+d)
d
;​​​​(*wecalculaterenormalizedpositionofthef-level,chemicalpotentialandslavebosonamplitudebysolvingthreenonlinearequations*)​​​​ζ[ϵ_,μ_]:=ϵ-μ;(*ϵ=
2
k
2m
-d*)​​Rk[ϵ_,μ_,ϵf_,a_]:=Sqrt[
2
(ϵf-ζ[ϵ,μ])
+
2
(2va)
];​​Enkp[ϵ_,μ_,ϵf_,a_]:=
1
2
(ϵf+ζ[ϵ,μ]+Rk[ϵ,μ,ϵf,a]);​​Enkm[ϵ_,μ_,ϵf_,a_]:=
1
2
(ϵf+ζ[ϵ,μ]-Rk[ϵ,μ,ϵf,a]);
In[]:=
fwp[ϵ_,μ_,ϵf_,a_]:=
1
2
1-Tanh
Enkp[ϵ,μ,ϵf,a]
2T
;​​fwm[ϵ_,μ_,ϵf_,a_]:=
1
2
1-Tanh
Enkm[ϵ,μ,ϵf,a]
2T
;​​diff[ϵ_,μ_,ϵf_,a_]:=1-UnitStepϵ-
(a^2+4ϵfμ)
(4ϵf)
;​​sum[ϵ_,μ_,ϵf_,a_]:=1-UnitStepϵ-
(a^2+4ϵfμ)
(4ϵf)
;​​eq1int[ϵ_,μ_,ϵf_,a_]:=
diff[ϵ,μ,ϵf,a]
Rk[ϵ,μ,ϵf,a]
;​​eq2int[ϵ_,μ_,ϵf_,a_]:=(ζ[ϵ,μ]-ϵf)eq1int[ϵ,μ,ϵf,a];​​​​eq3int[ϵ_,μ_,ϵf_,a_]:=sum[ϵ,μ,ϵf,a];​​​​fint1[ϵ_,μ_,ϵf_,a_]:=ds[ϵ]eq1int[ϵ,μ,ϵf,a];​​fint2[ϵ_,μ_,ϵf_,a_]:=ds[ϵ]eq2int[ϵ,μ,ϵf,a];​​fint3[ϵ_,μ_,ϵf_,a_]:=ds[ϵ]eq3int[ϵ,μ,ϵf,a];
In[]:=
eqn1[μ_?NumericQ,ϵf_?NumericQ,a_?NumericQ]:=ϵf-ϵf0-
2
v
NIntegrate[fint1[ϵ,μ,ϵf,a],{ϵ,-d,d},MaxRecursion20,Method{GlobalAdaptive,MaxErrorIncreases8000}]​​eqn2[μ_?NumericQ,ϵf_?NumericQ,a_?NumericQ]:=2(q0-
2
a
)-p-NIntegrate[fint2[ϵ,μ,ϵf,a],{ϵ,-d,d},MaxRecursion20,Method{GlobalAdaptive,MaxErrorIncreases8000}]​​eqn3[μ_?NumericQ,ϵf_?NumericQ,a_?NumericQ]:=p-NIntegrate[fint3[ϵ,μ,ϵf,a],{ϵ,-d,d},MaxRecursion20,Method{GlobalAdaptive,MaxErrorIncreases8000}]
In[]:=
mysol=FindRoot[{eqn1[μ,ϵf,a],eqn2[μ,ϵf,a],eqn3[μ,ϵf,a]},{{μ,-0.42628},{ϵf,0.035940},{a,-0.37342}},WorkingPrecision16,AccuracyGoal15,MaxIterationsInfinity]​​
FindRoot
:Encountered a singular Jacobian at the point {μ,ϵf,a} = {0.845381,0.196724,-0.730453}. Try perturbing the initial point(s).
Out[]=
{μ0.845381,ϵf0.196724,a-0.730453}
In[]:=
{eqn1[μ,ϵf,a],eqn2[μ,ϵf,a],eqn3[μ,ϵf,a]}/.mysol
Out[]=
{0.0630848,-0.292974,-0.125}