Møller-Plesset and density-fixed adiabatic connections for a model diatomic system at different correlation regimes

Summary

Authors:
Sara Giarrusso and Aurora Pribram-Jones

This Mathematica notebook provides the analytical expressions and numerical calculations in support of the accompanying manuscript: 
“Møller-Plesset and density-fixed adiabatic connections for a model diatomic system at different correlation regimes”
Readers are invited to reproduce, and extend, the results contained in this notebook. 
Each cell must be evaluated sequentially to ensure the correct functionality.

Supporting notebook

MaTeX

Needs["MaTeX`"]
SetOptions[MaTeX,"Preamble"->{"\\usepackage{amssymb,amsmath,latexsym,amsfonts,amsthm,xcolor,bm,mhchem,physics}"}];

Interacting Hamiltonian

H3D[U_,t_,v0_,v1_]={{U + 2 v0,0,-(Sqrt[2 ])t},{0,U + 2 v1,-(Sqrt[2 ])t},{-(Sqrt[2 ])t,-(Sqrt[2 ])t,v0+v1}};

We define (v1-v0)=Δv=V

Ground State quantities

GS energy

ϵ[U_,V_]=(2/3) (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]);

The above energy expression assumes that the gauge of the potential is v1+v0=0

GS solution

normpsi0squared[U_, V_] ={((-2 ϵ+(2/(-2 U-2 V+2 ϵ)))/(Sqrt[2])),((Sqrt[2])/(2 U+2 V-2 ϵ)),1}.{((-2 ϵ+(2/(-2 U-2 V+2 ϵ)))/(Sqrt[2])),((Sqrt[2])/(2 U+2 V-2 ϵ)),1}/. ϵ-> ϵ[U,V];

Psi0[U_, V_] =(1/(Sqrt[normpsi0squared[U,V]])){((-2 ϵ+(2/(-2 U-2 V+2 ϵ)))/(Sqrt[2])),((Sqrt[2])/(2 U+2 V-2 ϵ)),1}/. ϵ-> ϵ[U,V];

Interacting site occupation difference

nn[U_,V_]=(4 V (-U-2 (Sqrt[3+3 (V^2)+(U^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9-18 (V^2)+2 (U^2)))/(2 (((3+3 (V^2)+(U^2)))^(3/2))))])]))/(3 ((V^2)+(1/9) (((U+2 (Sqrt[3+3 (V^2)+(U^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9-18 (V^2)+2 (U^2)))/(2 (((3+3 (V^2)+(U^2)))^(3/2))))])]))^2) (1+(9/(4 (((U-(Sqrt[3+3 (V^2)+(U^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9-18 (V^2)+2 (U^2)))/(2 (((3+3 (V^2)+(U^2)))^(3/2))))])]))^2))))));

Møller-Plesset adiabatic connection

RHF quantities

TBV[x_]=- (((*2t*) x)/(Sqrt[4 -(x^2)])) (*Tigh-binding potential as a function of the occupation x=Δn*);

VvsNhf[U_,x_]=TBV[x]-U*(x/2)(*"VvsNhf" stands for external potential, V, given the HF density, "Nhf"=x=((Δn)^HF)*);

NhfvsVn[ U_,x_]= (*this is the negative branch of the inverse of VvsNhf*)(1/(3 U))(-3 x-(Sqrt[3]) √(3 (x^2)+2 (-1+(U^2)-(x^2))+(((1-(U^2)+(x^2)))^2)/(((-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-16 (x^2)+(x^4))-6 (Sqrt[3]) U x (Sqrt[-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-7 (x^2)+(x^4))])))^(1/3))+(((-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-16 (x^2)+(x^4))-6 (Sqrt[3]) U x (Sqrt[-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-7 (x^2)+(x^4))])))^(1/3)))+(Sqrt[3]) √(-4+4 (U^2)+2 (x^2)-(((1-(U^2)+(x^2)))^2)/(((-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-16 (x^2)+(x^4))-6 (Sqrt[3]) U x (Sqrt[-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-7 (x^2)+(x^4))])))^(1/3))-(((-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-16 (x^2)+(x^4))-6 (Sqrt[3]) U x (Sqrt[-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-7 (x^2)+(x^4))])))^(1/3))-(6 (Sqrt[3]) (1+(U^2)) x)/(√(3 (x^2)+2 (-1+(U^2)-(x^2))+(((1-(U^2)+(x^2)))^2)/(((-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-16 (x^2)+(x^4))-6 (Sqrt[3]) U x (Sqrt[-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-7 (x^2)+(x^4))])))^(1/3))+(((-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-16 (x^2)+(x^4))-6 (Sqrt[3]) U x (Sqrt[-(U^6)+3 (U^4) (1+(x^2))+(((1+(x^2)))^3)-3 (U^2) (1-7 (x^2)+(x^4))])))^(1/3))))));

NhfvsV[U_, V_]=Piecewise[{{NhfvsVn[U,V], V≤0}},-NhfvsVn[U,-V]] (*full inverse*)

Vλhfgen[U_,V_,λ_]=TBV[NhfvsV[U, V]]- λ U ((NhfvsV[U, V])/2);(*MP AC λ-dependent external potential for the Hamiltonian with interaction λ U*)

Møller-Plesset fluctuation potential and (Subscript[c, 0]^HF)

n03D=DiagonalMatrix[{2,0,1}];
n13D=DiagonalMatrix[{0,2,1}];

Umat3d[U_]=(({{U, 0, 0}, {0, U, 0}, {0, 0, 0}}));

MPFP3D[U_,x_] = Umat3d[U] - (U/2)((1-(x/2))n03D + (1+(x/2))n13D);

c0hf[U_,x_]=(U/2)(1+((x^2)/4));

AC integrand as a function of external potential

Wλhf[U_,V_,λ_]= Psi0[λ U,Vλhfgen[U,V,λ]].MPFP3D[U, NhfvsV[U,V]].Psi0[λ U,Vλhfgen[U,V,λ]]+c0hf[U, NhfvsV[U,V]];

AC integrand as a function of HF density ("Nhf"=x=((Δn)^HF))

WλhfvsNhf[U_,x_,λ_]=-(U/3)+(1/8) U (4+(x^2))+(24 U (Sqrt[3+(U^2) (λ^2)+3 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)]) (-64 (Sqrt[4-(x^2)])+(U^3) (x^8) (λ^3)+2 U (x^6) λ (3+U (Sqrt[4-(x^2)]) λ-6 (U^2) (λ^2))+16 U (x^4) λ (-6-2 U (Sqrt[4-(x^2)]) λ+3 (U^2) (λ^2))+16 (x^2) (3 (Sqrt[4-(x^2)])+18 U λ+6 (U^2) (Sqrt[4-(x^2)]) (λ^2)-4 (U^3) (λ^3))) Cos[(1/6) (π+2 ArcCos[-((U λ (9+2 (U^2) (λ^2)-18 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)))/(2 (((3+(U^2) (λ^2)+3 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)))^(3/2))))])])/((((48+12 U (x^2) (Sqrt[4-(x^2)]) λ+(U^2) (16+8 (x^2)-3 (x^4)) (λ^2)))^(5/2)) (Sqrt[1-(((U^2) (λ^2) (((9+2 (U^2) (λ^2)-18 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)))^2))/(4 (((3+(U^2) (λ^2)+3 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)))^3)))]))-((U (2 U λ+3 (x^2) ((1/(Sqrt[4-(x^2)]))+((U λ)/2))) Sin[(1/6) (π+2 ArcCos[-((U λ (9+2 (U^2) (λ^2)-18 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)))/(2 (((3+(U^2) (λ^2)+3 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)))^(3/2))))])])/(3 (Sqrt[3+(U^2) (λ^2)+3 (x^2) ((((1/(Sqrt[4-(x^2)]))+((U λ)/2)))^2)])));

Density-fixed adiabatic connection

Lieb Maximisation

VλLM[U_, Nn_, λ_]:=Return[V/.NMaximize[ϵ[λ U,V]-(V/2)Nn,V][[2]]]

AC integrand as a function of external potential

PsiλKS[U_,V_,λ_]:=Return[Psi0[λ U, VλLM[U, nn[U, V],λ]]]

c0KSvsV[U_,V_]= - (U/2)(1+(((((nn[U,V])/2)))^2));

WλKS[U_,V_,λ_]:= PsiλKS[U, V, λ].Umat3d[U].PsiλKS[U, V, λ]+c0KSvsV[U,V]

AC integrand as a function of interacting density ("N"=x=((Δn)^HF))

PsiλKSvsN[U_,x_,λ_]:=Return[Psi0[λ U, VλLM[U,x,λ]]]

c0KSvsN[U_,x_]=- (U/2)(1+((((x/2)))^2));

WλKSvsN[U_,x_,λ_]:=PsiλKSvsN[U, x, λ].Umat3d[U].PsiλKSvsN[U, x, λ]+c0KSvsN[U,x]

Correlation energies, their separate contributions and their weak- and strong-interaction expansions [see also J. Chem. Phys. 157, 054102 (2022)]

Note that, for each of the 3 following ingredients,  (Subscript[I, c]^HF)(λ U, x)=(Subscript[I, c,λ]^HF)(U,x)  with I=U,T,V

(Subscript[U, c]^HF)[eq (53a)]

UchfvsNhf[ U_,x_]=-(1/8) U (4+(x^2))+(U (9/(((U+(3/2) x (U+(2/(Sqrt[4-(x^2)])))+2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]))^2)+(((((4 U)/3)-(4/3) (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]+3/(U+((3 U x)/2)+((3 x)/(Sqrt[4-(x^2)]))+2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])])))^2)))/(2+9/(((U+(3/2) x (U+(2/(Sqrt[4-(x^2)])))+2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]))^2)+(((((4 U)/3)-(4/3) (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]+3/(U+((3 U x)/2)+((3 x)/(Sqrt[4-(x^2)]))+2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])])))^2));

small U [compare eqs (51) and (54)]

Assuming[U>0&& -2<x<2, FullSimplify[Series[UchfvsNhf[U,x],{U,0,3}]]]

large U [compare eq (80) and preceding line]

The expectation value of the interaction operator, (OverHat[U]), on the asymptotic wavefunction, {0,0,1}, is clearly zero, then what remains from the 2-body term of the AC integrand is -<(Φ^HF)| (OverHat[U])|(Φ^HF)> =- (Subscript[c, 0]^HF)- (Subscript[U, H]) [eq (68)]

Assuming[U>0&&0<x<2, Simplify[Psi0[0,TBV[x]].Umat3d[U].Psi0[0,TBV[x]]]]

(Subscript[T, c]^HF)[eq (53b)]

TchfvsN[ U_,x_]=(Sqrt[1-((x^2)/4)])+(4 (2 U-2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]))/(3 (2+9/(((U+(3/2) x (U+(2/(Sqrt[4-(x^2)])))+2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]))^2)+(((((4 U)/3)-(4/3) (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]+3/(U+((3 U x)/2)+((3 x)/(Sqrt[4-(x^2)]))+2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])])))^2)));

small U [eq (56a)]

Assuming[U>0&& -2<x<2, FullSimplify[Series[TchfvsN[U,x],{U,0,3}]]]

large U

Assuming[ -2<x<2,Limit[TchfvsN[U,x],U -> ∞]]

The leading-order coefficient does not depend on interaction parameter (therefore it is subleading with respect to the large-U expansions of (Subscript[U, c]^HF) and (Subscript[V, c]^HF))

(Subscript[V, c]^HF)[eq (53c)]

VchfvsNhf[U_,x_]=(1/2) (-((U x)/2)-(x/(Sqrt[4-(x^2)]))) (-x+(4 (-((U x)/2)-(x/(Sqrt[4-(x^2)]))) (-U-2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]))/(3 ((((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)+(1/9) (((U+2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]))^2) (1+9/(((2 U-2 (Sqrt[3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (((((U x)/2)+(x/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]))^2)))));

small U [eq (56b)]

Assuming[U>0&& -2<x<2, FullSimplify[Series[VchfvsNhf[U,x],{U,0,3}]]]

large U [compare eq (80) and preceding line]

Using the asymptotic wavefunction {0, 0,1} in the one-body term of eq (39), one obtains

Assuming[0<x<2&&U>0,FullSimplify[{0, 0,1}. H3D[0 ,0,-(U/2) (1-(x/2)),-(U/2) (1+(x/2))].{0, 0,1}]]

(Subscript[E, c]^SD)

(Subscript[E, c]^HF) as a function of HF site-occupation difference

EchfvsNhf[U_,x_]=(1/24) (4 U+3 U (x^2)+(48/(Sqrt[4-(x^2)]))-16 (Sqrt[3+(U^2)+3 (x^2) ((((U/2)+(1/(Sqrt[4-(x^2)]))))^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (x^2) ((((U/2)+(1/(Sqrt[4-(x^2)]))))^2)))/(2 (((3+(U^2)+3 (x^2) ((((U/2)+(1/(Sqrt[4-(x^2)]))))^2)))^(3/2))))])]);

(Subscript[E, c]^HF) as a function of external potential

Echf[U_,V_]= (*Tchf+Uchf+Vchf*)(U (1+(((1+(4/9) (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]])))^2)))/(1+(2/9) (((U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]))^2)+(((1+(4/9) (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]])))^2))+(8 (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (((U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]))^2))/(3 (9+2 (((U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]))^2)+9 (((1+(4/9) (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]])))^2)))+√(1-(1/4) (((({{\[Piecewise], ({{(1/(3 U))(-3 V+(Sqrt[3]) √(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(Sqrt[3]) √(6 (V^2)+4 (-1+(U^2)-(V^2))-((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3))+(6 (Sqrt[3]) (1+(U^2)) V)/(√(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))))), V≥0}, {(1/(3 U))(-3 V-(Sqrt[3]) √(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(Sqrt[3]) √(-4+4 (U^2)+2 (V^2)-((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3))-(6 (Sqrt[3]) (1+(U^2)) V)/(√(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))))), True}})}})))^2))-(1/8) U (4+(((({{\[Piecewise], ({{(1/(3 U))(-3 V+(Sqrt[3]) √(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(Sqrt[3]) √(6 (V^2)+4 (-1+(U^2)-(V^2))-((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3))+(6 (Sqrt[3]) (1+(U^2)) V)/(√(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))))), V≥0}, {(1/(3 U))(-3 V-(Sqrt[3]) √(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(Sqrt[3]) √(-4+4 (U^2)+2 (V^2)-((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3))-(6 (Sqrt[3]) (1+(U^2)) V)/(√(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))))), True}})}})))^2))+(1/2) V (-(({{\[Piecewise], ({{(1/(3 U))(-3 V-(Sqrt[3]) √(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(Sqrt[3]) √(-4+4 (U^2)+2 (V^2)-((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3))-(6 (Sqrt[3]) (1+(U^2)) V)/(√(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))-6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))))), V≤0}, {-(1/(3 U))(3 V-(Sqrt[3]) √(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(Sqrt[3]) √(-4+4 (U^2)+2 (V^2)-((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))-(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3))+(6 (Sqrt[3]) (1+(U^2)) V)/(√(3 (V^2)+2 (-1+(U^2)-(V^2))+((((1-(U^2)+(V^2)))^2)/(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))+(((-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-16 (V^2)+(V^4))+6 (Sqrt[3]) U V (Sqrt[-(U^6)+3 (U^4) (1+(V^2))+(((1+(V^2)))^3)-3 (U^2) (1-7 (V^2)+(V^4))])))^(1/3)))))), True}})}}))+((4 V (-U-2 (Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))/(3 ((V^2)+(1/9) (((U+2 (Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2) (1+(9/(4 (((U-(Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2))))))));

(Subscript[E, c]^KS) as a function of external potential

Ecks[U_,V_]=(U (1+(((1+(4/9) (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]])))^2)))/(1+(2/9) (((U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]))^2)+(((1+(4/9) (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]])))^2))+(8 (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (((U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]))^2))/(3 (9+2 (((U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]))^2)+9 (((1+(4/9) (U-(Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]]) (U+3 V+2 (Sqrt[3+(U^2)+3 (V^2)]) Cos[(1/3) ArcCos[((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))]])))^2)))+√(1-(4 (V^2) (((U+2 (Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2))/(9 ((((V^2)+(1/9) (((U+2 (Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2) (1+(9/(4 (((U-(Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2))))))^2)))-(1/2) U (1+(4 (V^2) (((U+2 (Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2))/(9 ((((V^2)+(1/9) (((U+2 (Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2) (1+(9/(4 (((U-(Sqrt[3+(U^2)+3 (V^2)]) Sin[(1/6) (π+2 ArcCos[-((U (9+2 (U^2)-18 (V^2)))/(2 (((3+(U^2)+3 (V^2)))^(3/2))))])]))^2))))))^2)));

(Subscript[E, c,(2)]^SD) and (Subscript[W, ∞]^SD)

EcPT2[U_,x_]=-(1/256) (U^2) (((4-(x^2)))^(5/2)); (* this is exactly the same for EcGL2 and EcMP2 *)
WinfKS[U_,x_]=-(( U)/2)(((1-(Sqrt[((((x/2)))^2)])))^2);
WinfHF[U_,x_]=(1/8) (-4+(x^2)) U;

Figure 3

Show[Plot[{Ecks[0.2,VλLM[0.2,x,1]]-EcPT2[0.2,x],Ecks[1,VλLM[1,x,1]]-EcPT2[1,x],Ecks[2,VλLM[2,x,1]]-EcPT2[2,x],Ecks[5,VλLM[5,x,1]]-EcPT2[5,x],Ecks[10,VλLM[10,x,1]]-EcPT2[10,x]},{x,-2,2},Frame->True,FrameStyle->Directive[Thick,24,Black],FrameLabel->{MaTeX["\\Delta n \\equiv \\Delta n^\\text{HF}",FontSize->24],None}, PlotLegends->Placed[LineLegend[{Directive[ColorData[97,"ColorList"][[1]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[2]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[3]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[4]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[5]],AbsoluteThickness[3]]},{ "U=0.2", "U=1","U=2","U=5", "U=10"},LegendLabel-> MaTeX["E_c -E_{c, (2)}^\\text{SD} ",FontSize->24]],Right], ImageSize->500,LabelStyle->{FontFamily->"Times New Roman",FontSize->22,Black}, AxesLabel->{"Δn",None}, PlotStyle->Thickness[0.007],ScalingFunctions->"Log"],Plot[{EchfvsNhf[0.2,x]-EcPT2[0.2,x],EchfvsNhf[1,x]-EcPT2[1,x],EchfvsNhf[2,x]-EcPT2[2,x],EchfvsNhf[5,x]-EcPT2[5,x],EchfvsNhf[10,x]-EcPT2[10,x]},{x,-2,2}, Frame->True,FrameStyle->Directive[Thick,18,Black], PlotStyle->{{Dashing[0.01],Thickness[0.007]}}, ImageSize->500,ScalingFunctions->"Log"]]//Quiet

Figure 6

Show[Plot[{WinfKS[0.2,n],WinfKS[1,n],WinfKS[2,n],WinfKS[5,n],WinfKS[10,n]},{n,-2,2}, PlotRange->All,Frame->True,FrameStyle->Directive[Thick,24,Black],FrameLabel->{MaTeX["\\Delta n \\equiv \\Delta n^\\text{HF}",FontSize->24],None}, PlotLegends->Placed[LineLegend[{Directive[ColorData[97,"ColorList"][[1]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[2]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[3]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[4]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[5]],AbsoluteThickness[3]]},{ "U=0.2", "U=1","U=2","U=5", "U=10"},LegendLabel-> MaTeX["W_\\infty ",FontSize->24]],Right], ImageSize->500,LabelStyle->{FontFamily->"Times New Roman",FontSize->22,Black}, AxesLabel->{"Δn",None}, PlotStyle->Thickness[0.007]],Plot[{WinfHF[0.2,n],WinfHF[1,n],WinfHF[2,n],WinfHF[5,n],WinfHF[10,n]},{n,-2,2}, PlotRange->All, Frame->True,FrameStyle->Directive[Thick,18,Black], PlotStyle->{{Dashing[0.01],Thickness[0.007]}}, ImageSize->500]]//Quiet

Figure 7

Show[Plot[{WinfKS[0.2,nn[0.2,V]],WinfKS[1,nn[1,V]],WinfKS[2,nn[2,V]],WinfKS[5,nn[5,V]],WinfKS[10,nn[10,V]]},{V,-20,20}, PlotRange->All,Frame->True,FrameStyle->Directive[Thick,24,Black],FrameLabel->{MaTeX["\\Delta v",FontSize->24],None}, PlotLegends->Placed[LineLegend[{Directive[ColorData[97,"ColorList"][[1]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[2]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[3]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[4]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[5]],AbsoluteThickness[3]]},{ "U=0.2", "U=1","U=2","U=5", "U=10"},LegendLabel-> MaTeX["\\tilde{W}_\\infty ",FontSize->24]],Right], ImageSize->500,LabelStyle->{FontFamily->"Times New Roman",FontSize->22,Black}, PlotStyle->Thickness[0.007]],Plot[{WinfHF[0.2,NhfvsV[0.2,V]],WinfHF[1,NhfvsV[1,V]],WinfHF[2,NhfvsV[2,V]],WinfHF[5,NhfvsV[5,V]],WinfHF[10,NhfvsV[10,V]]},{V,-20,20}, PlotRange->All, Frame->True,FrameStyle->Directive[Thick,18,Black], PlotStyle->{{Dashing[0.01],Thickness[0.007]}}, ImageSize->500]]//Quiet

Comparing MP and DFT AC integrands [Figure 1]

U1V1p1=Plot[{Wλhf[1,1,λ],WλKS[1,1,λ],2 EcPT2[1.,NhfvsV[1,1]]λ,2 EcPT2[1.,nn[1,1]]λ},{λ,0,1},PlotTheme->"Scientific",BaseStyle->20,
PlotStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red},{Thickness[0.01],Blue,Dashed},{Thickness[0.01],Red,Dashed}},PlotRange->All,ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{"λ",None}, PlotLabel-> MaTeX["(e) \\quad U=1,\, |\\Delta v|=1 \\quad\\quad",FontSize->20],PlotLegends->Placed[{None,None,MaTeX["2 E_\\text{c}^\\text{MP2}\\lambda",FontSize->20],MaTeX["2 E_\\text{c}^\\text{GL2}\\lambda",FontSize->20]},{0.75,0.75}],PlotRangePadding->{0.,0.}];

U1V1p2=Plot[{Wλhf[1,1,λ],WλKS[1,1,λ],2 EcPT2[1.,NhfvsV[1,1]]λ,2 EcPT2[1.,nn[1,1]]λ},{λ,0,1},PlotTheme->"Scientific",BaseStyle->20,
PlotStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red},{Thickness[0.01],Blue,Dashed},{Thickness[0.01],Red,Dashed}},PlotRange->All,ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{"λ",(*"(Subscript[W, λ])"*)None}, PlotLabel-> MaTeX["(e) \\quad U=1,\, |\\Delta v|=1 \\quad\\quad",FontSize->20],PlotLegends->Placed[{MaTeX["W_\\lambda^\\text{HF}",FontSize->20],MaTeX["W_\\lambda^\\text{KS}",FontSize->20],None,None},{0.25,0.25}],PlotRangePadding->{0.,0.}];

Show[U1V1p1,U1V1p2]

Plot[{Wλhf[1,1,λ],WλKS[1,1,λ],WinfHF[1,NhfvsV[1,1]],WinfKS[1,nn[1,1]]},{λ,0,20},PlotTheme->"Scientific",BaseStyle->20,
PlotStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red},{Thickness[0.01],Blue,Dashed},{Thickness[0.01],Red,Dashed}},PlotRange->{-0.43,0.},ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{"λ",(*"(Subscript[W, λ])"*)None}, PlotLabel-> MaTeX["(f)\\quad U=1,\, |\\Delta v|=1 \\quad \\quad ",FontSize->20],PlotLegends->Placed[{MaTeX["W_\\lambda^\\text{HF}",FontSize->20],MaTeX["W_\\lambda^\\text{KS}",FontSize->20],MaTeX["W_\\infty^\\text{HF}",FontSize->20],MaTeX["W_\\infty^\\text{KS}",FontSize->20]},{0.75,0.55}],PlotRangePadding->{0.,0.}]

Relative error vs accuracy predictor [Figure 4]

Calculation of the tables is not immediate (takes more than few minutes)

MAPHFvsNintU1tab=ParallelTable[{N[1-((Wλhf[1,VλLM[1,x,1],1])/(2EcPT2[1,NhfvsV[1,VλLM[1,x,1]]]))],N[Abs[(((EcPT2[1,NhfvsV[1,VλLM[1,x,1]]]-Echf[1,VλLM[1,x,1]]))/(Echf[1,VλLM[1,x,1]]))]]*100},{x,0.0001,2,0.0001}];

MAPHFVSNintU1=ListPlot[{Drop[MAPHFvsNintU1tab,-5],Drop[MAPHFvsNintU1tab,-1000],Drop[MAPHFvsNintU1tab,-10000]},
PlotTheme->"Scientific",BaseStyle->16,ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{MaTeX["|1-\\lambda^\\text{HF}_\\text{ext}|",FontSize->16],MaTeX["\\text{rel. err. HF / \%}",FontSize->16]},PlotMarkers->{Graphics[Triangle[]],.03},PlotRangePadding->{0.002,0.25},PlotRange->{{0,0.35},{0,35}}, PlotLabel->MaTeX["U=1",FontSize->16]]

MAPKSvsNU1tab=ParallelTable[{N[1-((WλKSvsN[1,x,1])/(2EcPT2[1,x]))],N[Abs[(((EcPT2[1,x]-Ecks[1,VλLM[1,x,1]]))/(Ecks[1,VλLM[1,x,1]]))]]*100},{x,0.0001,2,0.0001}];

MAPKSVSNU13=ListPlot[{Drop[MAPKSvsNU1tab,-25],Drop[MAPKSvsNU1tab,-1000],Drop[MAPKSvsNU1tab,-10000]},
PlotTheme->"Scientific",BaseStyle->16,ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{MaTeX["|1-\\lambda^\\text{KS}_\\text{ext}|",FontSize->16],MaTeX["\\text{rel. err. KS / \%}",FontSize->16]},PlotMarkers->{Graphics[Triangle[]],.03},PlotRangePadding->{0.001,0.3},PlotRange->{{0,0.42},{0,45}}, PlotLabel->MaTeX["U=1",FontSize->16]]

Appendix

AC integrand [Figure 9]

Wλlb[λ_]= Winf+(1/2) (W0-Winf) ((1/(((1+((4 W0pr λ)/(5 (-W0+Winf)))))^2))+(1/(Sqrt[1+((4 W0pr λ)/(5 (-W0+Winf)))])));

lbcompx10l1=Plot[{WλhfvsNhf[1, 1.,λ],WλKSvsN[1, 1.,λ],Wλlb[λ]/.{W0pr-> 2 EcPT2[1.,1.],W0->0,Winf->WinfHF[1,1]},Wλlb[λ]/.{W0pr-> 2 EcPT2[1.,1.],W0->0,Winf-> WinfKS[1.,1]}},{λ,0,1},PlotTheme->"Scientific",BaseStyle->20,
PlotStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red},{Thickness[0.01],Blue,Dashed},{Thickness[0.01],Red,Dashed}},ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{"λ",None}, PlotLabel-> MaTeX["(iii) \\quad \\Delta n \\equiv  \\Delta n^\\text{HF} =1.0 \\quad\\quad",FontSize->20],PlotLegends->Placed[{MaTeX["W_\\lambda^\\text{HF}",FontSize->20],MaTeX["W_\\lambda^\\text{KS}",FontSize->20],None,None},{0.75,0.75}],PlotRangePadding->{0.,0.}];

lbcompx10l1bis=Plot[{WλhfvsNhf[1, 1.,λ],WλKSvsN[1, 1.,λ],Wλlb[λ]/.{W0pr-> 2 EcPT2[1.,1.],W0->0,Winf->WinfHF[1,1]},Wλlb[λ]/.{W0pr-> 2 EcPT2[1.,1.],W0->0,Winf-> WinfKS[1.,1]}},{λ,0,1},PlotTheme->"Scientific",BaseStyle->20,
PlotStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red},{Thickness[0.01],Blue,Dashed},{Thickness[0.01],Red,Dashed}},ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{"λ",None}, PlotLabel-> MaTeX["(iii) \\quad \\Delta n \\equiv  \\Delta n^\\text{HF} =1.0 \\quad\\quad",FontSize->20],PlotLegends->Placed[{None,None,MaTeX["W_{\\lambda,\\text{LB}}^\\text{HF}",FontSize->20],MaTeX["W_{\\lambda,\\text{LB}}^\\text{KS}",FontSize->20]},{0.25,0.25}],PlotRangePadding->{0.,0.}];

Show[lbcompx10l1,lbcompx10l1bis]

lbcompx10ll=Plot[{WλhfvsNhf[1, 1.,λ],WλKSvsN[1, 1.,λ],Wλlb[λ]/.{W0pr-> 2 EcPT2[1.,1.],W0->0,Winf->WinfHF[1,1]},Wλlb[λ]/.{W0pr-> 2 EcPT2[1.,1.],W0->0,Winf-> WinfKS[1.,1]}},{λ,0,60},PlotTheme->"Scientific",BaseStyle->20,
PlotStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red},{Thickness[0.01],Blue,Dashed},{Thickness[0.01],Red,Dashed}},ImageSize->400,FrameStyle->Directive[Thick,18,Black],FrameLabel->{"λ",None}, PlotLabel-> MaTeX["(iv) \\quad \\Delta n \\equiv  \\Delta n^\\text{HF} =1.0 \\quad\\quad",FontSize->20],PlotLegends->Placed[{MaTeX["W_\\lambda^\\text{HF}",FontSize->20],MaTeX["W_\\lambda^\\text{KS}",FontSize->20],MaTeX["W_{\\lambda,\\text{LB}}^\\text{HF}",FontSize->20],MaTeX["W_{\\lambda,\\text{LB}}^\\text{KS}",FontSize->20]},{0.75,0.44}],PlotRangePadding->{0.,0.}, PlotRange->{-0.39,0}]

Global Correlation energy [Figure 10]

EclbdftHD[U_,x_]=((20 U (((1-((Abs[x])/2)))^4))/(80+20 (x^2)+U (((4-(x^2)))^(5/2))-80 Abs[x]))-(1/8) U (((-2+Abs[x]))^2)+(((-10+2 (Sqrt[25+((5 U (((4-(x^2)))^(5/2)))/(4 (((-2+Abs[x]))^2)))])) (((-2+Abs[x]))^4))/(4 (((4-(x^2)))^(5/2))));
EclbhfHD[U_,x_]=(1/8) (U (-4+(x^2))-((10 U (-4+(x^2)))/(20+U (((4-(x^2)))^(3/2))))+((2 (-10+(Sqrt[5]) (Sqrt[20+U (((4-(x^2)))^(3/2))])))/(Sqrt[4-(x^2)])));

Show[Plot[{((Echf[0.2,V]-EclbhfHD[0.2,NhfvsV[0.2,V]])/(Echf[0.2,V])),((Echf[1,V]-EclbhfHD[1,NhfvsV[1,V]])/(Echf[1,V])),((Echf[2,V]-EclbhfHD[2,NhfvsV[2,V]])/(Echf[2,V])),((Echf[5,V]-EclbhfHD[5,NhfvsV[5,V]])/(Echf[5,V])),((Echf[10,V]-EclbhfHD[10,NhfvsV[10,V]])/(Echf[10,V]))},{V,-20,20}, PlotRange->All, PlotLegends->Placed[LineLegend[{Directive[ColorData[97,"ColorList"][[1]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[2]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[3]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[4]],AbsoluteThickness[3]],Directive[ColorData[97,"ColorList"][[5]],AbsoluteThickness[3]]},{ "U=0.2", "U=1","U=2","U=5", "U=10"},LegendLabel-> MaTeX["\\text{rel. err. LB}",FontSize->24]],Right], ImageSize->500,LabelStyle->{FontFamily->"Times New Roman",FontSize->22,Black},FrameLabel->{MaTeX["\\Delta v",FontSize->24],None},PlotStyle->{{Dashing[0.01],AbsoluteThickness[3]}},Frame->True,FrameStyle->Directive[Thick,24,Black]],Plot[{((Ecks[0.2,V]-EclbdftHD[0.2,nn[0.2,V]])/(Ecks[0.2,V])),((Ecks[1,V]-EclbdftHD[1,nn[1,V]])/(Ecks[1,V])),((Ecks[2,V]-EclbdftHD[2,nn[2,V]])/(Ecks[2,V])),((Ecks[5,V]-EclbdftHD[5,nn[5,V]])/(Ecks[5,V])),((Ecks[10,V]-EclbdftHD[10,nn[10,V]])/(Ecks[10,V]))},{V,-20,20}, PlotRange->All, PlotStyle->{{AbsoluteThickness[3]}}], ImageSize->500]//Quiet