Initialization

HaToeV=27.21138602;

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

PlotOptions={
AspectRatio->1,
PlotTheme->"Scientific"
};

The Hubbard dimer

These notes are based on the paper of Carrascal et al. [J. Phys.: Cond. Mat. 27 (2015) 393001]

The two views of the Hubbard dimer

Many-body view

DFT view

The two-site Hubbard model

Notations

Exact treatment

This is the exact Hamiltonian of the system expanded in the basis |12⟩,|21⟩,|11⟩,|22⟩

(({{2(Subscript[v, 1])+U, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, 2 (Subscript[v, 2])+U}}))==((Subscript[v, 1])+(Subscript[v, 2]))(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))+(({{-Δv+U, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, Δv+U}}))/.Δv->(Subscript[v, 2])-(Subscript[v, 1])/.(Subscript[v, 2])->-(Subscript[v, 1])//Simplify

The eigenstates are three singlets and a triplet state

Eigenvalues[(({{-Δv+U, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, +Δv+U}}))]//Simplify

With[{t=(1/2)},Plot[Evaluate[Table[Root[4 (t^2) U+(-4 (t^2)+(U^2)-((Δv)^2)) #1-2 U ((#1)^2)+((#1)^3)&,1],{U,{0.2,1,2,5,10}}]],{Δv,0,20},PlotRange->{All,{-10,0}},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Top}]]]

Their exact energies is complicated but can be obtained in closed form (see Appendix A)

ϵ==(Ε/(2t));	ν== ((Δv)/(2t));	u==(U/(2t));	
ϵ==(2/3)(u-w Sin[θ+(π/6)]);	w==(Sqrt[3(1+(ν^2))+(u^2)]);	Cos[3θ]==(9((ν^2)-(1/2))-(u^2))(u/(w^3));

\[LeftBracketingBar]ψ⟩==α(\[LeftBracketingBar]12⟩+\[LeftBracketingBar]21⟩)+(Subscript[β, 1])\[LeftBracketingBar]11⟩+(Subscript[β, 2])\[LeftBracketingBar]22⟩;	2(α^2)+(Subscript[β, 1]^2)+(Subscript[β, 2]^2)==1

α==c(1-(u/ϵ));	(Subscript[β, 1])==c(u-ϵ+ν);	(Subscript[β, 2])==c(u-ϵ-ν);

Δn==4(c^2)ν(ϵ-u);	(c^(-2))==2((ν^2)+(((ϵ-u))^2)(1+(ϵ^(-2))))

T==-4t α((Subscript[β, 1])+(Subscript[β, 2]))==((4(c^2)(((ϵ-u))^2))/ϵ);	(Subscript[V, ee])==U((Subscript[β, 1]^2)+(Subscript[β, 2]^2))==4(c^2)t u((((ϵ-u))^2)+(ν^2));	V==-Δv ((Subscript[β, 1]^2)-(Subscript[β, 2]^2))==((Δv Δn)/2);	Ε==T+(Subscript[V, ee])+V;

With[{t=(1/2)},Plot[Evaluate[Table[Abs[4(c^2)ν(ϵ-u)]/.c->(((2((ν^2)+(((ϵ-u))^2)(1+(ϵ^(-2))))))^(-1/2))/.ϵ->((Root[4 (t^2) U+(-4 (t^2)+(U^2)-((Δv)^2)) #1-2 U ((#1)^2)+((#1)^3)&,1])/(2t))/.{u->(U/(2t)),ν->((Δv)/(2t))},{U,{0.2,1,2,5,10}}]],{Δv,0,20},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","|Δn|"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Center}]]]

There’s a problem here as the curves reach 1 and not 2 as depicted in the review!

Non-interacting treatment

When U = 0, we have the simple tight-binding result, for which energies are

Eigenvalues[(({{-Δv, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, +Δv}}))]=={0,0,-(Sqrt[(((2t))^2)+((Δv)^2)]),+(Sqrt[(((2t))^2)+((Δv)^2)])}//Simplify

The difference in site occupation is given by

Δn==-2((Δv)/(Sqrt[(((2t))^2)+((Δv)^2)]))

Hartree-Fock treatment

The RHF energy (paramagnetic state) is

(Subscript[E, HF])==(U/2)(1-(((((Δn)/2)))^2))-(Sqrt[(((2t))^2)+(((Δv+(U/2)Δn))^2)])

The equation giving the occupation difference for the RHF case is

Δn==((-2Δv-U Δn)/(Sqrt[(((2t))^2)+(((Δv+U Δn/2))^2)]))

Show[{
With[{t=1},Plot[Evaluate[Table[(U/2)(1-(((((Δn)/2)))^2))-(Sqrt[(((2t))^2)+(((Δv+(U/2)Δn))^2)])/.Solve[Δn==((-2Δv-U Δn)/(Sqrt[(((2t))^2)+(((Δv+U Δn/2))^2)])),Δn][[3]],{U,{0.2,1,2,5,10}}]],{Δv,0,20},PlotRange->{All,{-10,5}},AspectRatio->1,PlotTheme->"Web",PlotStyle->Dashed,BaseStyle->16,FrameLabel->{"Δv","E"}]],
With[{t=1},Plot[Evaluate[Table[Root[4 (t^2) U+(-4 (t^2)+(U^2)-((Δv)^2)) #1-2 U ((#1)^2)+((#1)^3)&,1],{U,{0.2,1,2,5,10}}]],{Δv,0,20},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Top}]]]
}]

Correlation energy

This is the good, usual definition of the correlation energy.

(Subscript[E, c])=E-(E^HF)

Show[{
With[{t=(1/2)},Plot[Evaluate[Table[Root[4 (t^2) U+(-4 (t^2)+(U^2)-((Δv)^2)) #1-2 U ((#1)^2)+((#1)^3)&,1]-(U/2)(1-(((((Δn)/2)))^2))+(Sqrt[(((2t))^2)+(((Δv+(U/2)Δn))^2)])/.Solve[Δn==((-2Δv-U Δn)/(Sqrt[(((2t))^2)+(((Δv+U Δn/2))^2)])),Δn][[3]],{U,{0.2,1,2,5,10}}]],{Δv,0,20},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E"},PlotRange->All,PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Bottom}]]]
}]

Kohn-Sham density-functional theory

Generalities

Euler equation

((\[DifferentialD]F[(Subscript[n, 1])])/(\[DifferentialD](Subscript[n, 1])))-((Δv)/2)==0	\[Implies]	(Subscript[n, 1])[Δv]	\[Implies]	E[Δv]==F[(Subscript[n, 1]),Δv]+((Δv Δn)/2)

(Subscript[U, H])[Δn]==(U/2)((Subscript[n, 1]^2)+(Subscript[n, 2]^2))	(Subscript[E, x])[Δn]==-(U/4)((Subscript[n, 1]^2)+(Subscript[n, 2]^2))	(Subscript[E, Hx])[Δn]==(U/4)((Subscript[n, 1]^2)+(Subscript[n, 2]^2))==(U/2)(1+(((((Δn)/2)))^2))

(Subscript[T, s])[(Subscript[n, 1])]==-2t (Sqrt[(Subscript[n, 1])(Subscript[n, 2])])	(Subscript[E, c])[(Subscript[n, 1])]==F[(Subscript[n, 1])]-(Subscript[T, s])[(Subscript[n, 1])]-(Subscript[E, Hx])[(Subscript[n, 1])]

(Subscript[Δv, s])[Δn]==Δv+U ((Δn)/2)+(Subscript[Δv, c])[Δn]	(Subscript[Δv, c])[(Subscript[n, 1])]=-2((\[DifferentialD](Subscript[E, c])[(Subscript[n, 1])])/(\[DifferentialD](Subscript[n, 1])))

((\[DifferentialD]F[(Subscript[n, 1])])/(\[DifferentialD](Subscript[n, 1])))-((Δv)/2)==0	\[Implies]	F[(Subscript[n, 1])]==(1/2)∫Δv[(Subscript[n, 1])]\[DifferentialD](Subscript[n, 1])

Exact functional via Levy-Lieb constraint search

ϵ==(Ε/(2t)) ;	u==(U/(2t));	ν==((Δv)/(2t))

α==c(1-(u/ϵ));	(Subscript[β, 1])==c(u-ϵ+ν);	(Subscript[β, 2])==c(u-ϵ-ν);

ρ==((Abs[Δn])/2);	Δn==4(c^2)ν(ϵ-u);	(c^(-2))==2((ν^2)+(((ϵ-u))^2)(1+(ϵ^(-2))))

\[LeftBracketingBar]ψ⟩==α(\[LeftBracketingBar]12⟩+\[LeftBracketingBar]21⟩)+(Subscript[β, 1])\[LeftBracketingBar]11⟩+(Subscript[β, 2])\[LeftBracketingBar]22⟩;	2(α^2)+(Subscript[β, 1]^2)+(Subscript[β, 2]^2)==1

T==-4t α((Subscript[β, 1])+(Subscript[β, 2]));	(Subscript[V, ee])==U((Subscript[β, 1]^2)+(Subscript[β, 2]^2));	V==-Δv ((Subscript[β, 1]^2)-(Subscript[β, 2]^2));	Ε==T+(Subscript[V, ee])+V;

Ε==(Underscript[Min, α,(Subscript[β, 1]),(Subscript[β, 2])])Ε[α,(Subscript[β, 1]),(Subscript[β, 2])]

F[Δn]=(Underscript[Min, α,(Subscript[β, 1]),(Subscript[β, 2])])(T[α,(Subscript[β, 1]),(Subscript[β, 2])]+(Subscript[V, ee])[α,(Subscript[β, 1]),(Subscript[β, 2])])

Δn==2((Subscript[β, 1]^2)-(Subscript[β, 2]^2));	2(α^2)+(Subscript[β, 1]^2)+(Subscript[β, 2]^2)==1

F[Δn]=(Underscript[Min, (α^2)+(Subscript[β, 2]^2)=(1/2)(1+((\[LeftBracketingBar]Δn\[RightBracketingBar])/2))])(T[α,(Subscript[β, 2]),Δn]+(Subscript[V, ee])[α,(Subscript[β, 2]),Δn])

g==2α((((Subscript[β, 1])+β))^2)	ρ=((\[LeftBracketingBar]Δn\[RightBracketingBar])/2)

F[ρ]==(Underscript[Min, g]) f[ρ,g]	f[ρ,g]==-2t g+U(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2))))

Exact functional

The exact KS functional is defined as such

gmin[t_,U_,Δn_]:=Module[{ρ,g},
ρ=((Abs[Δn])/2);	Return[g/.N[Minimize[-2t g+U(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2)))),g],15][[2]]];
];

F[t_,U_,Δn_]:=Module[{ρ,g},
ρ=((Abs[Δn])/2);	Return[N[Minimize[-2t g+U(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2)))),g],15][[1]]];
];

We will also need its derivatives wrt Δn in order to derive the potential Δv

dFdΔn[t_,U_,Δn_]:=Module[{ρ,g0},
ρ=((Abs[Δn])/2);	g0=N[Minimize[-2t g+U(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2)))),g],15][[2]];
Return[(1/2)D[-2t g0+U(((g0^2)(1-(Sqrt[1-(g0^2)-(r^2)]))+2(r^2))/(2((g0^2)+(r^2)))),r]/.r->((Δn)/2)];
];

How does the exact functional looks like?

With[{t=(1/2)},Plot[Evaluate[Table[F[t,U,2(1-(Subscript[n, 1]))],{U,{(2/10),1,2,5,10}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,PlotRange->{{0,2},{-1,10}},FrameLabel->{"(Subscript[n, 1])","F(U,(Subscript[n, 1]))"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Center,Top}]]]

Correlation part

The correlation part looks like this

With[{t=(1/2)},Plot[Evaluate[Table[F[t,U,2(1-(Subscript[n, 1]))]+2t (Sqrt[(Subscript[n, 1])(2-(Subscript[n, 1]))])-(U/4)((Subscript[n, 1]^2)+(((2-(Subscript[n, 1])))^2)),{U,{(2/10),1,2,5,10}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Bottom}]]]

The two limiting forms of the correlation energy are

(Subscript[E, c])~-((U^2)/8)(((1-((((Subscript[n, 1])-1))^2)))^(5/2))	(U≪ 2t)		(Subscript[E, c])~-(U/2)(((1-\[LeftBracketingBar](Subscript[n, 1])-1\[RightBracketingBar]))^2)	(U≫2t)

KS potential

How does the exact KS potential looks like?

With[{t=(1/2)},Plot[Evaluate[Table[dFdΔn[t,U,2(1-(Subscript[n, 1]))],{U,{(2/10),1,2,5,10}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,PlotRange->{{0,2},{-1,10}},FrameLabel->{"(Subscript[n, 1])","((δF)/(δΔn))"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Center,Top}]]]

Small-U limit of the correlation part

Show[{
With[{t=(1/2)},Plot[Evaluate[Table[F[t,U,2(1-(Subscript[n, 1]))]+2t (Sqrt[(Subscript[n, 1])(2-(Subscript[n, 1]))])-(U/4)((Subscript[n, 1]^2)+(((2-(Subscript[n, 1])))^2)),{U,{(2/10),(1/2),1}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotLegends->Placed[{"U=0.2","U=0.5","U=1"},{Right,Bottom}]]]
,
With[{t=(1/2)},Plot[Evaluate[Table[-((U^2)/8)(((1-((((Subscript[n, 1])-1))^2)))^(5/2))(1-(U/(2t))((((Subscript[n, 1])-1))^2)(Sqrt[1-((((Subscript[n, 1])-1))^2)])),{U,{(2/10),(1/2),1}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotStyle->Dashed]]
}]

Large-U limit of the correlation part

Show[{
With[{t=(1/2)},Plot[Evaluate[Table[F[t,U,2(1-(Subscript[n, 1]))]+2t (Sqrt[(Subscript[n, 1])(2-(Subscript[n, 1]))])-(U/4)((Subscript[n, 1]^2)+(((2-(Subscript[n, 1])))^2)),{U,{5,10,20}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotLegends->Placed[{"U=5","U=10","U=20"},{Right,Bottom}]]]
,
With[{t=(1/2)},Plot[Evaluate[Table[-(U/2)(((1-Abs[(Subscript[n, 1])-1]))^2),{U,{5,10,20}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotStyle->Dashed]]
}]

Different form of the exact functional

Hopping term

F[ρ]==(Underscript[Min, g]) f[g,ρ]	(Subscript[f, 1])[g,ρ]==-2t g+2t u(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2))))

Minimizing F wrt g is equivalent to finding one of the roots of the following polynomial:

(4+(u^2))(g^6)/4+((ρ^2)(3+(u^2))-1)(g^4)+2u (ρ^2)(g^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(g^2)-2u (ρ^2)(1-(ρ^2))g-(ρ^4)(1-(ρ^2))==0

With[{t=1/2,u=1,ρ=2/10},Minimize[-2t g+2t u(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2)))),g]]//N

With[{t=1/2,u=1,ρ=2/10},NSolve[(4+(u^2))(g^6)/4+((ρ^2)(3+(u^2))-1)(g^4)+2u (ρ^2)(g^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(g^2)-2u (ρ^2)(1-(ρ^2))g-(ρ^4)(1-(ρ^2))==0,g]]//N

Electron-electron term

F[ρ]==(Underscript[Min, f]) (Subscript[f, 2])[f,ρ]	(Subscript[f, 2])[f,ρ]==-2 t(Sqrt[1-f])  ((Sqrt[f+ρ])+(Sqrt[f-ρ]))+ 2 t u f

With[{t=1/2,u=1,ρ=2/10},Minimize[-2 t(Sqrt[1-f])  ((Sqrt[f+ρ])+(Sqrt[f-ρ]))+ 2 t u f,f]]//N

With[{t=1/2,u=1,ρ=2/10},NSolve[8 (f^5) (u^2) (4+(u^2))-4 (f^6) (u^2) (4+(u^2))-4 (f^3) (u^2) (-1+(13+4 (u^2)) (ρ^2))-(ρ^2) (1+(2+8 (u^2)+4 (u^4)) (ρ^2)+(ρ^4))+4 (f^4) (u^2) (-5+7 (ρ^2)+(u^2) (-1+2 (ρ^2)))-4 (f^2) (ρ^2) (1+(u^4) (-2+(ρ^2))+(u^2) (-7+3 (ρ^2)))+4 f (ρ^2) (1+(ρ^2)+(u^2) (-1+(5+2 (u^2)) (ρ^2)))==0,f]]//N

Asymmetric limit

F[ρ]==(Underscript[Min, l]) (Subscript[f, 3])[l,ρ]	(Subscript[f, 3])[l,ρ]==-2 t (Sqrt[2l-(l^2)-(ρ^2)])+ 2u t (((l^2)+(ρ^2))/(2l))

With[{t=1/2,u=1,ρ=2/10},Minimize[-2 t (Sqrt[2l-(l^2)-(ρ^2)])+ 2u t (((l^2)+(ρ^2))/(2l)),l]]//N

With[{t=1/2,u=1,ρ=2/10},NSolve[(l^6) (4+(u^2))-2 (l^5) (4+(u^2))+(l^4) (4-(u^2) (ρ^2))+4 (l^3) (u^2) (ρ^2)-(l^2) (u^2) (ρ^4)-2 l (u^2) (ρ^4)+(u^2) (ρ^6)==0,l]]//N

Self-consistent Kohn-Sham module

KS[t_,U_,Δv_]:=Module[{Ts,EHx,EKS,Vext,Δn,λ=0,δ=(10^(-4))},
Δn=Chop[Δn/.FindRoot[((t/2)/(Sqrt[1-(((((Δn)/2)))^2)]))Δn+(U/4)Δn+λ ((F0[t,U,Δn+δ]-F0[t,U,Δn-δ])/(2δ))+((Δv)/2),{Δn,0.9}]];
Print[Δn];
Ts=-2t(Sqrt[(1-((Δn)/2)) (1+((Δn)/2))]);	EHx=(U/2)(1+(((((Δn)/2)))^2));	Vext=((Δv Δn)/2);
EKS=Ts+EHx+Vext+λ F0[t,U,Δn];
Return[Chop[EKS]];
];

With[{t=1,U=1/2,Δv=-1},{KS[t,U,Δv],Root[4 (t^2) U+(-4 (t^2)+(U^2)-((Δv)^2)) #1-2 U ((#1)^2)+((#1)^3)&,1],(U/2)(1-(((((Δn)/2)))^2))-(Sqrt[(((2t))^2)+(((Δv+(U/2)Δn))^2)])/.Solve[Δn==((-2Δv-U Δn)/(Sqrt[(((2t))^2)+(((Δv+U Δn/2))^2)])),Δn][[1]]}]//N

There’s still something wrong with this module!! I’m working on it.

Fundamental gap

Generalities

I=E[N-1]-E[N]	A=E[N]-E[N+1]	(Subscript[E, g])=I-A

(ϵ^HOMO)=-I	(ϵ^LUMO)≠-A	(Subscript[E, gs])=(ϵ^LUMO)-(ϵ^HOMO)	(Subscript[E, g])=(Subscript[E, gs])+(Subscript[Δ, xc])

(Subscript[Δ, xc]) is the infamous derivative discontinuity

Energies of the 1-, 2- and 3-electron systems

One-electron system

Configurations: |1↑⟩, |2↑⟩

(Subscript[ϵ, 1])[t_,U_,Δv_]=Eigenvalues[(({{-((Δv)/2), -t}, {-t, ((Δv)/2)}}))]//FullSimplify

With[{t=(1/2),U=1},Plot[Evaluate[(Subscript[ϵ, 1])[t,U,Δv]],{Δv,0,3},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E(N=1)"},PlotLegends->Placed[Automatic,{Right,Center}]]]

Two-electron system

Configurations: |1↑1↓⟩, |1↑2↓⟩,|1↓2↑⟩,|2↑2↓⟩

(Subscript[ϵ, 2])[t_,U_,Δv_]=Eigenvalues[(({{-Δv+U, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, Δv+U}}))]//FullSimplify

With[{t=(1/2),U=1},Plot[Evaluate[(Subscript[ϵ, 2])[t,U,Δv]],{Δv,0,3},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E(N=2)"},PlotLegends->Placed[Automatic,{Right,Center}]]]

Three-electron system

Configurations: |1↑1↓2↑⟩, |1↑2↑2↓⟩

(Subscript[ϵ, 3])[t_,U_,Δv_]=Eigenvalues[(({{-((Δv)/2)+U, -t}, {-t, ((Δv)/2)+U}}))]//FullSimplify

With[{t=(1/2),U=1},Plot[Evaluate[(Subscript[ϵ, 3])[t,U,Δv]],{Δv,0,3},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E(N=3)"},PlotLegends->Placed[Automatic,{Right,Center}]]]

IP & EA

With[{t=(1/2),U=1},Plot[{Evaluate[(Subscript[ϵ, 2])[t,U,Δv][[2]]-(Subscript[ϵ, 1])[t,U,Δv][[1]]],Evaluate[(Subscript[ϵ, 3])[t,U,Δv][[1]]-(Subscript[ϵ, 2])[t,U,Δv][[2]]]},{Δv,0,3},PlotRange->{{0,3},{-1,2}},AspectRatio->1,PlotTheme->"Business",BaseStyle->16,FrameLabel->{"Δv",None},PlotLegends->Placed[{"-I","-A"},{Right,Center}]]]

With[{t=(1/2),U=5},Plot[{Evaluate[(Subscript[ϵ, 2])[t,U,Δv][[2]]-(Subscript[ϵ, 1])[t,U,Δv][[1]]],Evaluate[(Subscript[ϵ, 3])[t,U,Δv][[1]]-(Subscript[ϵ, 2])[t,U,Δv][[2]]]},{Δv,0,10},PlotRange->{{0,10},{-2,6}},AspectRatio->1,PlotTheme->"Business",BaseStyle->16,FrameLabel->{"Δv",None},PlotLegends->Placed[{"-I","-A"},{Right,Center}]]]

Correlation

Generalities

(Subscript[E, c])=(Subscript[T, c])+(Subscript[U, c])	(Subscript[T, c])=T-(Subscript[T, s])	(Subscript[U, c])=(Subscript[V, ee])-(Subscript[E, Hx])	(Subscript[U, x])=(Subscript[E, x])	(Subscript[E, c])<0	(Subscript[T, c])>0	(Subscript[U, c])<0

Δn==4(c^2)ν(ϵ-u);	(c^(-2))==2((ν^2)+(((ϵ-u))^2)(1+(ϵ^(-2))))

T==((4(c^2)(((ϵ-u))^2))/ϵ)	(Subscript[V, ee])==4(c^2)t u((((ϵ-u))^2)+(ν^2))	V==((Δv Δn)/2)	E==T+(Subscript[V, ee])+V=F+V

(Subscript[T, s])==-2t(Sqrt[(Subscript[n, 1])(Subscript[n, 2])])	(Subscript[U, H])[Δn]==(U/2)((Subscript[n, 1]^2)+(Subscript[n, 2]^2))	(Subscript[E, x])[Δn]==-(U/4)((Subscript[n, 1]^2)+(Subscript[n, 2]^2))	(Subscript[E, Hx])[Δn]==(U/4)((Subscript[n, 1]^2)+(Subscript[n, 2]^2))==(U/2)(1+(((((Δn)/2)))^2))

((\[DifferentialD]F[(Subscript[n, 1])])/(\[DifferentialD](Subscript[n, 1])))-((Δv)/2)==0	\[Implies]	Δv=2((\[DifferentialD]F[(Subscript[n, 1])])/(\[DifferentialD](Subscript[n, 1])))	((\[DifferentialD]F)/(\[DifferentialD](Subscript[n, 1])))=((\[DifferentialD]Δn)/(\[DifferentialD](Subscript[n, 1])))((\[DifferentialD]F)/(\[DifferentialD]Δn))=-2((\[DifferentialD]F)/(\[DifferentialD]Δn))	Δv=-4((\[DifferentialD]F[Δn])/(\[DifferentialD]Δn))

Accurate parametrization of correlation energy

“Good” functional

F0[t_,U_,Δn_]:=Module[{u,ρ,a1,a2,a11,a12,a21,a22,g0},
u=(U/(2t));	ρ=((Abs[Δn])/2);
a21=(1/2)(Sqrt[(1-ρ)ρ/2]);	a11=a21(1+(ρ^(-1)));
a12=(1/2)(1-ρ);	a22=a12/2;
a1=a11+u a12;	a2=a21+u a22;
g0=(Sqrt[(((1-ρ)(1+ρ(1+(((1+ρ))^3)u a1)))/(1+(((1+ρ))^3)u a2))]);
Return[-2t g0+U(((g0^2)(1-(Sqrt[1-(g0^2)-(ρ^2)]))+2(ρ^2))/(2((g0^2)+(ρ^2))))];
];

With[{t=(1/2)},Plot[Evaluate[Table[((F0[t,U,2(1-(Subscript[n, 1]))]-F[t,U,2(1-(Subscript[n, 1]))])/U),{U,{(2/10),1,2,5,10}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,PlotRange->All,FrameLabel->{"(Subscript[n, 1])","(((Subscript[F, 0]) (U, (Subscript[n, 1])) - F (U, (Subscript[n, 1])))/U)"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Top}]]]

“Very good” functional

F1[t_,U_,Δn_]:=Module[{u,ρ,a1,a2,a11,a12,a21,a22,g0,dhdg,dgdλ,g1},
u=(U/(2t));	ρ=((Δn)/2);
a21=(1/2)(Sqrt[(1-ρ)ρ/2]);	a11=a21(1+(ρ^(-1)));
a12=(1/2)(1-ρ);	a22=a12/2;
a1=a11+u a12;	a2=a21+u a22;
g0=(Sqrt[(((1-ρ)(1+ρ(1+(((1+ρ))^3)u a1)))/(1+(((1+ρ))^3)u a2))]);
dhdg=D[(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2)))),g]/.g->g0;
dgdλ=(((1-ρ)(((1+ρ))^3)(u^2)((3ρ/2-1+ρ(((1+ρ))^3)u a2)a12-ρ(1+(((1+ρ))^3)u a1)a22))/(2g0(((1+(((1+ρ))^3)u a2))^2)));
g1=g0+(u dhdg-1)dgdλ;
Return[-2t g1+U(((g1^2)(1-(Sqrt[1-(g1^2)-(ρ^2)]))+2(ρ^2))/(2((g1^2)+(ρ^2))))];
]

dF1dΔn[t_,U_,Δn_]:=Evaluate[(∂_(Δn))F1[t,U,Δn]]

With[{t=(1/2)},Plot[Evaluate[Table[((F1[t,U,2(1-(Subscript[n, 1]))]-F[t,U,2(1-(Subscript[n, 1]))])/U),{U,{(2/10),1,2,5,10}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,PlotRange->All,FrameLabel->{"(Subscript[n, 1])","(((Subscript[F, 1]) (U, (Subscript[n, 1])) - F (U, (Subscript[n, 1])))/U)"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Top}]]]

Hartree-Fock

Restricted Hartree-Fock (paramagnetic state)

The one-electron effective potential is

(Subscript[v, i]^RHF)[(Subscript[n, i])]==(Subscript[v, i])+U (Subscript[n, i])/2

and the RHF eigenvalues are

(ϵ^RHF)=((U-(Sqrt[(((((Δv)^HF)))^2)+(((2t))^2)]))/2)

The RHF energy is

(Subscript[E, HF])==(U/2)(1-(((((Δn)/2)))^2))-(Sqrt[(((2t))^2)+(((Δv+(U/2)Δn))^2)])

The equation giving the occupation difference for the RHF case is

Δn==((-2Δv-U Δn)/(Sqrt[(((2t))^2)+(((Δv+U Δn/2))^2)]))

Unrestricted Hartree-Fock

σ==+1== ↑	σ==-1== ↓	(Overscript[σ, -])=-σ

(Subscript[n, i])=(Subscript[n, i,↑])+(Subscript[n, i,↓])	(* Occupation *)	(Subscript[m, i])=(Subscript[n, i,↑])-(Subscript[n, i,↓])	(* Spin occupation *)

(Subscript[v, iσ]^HF)==(Subscript[v, i])+U (Subscript[n, i(Overscript[σ, –])])	\[Implies]	(Subscript[Δv, σ]^HF)=Δv+(U/2)(Δn-σ Δm)	(Subscript[t, σ]^HF)=(Sqrt[(((2t))^2)+((((Subscript[Δv, σ]^HF)))^2)])	(* this formula has been corrected *)

(Subscript[ϵ, ±,σ]^HF)==(U/4)(N-σ M)±((Subscript[t, (Overscript[σ, –])]^HF)/2)	E==(U/2)(1-((((Δn)^2)-((Δm)^2))/4))-(((Subscript[t, ↑]^HF)+(Subscript[t, ↓]^HF))/2)	Δn==-(∑_σ)((Subscript[Δv, σ]^HF)/(Subscript[t, σ]^HF))	Δm==-(∑_σ)σ((Subscript[Δv, σ]^HF)/(Subscript[t, σ]^HF))

Δn==-((Subscript[Δv, ↑]^HF)/(Subscript[t, ↑]^HF))-((Subscript[Δv, ↓]^HF)/(Subscript[t, ↓]^HF))==-((Subscript[Δv, ↑]^HF)/(Sqrt[(((2t))^2)+((((Subscript[Δv, ↑]^HF)))^2)]))-((Subscript[Δv, ↓]^HF)/(Sqrt[(((2t))^2)+((((Subscript[Δv, ↓]^HF)))^2)]))==-((Δv+(U/2)(Δn-Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn-Δm)))^2)]))-((Δv+(U/2)(Δn+Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn+Δm)))^2)]))
Δm==-((Subscript[Δv, ↑]^HF)/(Subscript[t, ↑]^HF))+((Subscript[Δv, ↓]^HF)/(Subscript[t, ↓]^HF))==-((Subscript[Δv, ↑]^HF)/(Sqrt[(((2t))^2)+((((Subscript[Δv, ↑]^HF)))^2)]))+((Subscript[Δv, ↓]^HF)/(Sqrt[(((2t))^2)+((((Subscript[Δv, ↓]^HF)))^2)]))==-((Δv+(U/2)(Δn-Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn-Δm)))^2)]))+((Δv+(U/2)(Δn+Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn+Δm)))^2)]))

Δm=0	=>	paramagnetic state i.e. (singlet) RHF state
Δm≠0	=>	ferromagnetic or antiferromagnetic states i.e. (singlet or triplet) UHF states

Show[{
With[{t=(1/2),U=5},Plot[Abs[Δn/.Solve[(((((Δn)/2)))^2)==((((Δv+(U/2) Δn))^2)/((((2t))^2)+(((Δv+U Δn/2))^2))),Δn][[3]]],{Δv,0,10},PlotRange->{All,{0,2}},AspectRatio->1,PlotTheme->"Web",PlotStyle->Dashed,BaseStyle->16,FrameLabel->{"Δv","|Δn|"}]],
With[{t=(1/2),U=5},Plot[Abs[Δn/.NSolve[{Δn==-((Δv+(U/2)(Δn-Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn-Δm)))^2)]))-((Δv+(U/2)(Δn+Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn+Δm)))^2)])),Δm==-((Δv+(U/2)(Δn-Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn-Δm)))^2)]))+((Δv+(U/2)(Δn+Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn+Δm)))^2)]))},{Δn,Δm}]//Last],{Δv,0,10},PlotRange->{All,{0,2}},AspectRatio->1,PlotTheme->"Web",PlotStyle->DotDashed,BaseStyle->16,FrameLabel->{"Δv","|Δn|"}]],
With[{t=(1/2),U=5},Plot[2Abs[4(c^2)ν(ϵ-u)]/.c->(((2((ν^2)+(((ϵ-u))^2)(1+(ϵ^(-2))))))^(-1/2))/.ϵ->((Root[4 (t^2) U+(-4 (t^2)+(U^2)-((Δv)^2)) #1-2 U ((#1)^2)+((#1)^3)&,1])/(2t))/.{u->(U/(2t)),ν->((Δv)/(2t))},{Δv,0,10},PlotRange->{All,{0,2}},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","|Δn|"}]]
}]

Show[{
With[{t=1},Plot[Evaluate[Table[(U/2)(1-(((((Δn)/2)))^2))-(Sqrt[(((2t))^2)+(((Δv+(U/2)Δn))^2)])/.Solve[(((((Δn)/2)))^2)==((((Δv+(U/2) Δn))^2)/((((2t))^2)+(((Δv+U Δn/2))^2))),Δn][[3]],{U,{0.2,1,2,5,10}}]],{Δv,0,20},PlotRange->{All,{-10,5}},AspectRatio->1,PlotTheme->"Web",PlotStyle->Dashed,BaseStyle->16,FrameLabel->{"Δv","E"}]],
With[{t=1},Plot[Evaluate[Table[(U/2)(1-((((Δn)^2)-((Δm)^2))/4))-(((Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn-Δm)))^2)])+(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn+Δm)))^2)]))/2)/.Solve[{Δn==-((Δv+(U/2)(Δn-Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn-Δm)))^2)]))-((Δv+(U/2)(Δn+Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn+Δm)))^2)])),Δm==-((Δv+(U/2)(Δn-Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn-Δm)))^2)]))+((Δv+(U/2)(Δn+Δm))/(Sqrt[(((2t))^2)+(((Δv+(U/2)(Δn+Δm)))^2)]))},{Δn,Δm}][[3]],{U,{0.2,1,2,5,10}}]],{Δv,0,20},PlotRange->{All,{-10,5}},AspectRatio->1,PlotTheme->"Web",PlotStyle->DotDashed,BaseStyle->16,FrameLabel->{"Δv","E"}]],
With[{t=1},Plot[Evaluate[Table[Root[4 (t^2) U+(-4 (t^2)+(U^2)-((Δv)^2)) #1-2 U ((#1)^2)+((#1)^3)&,1],{U,{0.2,1,2,5,10}}]],{Δv,0,20},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E"},PlotLegends->Placed[{"U=0.2","U=1","U=2","U=5","U=10"},{Right,Top}]]]
}]

There’s still something wrong with this part.

Fractional particle number

Many limits of the exact functional

Small-U limit (work in progress)

(Subscript[g, n_])[ρ_,u_]:=(∑_(k=0)^n)(Subscript[c, k])((u^k)/(k!))

(4+(u^2))(g^6)/4+((ρ^2)(3+(u^2))-1)(g^4)+2u (ρ^2)(g^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(g^2)-2u (ρ^2)(1-(ρ^2))g-(ρ^4)(1-(ρ^2))==0

With[{n=0},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0,(Subscript[c, n])]//Last
]

With[{n=1},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.(Subscript[c, 0])->(Sqrt[1-(ρ^2)])//FullSimplify,(Subscript[c, n])]//Last
]

With[{n=2},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0}//FullSimplify,(Subscript[c, n])]
]

With[{n=3},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)])))}//FullSimplify,(Subscript[c, n])]
]

With[{n=4},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3)}//FullSimplify,(Subscript[c, n])]
]

With[{n=5},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)])))}//FullSimplify,(Subscript[c, n])]
]

With[{n=6},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2))}//FullSimplify,(Subscript[c, n])]
]

With[{n=7},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)])))}//FullSimplify,(Subscript[c, n])]
]

(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))

With[{n=8},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)]))),(Subscript[c, 7])->-945 (ρ^4) (((-1+(ρ^2)))^5) (-1+24 (ρ^2)-70 (ρ^4)+51 (ρ^6))}//FullSimplify,(Subscript[c, n])]
]

{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)]))),(Subscript[c, 7])->-945 (ρ^4) (((-1+(ρ^2)))^5) (-1+24 (ρ^2)-70 (ρ^4)+51 (ρ^6)),-((2205 (((-1+(ρ^2)))^6) (-5-33 (ρ^2)+65 (ρ^4)-16875 (ρ^6)+272 (ρ^8) (419-798 (ρ^2)+456 (ρ^4))))/(256 (Sqrt[1-(ρ^2)])))}

Evaluate[Table[(Subscript[c, n]),{n,0,8}]]/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)]))),(Subscript[c, 7])->-945 (ρ^4) (((-1+(ρ^2)))^5) (-1+24 (ρ^2)-70 (ρ^4)+51 (ρ^6)),(Subscript[c, 8])->-((2205 (((-1+(ρ^2)))^6) (-5-33 (ρ^2)+65 (ρ^4)-16875 (ρ^6)+272 (ρ^8) (419-798 (ρ^2)+456 (ρ^4))))/(256 (Sqrt[1-(ρ^2)])))}//FullSimplify

Coef={(Sqrt[1-(ρ^2)]),0,-(1/4) (((1-(ρ^2)))^(5/2)),(3/2)  (ρ^2)(((1-(ρ^2)))^3),(9/16) (((1-(ρ^2)))^(7/2)) (1+7 (ρ^2)-24 (ρ^4)),15 (ρ^4)(((1-(ρ^2)))^4) (-7+11 (ρ^2)),-(225/64) (((1-(ρ^2)))^(9/2)) (1+6 (ρ^2)+121 (ρ^4)-728 (ρ^6)+728 (ρ^8)),-945 (ρ^4) (((1-(ρ^2)))^5) (1-24 (ρ^2)+70 (ρ^4)-51 (ρ^6)),(2205/256) (((1-(ρ^2)))^(11/2)) (5+33 (ρ^2)-65 (ρ^4)+16875 (ρ^6)-113968 (ρ^8)+217056 (ρ^10)-124032 (ρ^12))}//FullSimplify

(∑_(k=0)^8)Coef[[k+1]]((u^k)/(k!))

Show[{
With[{t=(1/2)},Plot[Evaluate[Table[gmin[t,U,2(1-(Subscript[n, 1]))],{U,{5}}]],{(Subscript[n, 1]),0,2},AspectRatio->1/3,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotLegends->Placed[{"U=0.2","U=1","U=2"},{Right,Bottom}]]]
,
With[{t=(1/2)},Plot[Evaluate[Table[PadeApproximant[(∑_(k=0)^8)Coef[[k+1]]((u^k)/(k!)),{u,0,{2,2}}]/.ρ->1-(Subscript[n, 1])/.u->(U/(2t)),{U,{5}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotStyle->Dashed]]
}]

Show[{
With[{t=(1/2)},Plot[Evaluate[Table[gmin[t,U,2(1-(Subscript[n, 1]))],{U,{2}}]],{(Subscript[n, 1]),0,2},AspectRatio->1/3,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotLegends->Placed[{"U=0.2","U=1","U=2"},{Right,Bottom}]]]
,
With[{t=(1/2)},Plot[Evaluate[Table[(∑_(k=0)^8)Coef[[k+1]]((u^k)/(k!))/.ρ->1-(Subscript[n, 1])/.u->(U/(2t)),{U,{2}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotStyle->Dashed]]
}]

Small-ρ limit (work in progress)

(Subscript[g, n_])[ρ_,u_]:=(∑_(k=0)^n)(Subscript[c, k])((ρ^k)/(k!))

(4+(u^2))(g^6)/4+((ρ^2)(3+(u^2))-1)(g^4)+2u (ρ^2)(g^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(g^2)-2u (ρ^2)(1-(ρ^2))g-(ρ^4)(1-(ρ^2))==0

(4+(u^2))(g^6)/4+((ρ^2)(3+(u^2))-1)(g^4)+2u (ρ^2)(g^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(g^2)-2u (ρ^2)(1-(ρ^2))g-(ρ^4)(1-(ρ^2))/.u->2 (Sqrt[-1+(r^2)])//FullSimplify

With[{n=0},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,ρ][[n+1]]==0
//FullSimplify,(Subscript[c, n])]
]

With[{n=1},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,ρ][[n+1]]==0/.(Subscript[c, 0])->(2/(Sqrt[4+(u^2)]))//FullSimplify,(Subscript[c, n])]
]

With[{n=2},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,ρ][[n+1]]==0/.{(Subscript[c, 0])->(2/(Sqrt[4+(u^2)])),(Subscript[c, 1])->0}//FullSimplify,(Subscript[c, n])]//FullSimplify
]

With[{n=3},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,ρ][[n+1]]==0/.{(Subscript[c, 0])->(2/(Sqrt[4+(u^2)])),(Subscript[c, 1])->0,(Subscript[c, 2])->((-4+(u^2) (2+u (u+(Sqrt[4+(u^2)]))))/(2 (Sqrt[4+(u^2)])))}//FullSimplify,(Subscript[c, n])]
]

With[{n=4},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,ρ][[n+1]]==0/.{(Subscript[c, 0])->(2/(Sqrt[4+(u^2)])),(Subscript[c, 1])->0,(Subscript[c, 2])->((-4+(u^2) (2+u (u+(Sqrt[4+(u^2)]))))/(2 (Sqrt[4+(u^2)]))),(Subscript[c, 3])->0}//FullSimplify,(Subscript[c, n])]
]

With[{n=5},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,ρ][[n+1]]==0/.{(Subscript[c, 0])->(2/(Sqrt[4+(u^2)])),(Subscript[c, 1])->0,(Subscript[c, 2])->((-4+(u^2) (2+u (u+(Sqrt[4+(u^2)]))))/(2 (Sqrt[4+(u^2)]))),(Subscript[c, 3])->0,(Subscript[c, 4])->-((3 (4+8 (u^2)+34 (u^4)+20 (u^6)+3 (u^8)+12 (u^3) (Sqrt[4+(u^2)])+14 (u^5) (Sqrt[4+(u^2)])+3 (u^7) (Sqrt[4+(u^2)])))/(2 (Sqrt[4+(u^2)])))}//FullSimplify,(Subscript[c, n])]
]

With[{n=6},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2))}//FullSimplify,(Subscript[c, n])]
]

With[{n=7},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)])))}//FullSimplify,(Subscript[c, n])]
]

(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))

With[{n=8},

Solve[CoefficientList[(4+(u^2))(((Subscript[g, n])[ρ,u])^6)/4+((ρ^2)(3+(u^2))-1)(((Subscript[g, n])[ρ,u])^4)+2u (ρ^2)(((Subscript[g, n])[ρ,u])^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(((Subscript[g, n])[ρ,u])^2)-2u (ρ^2)(1-(ρ^2))(Subscript[g, n])[ρ,u]-(ρ^4)(1-(ρ^2))//Expand,u][[n+1]]==0/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)]))),(Subscript[c, 7])->-945 (ρ^4) (((-1+(ρ^2)))^5) (-1+24 (ρ^2)-70 (ρ^4)+51 (ρ^6))}//FullSimplify,(Subscript[c, n])]
]

{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)]))),(Subscript[c, 7])->-945 (ρ^4) (((-1+(ρ^2)))^5) (-1+24 (ρ^2)-70 (ρ^4)+51 (ρ^6)),-((2205 (((-1+(ρ^2)))^6) (-5-33 (ρ^2)+65 (ρ^4)-16875 (ρ^6)+272 (ρ^8) (419-798 (ρ^2)+456 (ρ^4))))/(256 (Sqrt[1-(ρ^2)])))}

Evaluate[Table[(Subscript[c, n]),{n,0,8}]]/.{(Subscript[c, 0])->(Sqrt[1-(ρ^2)]),(Subscript[c, 1])->0,(Subscript[c, 2])->((((-1+(ρ^2)))^3)/(4 (Sqrt[1-(ρ^2)]))),(Subscript[c, 3])->-(3/2) (ρ^2) (((-1+(ρ^2)))^3),(Subscript[c, 4])->-((9 (((-1+(ρ^2)))^4) (-1-7 (ρ^2)+24 (ρ^4)))/(16 (Sqrt[1-(ρ^2)]))),(Subscript[c, 5])->15 (ρ^4) (((-1+(ρ^2)))^4) (-7+11 (ρ^2)),(Subscript[c, 6])->((225 (((-1+(ρ^2)))^5) (1+6 (ρ^2)+121 (ρ^4)+728 (ρ^6) (-1+(ρ^2))))/(64 (Sqrt[1-(ρ^2)]))),(Subscript[c, 7])->-945 (ρ^4) (((-1+(ρ^2)))^5) (-1+24 (ρ^2)-70 (ρ^4)+51 (ρ^6)),(Subscript[c, 8])->-((2205 (((-1+(ρ^2)))^6) (-5-33 (ρ^2)+65 (ρ^4)-16875 (ρ^6)+272 (ρ^8) (419-798 (ρ^2)+456 (ρ^4))))/(256 (Sqrt[1-(ρ^2)])))}//FullSimplify

Coef={(Sqrt[1-(ρ^2)]),0,-(1/4) (((1-(ρ^2)))^(5/2)),(3/2)  (ρ^2)(((1-(ρ^2)))^3),(9/16) (((1-(ρ^2)))^(7/2)) (1+7 (ρ^2)-24 (ρ^4)),15 (ρ^4)(((1-(ρ^2)))^4) (-7+11 (ρ^2)),-(225/64) (((1-(ρ^2)))^(9/2)) (1+6 (ρ^2)+121 (ρ^4)-728 (ρ^6)+728 (ρ^8)),-945 (ρ^4) (((1-(ρ^2)))^5) (1-24 (ρ^2)+70 (ρ^4)-51 (ρ^6)),(2205/256) (((1-(ρ^2)))^(11/2)) (5+33 (ρ^2)-65 (ρ^4)+16875 (ρ^6)-113968 (ρ^8)+217056 (ρ^10)-124032 (ρ^12))}//FullSimplify

(∑_(k=0)^8)Coef[[k+1]]((u^k)/(k!))

Show[{
With[{t=(1/2)},Plot[Evaluate[Table[gmin[t,U,2(1-(Subscript[n, 1]))],{U,{5}}]],{(Subscript[n, 1]),0,2},AspectRatio->1/3,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotLegends->Placed[{"U=0.2","U=1","U=2"},{Right,Bottom}]]]
,
With[{t=(1/2)},Plot[Evaluate[Table[PadeApproximant[(∑_(k=0)^8)Coef[[k+1]]((u^k)/(k!)),{u,0,{2,2}}]/.ρ->1-(Subscript[n, 1])/.u->(U/(2t)),{U,{5}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotStyle->Dashed]]
}]

Show[{
With[{t=(1/2)},Plot[Evaluate[Table[gmin[t,U,2(1-(Subscript[n, 1]))],{U,{2}}]],{(Subscript[n, 1]),0,2},AspectRatio->1/3,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotLegends->Placed[{"U=0.2","U=1","U=2"},{Right,Bottom}]]]
,
With[{t=(1/2)},Plot[Evaluate[Table[(∑_(k=0)^8)Coef[[k+1]]((u^k)/(k!))/.ρ->1-(Subscript[n, 1])/.u->(U/(2t)),{U,{2}}]],{(Subscript[n, 1]),0,2},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"(Subscript[n, 1])","(Subscript[E, c])"},PlotStyle->Dashed]]
}]

Can I find the exact functional in closed form?

ϵ==(Ε/(2t)) ;	u==(U/(2t));	ν==((Δv)/(2t))

α==c(1-(u/ϵ));	(Subscript[β, 1])==c(u-ϵ+ν);	(Subscript[β, 2])==c(u-ϵ-ν);

ρ==((Abs[Δn])/2);	Δn==4(c^2)ν(ϵ-u);	(c^(-2))==2((ν^2)+(((ϵ-u))^2)(1+(ϵ^(-2))))

\[LeftBracketingBar]ψ⟩==α(\[LeftBracketingBar]12⟩+\[LeftBracketingBar]21⟩)+(Subscript[β, 1])\[LeftBracketingBar]11⟩+(Subscript[β, 2])\[LeftBracketingBar]22⟩;	2(α^2)+(Subscript[β, 1]^2)+(Subscript[β, 2]^2)==1

T==-4t α((Subscript[β, 1])+(Subscript[β, 2]));	(Subscript[V, ee])==U((Subscript[β, 1]^2)+(Subscript[β, 2]^2));	V==-Δv ((Subscript[β, 1]^2)-(Subscript[β, 2]^2));	Ε==T+(Subscript[V, ee])+V;

Ε==(Underscript[Min, α,(Subscript[β, 1]),(Subscript[β, 2])])Ε[α,(Subscript[β, 1]),(Subscript[β, 2])]

Δn==2((Subscript[β, 1]^2)-(Subscript[β, 2]^2));	2(α^2)+(Subscript[β, 1]^2)+(Subscript[β, 2]^2)==1

F[Δn]=(Underscript[Min, α,(Subscript[β, 1]),(Subscript[β, 2])])(T[α,(Subscript[β, 1]),(Subscript[β, 2])]+(Subscript[V, ee])[α,(Subscript[β, 1]),(Subscript[β, 2])])

Hopping term parametrization

f[g_]:=-2 t g+ 2u t(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2))))

min=Table[
res=Minimize[f[g]/.{t->1,u->1},g]//N;
{g/.res[[2]],res[[1]]//N},{ρ,0,1,1/50}]

Levy=Table[
res=Minimize[f[g]/.{t->1,u->1},g]//N;
{ρ,res[[1]]},{ρ,0,1,1/50}]

Show[
Flatten[
{
With[{t=1,u=1},
Plot[Evaluate[Table[-2 t g+ 2u t(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2))))
,{ρ,-1,1,1/20}]],{g,0,1},PlotRange->{-2,2},ImageSize->Large
]
]
,
ListPlot[min,PlotStyle->Black,PlotMarkers->{Automatic,Medium},Joined->True],
ListPlot[Levy,PlotStyle->Black,PlotMarkers->{Automatic,Medium},Joined->True],
Plot[Root[-432 (ρ^2)+864 (ρ^4)+80 (ρ^6)-768 (ρ^4) #1+(240 (ρ^2)+144 (ρ^4)) ((#1)^2)+(8-72 (ρ^2)) ((#1)^3)+4 (ρ^2) ((#1)^4)-4 ((#1)^5)+((#1)^6)&,1],{ρ,0,1},PlotStyle->{Thick,Red}]
}]
]

Asymmetric limit parametrization

f[l_]=-2 t (Sqrt[2l-(l^2)-(ρ^2)])+ 2u t (((l^2)+(ρ^2))/(2l))/.{t->1,u->1};

min=Table[
res=Minimize[f[l]/.{t->1,u->1},l]//N;
{l/.res[[2]],res[[1]]},{ρ,0,1,1/50}]

Show[
Flatten[
{
With[{t=1,u=1},
Plot[Evaluate[Table[f[l]
,{ρ,0,1,1/50}]],{l,0,1},PlotRange->{-2,2},ImageSize->Large
]
]
,
ListPlot[min,PlotStyle->Black,PlotMarkers->{Automatic,Medium}]
,
Plot[Root[-432 (ρ^2)+864 (ρ^4)+80 (ρ^6)-768 (ρ^4) #1+(240 (ρ^2)+144 (ρ^4)) ((#1)^2)+(8-72 (ρ^2)) ((#1)^3)+4 (ρ^2) ((#1)^4)-4 ((#1)^5)+((#1)^6)&,1],{ρ,0,1},PlotStyle->{Thick,Red}]
}]
]

Minimize[f[l],l]//FullSimplify

electron-electron term parametrization

Let’s switch to spherical coordinates which naturally enforce normalization

2(α^2)+(Subscript[β, 1]^2)+(Subscript[β, 2]^2)==1	⇔	(x^2)+(y^2)+(z^2)==1	⇔	x==Cos[ϕ]Sin[θ]	y==Sin[ϕ]Sin[θ]	z==Cos[θ]

⇔	(Subscript[β, 1])==Cos[ϕ]Sin[θ]	(Subscript[β, 2])==Sin[ϕ]Sin[θ]	α==(1/(Sqrt[2]))Cos[θ]

Within this coordinate system, the kinetic energy, electron-electron repulsion and external potential read

T==-(Sqrt[2]) t Sin[2 θ] (Cos[ϕ]+Sin[ϕ]);	(Subscript[V, ee])==U ((Sin[θ])^2);	V==-Δv Cos[2 ϕ] ((Sin[θ])^2);

{-4t α((Subscript[β, 1])+(Subscript[β, 2])),	U((Subscript[β, 1]^2)+(Subscript[β, 2]^2)),	Δn==2 ((Subscript[β, 1]^2)-(Subscript[β, 2]^2))}/.{(Subscript[β, 1])->Cos[ϕ]Sin[θ],(Subscript[β, 2])->Sin[ϕ]Sin[θ],α->(1/(Sqrt[2]))Cos[θ]}//FullSimplify
FullSimplify[{%[[1]]+%[[2]],%[[3]]}/.{U->2t u,Δn->2ρ}]//TrigToExp//FullSimplify

Solve[ρ==Cos[2 ϕ] ((Sin[θ])^2),ϕ]/.C[1]->0//FullSimplify
Solve[ρ==Cos[2 ϕ] ((Sin[θ])^2),θ]/.C[1]->0//FullSimplify

-2 t Sin[θ] (-u Sin[θ]+(Sqrt[2]) Cos[θ] (Cos[ϕ]+Sin[ϕ]))/.{ϕ->(1/2) ArcCos[ρ ((Csc[θ])^2)]}//FullSimplify
-2 t Sin[θ] (-u Sin[θ]+(Sqrt[2]) Cos[θ] (Cos[ϕ]+Sin[ϕ]))/.{θ->ArcSin[(Sqrt[ρ]) (Sqrt[Sec[2 ϕ]])]}//FullSimplify

Minimize[2 t (u+Cot[θ] ((Sqrt[1-ρ ((Csc[θ])^2)])-(Sqrt[1+ρ ((Csc[θ])^2)]))) ((Sin[θ])^2)/.{t->1,u->1,ρ->1/4},θ]//N
Minimize[-2 t Sin[θ] (Cos[θ] ((Sqrt[1-ρ ((Csc[θ])^2)])+(Sqrt[1+ρ ((Csc[θ])^2)]))-u Sin[θ])/.{t->1,u->1,ρ->1/4},θ]//N

min=Table[
res=Minimize[-2 t Sin[θ] (Cos[θ] ((Sqrt[1-ρ ((Csc[θ])^2)])+(Sqrt[1+ρ ((Csc[θ])^2)]))-u Sin[θ])/.{t->1,u->1},θ]//N;{Mod[θ/.res[[2]],π/2],res[[1]]//N},{ρ,0,1,1/20}]

Show[
Flatten[
{
Table[
With[{t=1,u=1},
Plot[{2 t (u+Cot[θ] ((Sqrt[1-ρ ((Csc[θ])^2)])-(Sqrt[1+ρ ((Csc[θ])^2)]))) ((Sin[θ])^2),-2 t Sin[θ] (Cos[θ] ((Sqrt[1-ρ ((Csc[θ])^2)])+(Sqrt[1+ρ ((Csc[θ])^2)]))-u Sin[θ])},{θ,0,π/2},ImageSize->Large,PlotRange->All]
]
,{ρ,0,1,1/20}]
,
ListPlot[min,PlotStyle->Black,PlotMarkers->{Automatic,Medium}]
}
]
]

f[θ_]=-2 t Sin[θ] (Cos[θ] ((Sqrt[1-ρ ((Csc[θ])^2)])+(Sqrt[1+ρ ((Csc[θ])^2)]))-u Sin[θ])/.{t->1,u->1}//FullSimplify
f'[θ]//FullSimplify
f''[θ]//FullSimplify

Other transformation

Solve[{Δn==2((Subscript[β, 1]^2)-(Subscript[β, 2]^2)),2(α^2)+(Subscript[β, 1]^2)+(Subscript[β, 2]^2)==1},{(Subscript[β, 1]),(Subscript[β, 2])}]//Last
{-4t α((Subscript[β, 1])+(Subscript[β, 2])),	U((Subscript[β, 1]^2)+(Subscript[β, 2]^2)),	2 ((Subscript[β, 1]^2)-(Subscript[β, 2]^2))}/.%//FullSimplify
FullSimplify[%/.{U->2t u,Δn->2ρ}]

Minimize[-2 t α ((Sqrt[2-4 (α^2)-2 ρ])+(Sqrt[2-4 (α^2)+2 ρ]))+2 t u (1-2 (α^2))/.{t->1/3,u->1,ρ->1/2},α]//N

Minimize[(Subscript[f, 1/3,1,1/2])[g],g]//N

Solve[{a==(Subscript[β, 1]^2)-(Subscript[β, 2]^2),b==2(Subscript[β, 1])(Subscript[β, 2])},{(Subscript[β, 1]),(Subscript[β, 2])}]//Last
{-4t α((Subscript[β, 1])+(Subscript[β, 2])),	U((Subscript[β, 1]^2)+(Subscript[β, 2]^2)),	2 ((Subscript[β, 1]^2)-(Subscript[β, 2]^2))}/.%//FullSimplify
FullSimplify[%/.{a->r Cos[θ],b->r Sin[θ]},r>0]

Solve[{a==(Subscript[β, 1]^2)-(Subscript[β, 2]^2),b==2(Subscript[β, 1])(Subscript[β, 2])},{(Subscript[β, 1]),(Subscript[β, 2])}]//Last
{-4t α((Subscript[β, 1])+(Subscript[β, 2]))/.α->(Sqrt[((1-(Subscript[β, 1]^2)-(Subscript[β, 2]^2))/2)]),	U((Subscript[β, 1]^2)+(Subscript[β, 2]^2)),	2 ((Subscript[β, 1]^2)-(Subscript[β, 2]^2))}/.%//FullSimplify
FullSimplify[%/.{a->r Cos[θ],b->r Sin[θ]},r>0]

FullSimplify[{-2 r (Sqrt[2-(r^2)]) t Cos[θ],(((r^2) U)/2)}/.{θ->(1/2) ArcSin[((Δn)/(r^2))]}/.{U->2t u,Δn->2ρ}]
FullSimplify[{2 (Sqrt[1-r]) t (Sqrt[r (1+Cos[θ])]) (-1+Cot[θ]-Csc[θ]),r U,2 r Cos[θ]}/.{θ->ArcSec[((2 r)/(Δn))]}/.{U->2t u,Δn->2ρ}]
FullSimplify[{2 (Sqrt[1-r]) t (Sqrt[r (1+Cos[θ])]) (-1+Cot[θ]-Csc[θ]),r U,2 r Cos[θ]}/.{r->(1/2) Δn Sec[θ]}/.{U->2t u,Δn->2ρ}]

Minimize[(Subscript[f, 1/3,1,1/2])[g],g]//N

min=Table[
res=Minimize[-2 (Sqrt[1-r]) t ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{t->1,u->1},r]//N;
{r/.res[[2]],res[[1]]//N},{ρ,0,1,1/50}]

Show[
Flatten[
{
With[{t=1,u=1},
Table[Plot[{2 (Sqrt[1-r]) t (-(Sqrt[r+ρ])+(Sqrt[r-ρ]))+2 r t u,-2 (Sqrt[1-r]) t ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u},{r,0,1},PlotRange->{-2,2},ImageSize->Large]
,{ρ,0,1,1/50}]]
,
ListPlot[min,PlotStyle->Black,PlotMarkers->{Automatic,Medium}],
Plot[{-(Sqrt[28-10 (Sqrt[5])+4 (Sqrt[60-26 (Sqrt[5])])]) (Sqrt[1-ρ])+(Sqrt[2]) (Sqrt[1+ρ]),-(Sqrt[28-10 (Sqrt[5])+4 (Sqrt[60-26 (Sqrt[5])])]) (Sqrt[1-((5 (7+5 (Sqrt[5])) (ρ^2))/(2 (5+(Sqrt[5])) (15+7 (Sqrt[5]))))])+(Sqrt[2]) (Sqrt[1+((5 (7+5 (Sqrt[5])) (ρ^2))/(2 (5+(Sqrt[5])) (15+7 (Sqrt[5]))))])},{ρ,(1/2)-(Sqrt[(1/20)]),1},PlotStyle->{Black,Red}]
}]
]

(1/((1/10) (5+(Sqrt[5]))))(((7+5 (Sqrt[5])) (ρ^2))/(4 (15+7 (Sqrt[5]))))

(((7+5 (Sqrt[5])) )/(4 (15+7 (Sqrt[5]))))//N

(1/10) (5+(Sqrt[5]))+(((7+5 (Sqrt[5])) (ρ^2))/(4 (15+7 (Sqrt[5]))))

f[r_]=-2 (Sqrt[1-r]) t ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u//Expand

(1-r)(Sqrt[r-ρ])(Sqrt[r+ρ])r f'[r]//Expand

D[-2 (Sqrt[1-r]) t ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{t->1,u->1},r]

Series[r/.Solve[Series[2-2 (Sqrt[1-r]) ((1/(2 (Sqrt[r-ρ])))+(1/(2 (Sqrt[r+ρ]))))+(((Sqrt[r-ρ])+(Sqrt[r+ρ]))/(Sqrt[1-r])),{ρ,0,2}]==0//Normal,r],{ρ,0,2}]//Last

Table[Root[9 (ρ^4)-12 (ρ^4) #1+(48 (ρ^2)+4 (ρ^4)) ((#1)^2)-128 (ρ^2) ((#1)^3)+(64+64 (ρ^2)) ((#1)^4)-320 ((#1)^5)+320 ((#1)^6)&,3],{ρ,0.,1.,0.1}]//N

Plot[]

Solve[(2-((2 (Sqrt[1-r]))/(Sqrt[r]))+((2 (Sqrt[r]))/(Sqrt[1-r])))+(((-3+2 r) (ρ^2))/(4 (Sqrt[1-r]) (r^(5/2))))==0,r]

Series[2 t u-2 (Sqrt[1-r]) t ((1/(2 (Sqrt[r-ρ])))+(1/(2 (Sqrt[r+ρ]))))+((t ((Sqrt[r-ρ])+(Sqrt[r+ρ])))/(Sqrt[1-r])),{}

Plot[(1/10) (5+(Sqrt[5]))+(((7+5 (Sqrt[5])) (ρ^2))/(4 (15+7 (Sqrt[5])))),{ρ,

f[r_]=-2 (Sqrt[1-r]) t ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{t->1,u->1}
f'[r]
f''[r]

Series[-2 (Sqrt[1-r]) t ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{u->1,t->1},{r,1-(Sqrt[5]),2}]//Normal

Minimize[(-1+(Sqrt[5])+r) (2-((5^(1/4))/(Sqrt[1-(Sqrt[5])-ρ]))+((Sqrt[1-(Sqrt[5])-ρ])/(5^(1/4)))-((5^(1/4))/(Sqrt[1-(Sqrt[5])+ρ]))+((Sqrt[1-(Sqrt[5])+ρ])/(5^(1/4))))-2 (-1+(Sqrt[5])+(5^(1/4)) (Sqrt[1-(Sqrt[5])-ρ])+(5^(1/4)) (Sqrt[1-(Sqrt[5])+ρ])),ρ]

D[-4 (Sqrt[1-r]) (Sqrt[r])+2 r+(((Sqrt[1-r]) (ρ^2))/(2 (r^(3/2)))),r]//FullSimplify

Solve[-3 (ρ^2)+2 r (4 r (-1+2 r+(Sqrt[-(-1+r) r]))+(ρ^2))==0,r]//RootReduce

Series[r/.Solve[ -3 (ρ^2)+2 r (4 r (-1+2 r+(Sqrt[-(-1+r) r]))+(ρ^2))==0,r]//Last,{ρ,0,2}]

Minimize[-2 t(Sqrt[1-r]) ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{t->1,u->1,ρ->0},r]//RootReduce
Minimize[-2 t(Sqrt[1-r]) ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{t->1,u->1,ρ->1},r]//RootReduce

Series[(Sqrt[1-r]) ((Sqrt[r+ρ])+(Sqrt[r-ρ])),{r,0,0}]

D[(1/10) (5+(Sqrt[5]))+(((7+5 (Sqrt[5])) (ρ^2))/(4 (15+7 (Sqrt[5])))),ρ]/.ρ->0

Minimize[-2 t (2 (Sqrt[1-r]) (Sqrt[r])-(((Sqrt[1-r]) (ρ^2))/(4 (r^(3/2)))))+ 2 r t u/.{t->1,u->1},r]
Minimize[-2 t(Sqrt[1-r]) ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{t->1,u->1,ρ->1},r]//RootReduce

Solve[((2-(1-(Sqrt[5])))/(1-(1/10) (5-(Sqrt[5]))))2+b==1]

Plot[{-(Sqrt[1-ρ]),-(Sqrt[1-ρ])},{ρ,0,1}]

f[ρ_]:=2 (1-(Sqrt[5]))+(1/2)((4-2(1-(Sqrt[5])))/(1-(1/10) (5-(Sqrt[5]))))ρ
{f[(1/10) (5-(Sqrt[5]))],f[1]}//Simplify

10ρ (5/(5-(Sqrt[5])))//FullSimplify

(1/10) (5+(Sqrt[5]))+(((7+5 (Sqrt[5])) (ρ^2))/(4 (15+7 (Sqrt[5]))))

g[ρ_]:=(Sqrt[a-ρ])+(Sqrt[b-ρ])
{g[(1/10) (5-(Sqrt[5]))],g[1]}//Simplify

g[ρ_]:=a(Sqrt[1-ρ])+b(Sqrt[1+ρ]);
res=Solve[{g[(1/10) (5-(Sqrt[5]))]==1-(Sqrt[5]),g[1]==2},{a,b}]//FullSimplify

g[ρ]/.{a->-(Sqrt[28-10 (Sqrt[5])+4 (Sqrt[60-26 (Sqrt[5])])]),b->(Sqrt[2])}

D[(1/10) (5+(Sqrt[5]))+(((7+5 (Sqrt[5])) (ρ^2))/(4 (15+7 (Sqrt[5])))),ρ]/.ρ->(1/10) (5-(Sqrt[5]))

Series[D[a(Sqrt[1-ρ])+b(Sqrt[1+ρ])+c,ρ],{ρ,0,2}]

Solve[-(a/4)-(b/4)==0]

g[ρ_]

g'[ρ]

-a(Sqrt[1-ρ])+b(Sqrt[1+ρ])/.{a->5 (Sqrt[(2/(5+(Sqrt[5])))])-(Sqrt[(10/(5+(Sqrt[5])))])+(Sqrt[((2 (15-(Sqrt[5])))/(5+(Sqrt[5])))]),b->(Sqrt[2])}//FullSimplify

c-(Sqrt[a+(Sqrt[ρ b])])/.{a->14-6 (Sqrt[5])-3 (Sqrt[2 (5-(Sqrt[5]))])+(Sqrt[10 (5-(Sqrt[5]))])-4 c+4 (Sqrt[5]) c+2 (Sqrt[2 (5-(Sqrt[5]))]) c+(c^2),b->8 (60-26 (Sqrt[5])-15 (Sqrt[2 (5-(Sqrt[5]))])+7 (Sqrt[10 (5-(Sqrt[5]))])-50 c+18 (Sqrt[5]) c+10 (Sqrt[2 (5-(Sqrt[5]))]) c-6 (Sqrt[10 (5-(Sqrt[5]))]) c+15 (c^2)-(Sqrt[5]) (c^2)+2 (Sqrt[10 (5-(Sqrt[5]))]) (c^2))}

c-(Sqrt[a+ρ b])/.{a->14-6 (Sqrt[5])-4 c+4 (Sqrt[5]) c+(c^2),b->-2 (5-3 (Sqrt[5])+2 (Sqrt[5]) c)}//FullSimplify

{c-(Sqrt[a+ρ b])==1-(Sqrt[5])/.ρ->(1/10) (5-(Sqrt[5])),c-(Sqrt[a+ρ b])==2/.ρ->1}

1-(Sqrt[((5(a+ρ b))/2)])+(Sqrt[(5+(Sqrt[5])) ρ])/.{a->(2/5) (11-5 (Sqrt[2])-2 (Sqrt[5])+5 (Sqrt[10])+3 (Sqrt[5+(Sqrt[5])])-(Sqrt[5 (5+(Sqrt[5]))])),b->(2/5) (-5+5 (Sqrt[2])+3 (Sqrt[5])-5 (Sqrt[10])-5 (Sqrt[5+(Sqrt[5])])+(Sqrt[5 (5+(Sqrt[5]))]))}//FullSimplify

Solve[{((((γ12)^2)(1-(Sqrt[1-((γ12)^2)])))/(2(((γ12)^2))))==2},{γ12}]

Solve[-(-5+(Sqrt[5])) a==5]

-2 (-1+(Sqrt[5]))//N

Solve[r==1-(Sqrt[(5-ρ)])+3(Sqrt[ρ]),ρ]

FullSimplify[SeriesCoefficient[(Sqrt[1-r]) ((Sqrt[r+ρ])+(Sqrt[r-ρ])),{r,0,n}],0<ρ<1]

(((-ρ))^(-n)) (Sqrt[ρ]) Binomial[(1/2),n] ((E^(I n π)) Hypergeometric2F1[-(1/2),-n,(3/2)-n,-ρ]+I Hypergeometric2F1[-(1/2),-n,(3/2)-n,ρ])

Solve[{((((γ12)^2)(1-(Sqrt[1-((γ12)^2)-(((γ11-1))^2)]))+2(((γ11-1))^2))/(2(((γ12)^2)+(((γ11-1))^2))))==1-(Sqrt[5])/.γ11->1}]

Limit[((((γ12)^2)(1-(Sqrt[1-((γ12)^2)-(((γ11-1))^2)]))+2(((γ11-1))^2))/(2(((γ12)^2)+(((γ11-1))^2)))),γ11->0]
Limit[((((γ12)^2)(1-(Sqrt[1-((γ12)^2)-(((γ11-1))^2)]))+2(((γ11-1))^2))/(2(((γ12)^2)+(((γ11-1))^2)))),γ11->1]

((2+((γ12)^2)+(((-((γ12)^2)))^(3/2)))/(2 (1+((γ12)^2))))

Inverse[

Plot[(1/10) (4+ρ-2 (Sqrt[4+2 ρ-(ρ^2)])),{ρ,0,1}]

Minimize[-2 t(Sqrt[1-r]) ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+ 2 r t u/.{t->1/3,u->1,ρ->1/5},r]//N
Minimize[2 t (-1+Cot[θ]-Csc[θ]) (Sqrt[ρ (1+Sec[θ])]) (Sqrt[1-ρ Sec[θ]])+2 t u ρ Sec[θ]/.{t->1/3,u->1,ρ->1/5},θ]//N

P[u_,ρ_,g_]:=(4+(u^2))(g^6)/4+((ρ^2)(3+(u^2))-1)(g^4)+2u (ρ^2)(g^3)+(ρ^2)((ρ^2)(3+(u^2))-(2+(u^2)))(g^2)-2u (ρ^2)(1-(ρ^2))g-(ρ^4)(1-(ρ^2));

r[u_,ρ_]:=g/.Solve[P[u,ρ,g]==0,g]

r[u,ρ][[4]]/.{u->1,ρ->1/2}//N

NSolve[P[1,1,g]==0]

Solve[P[u,ρ,g]==0/.{u->1,ρ->1}]//N

Table[r[[k]],{k,Length[r]}]

r[u,ρ][[1]]

Table[Evaluate[Minimize[-2 t(Sqrt[1-r])  ((Sqrt[r+ρ])+(Sqrt[r-ρ]))+2 r t u/.t->1,r]],{u,1,2,1},{ρ,1,2,1}]

With[{t=1},Table[{u,ρ,Minimize[2 (q^2) t u-2 q (Sqrt[1-(q^2)]) t ((Sqrt[1-(ρ/(q^2))])+(Sqrt[1+(ρ/(q^2))])),q][[1]]},{u,1,2,1},{ρ,1,2,1}]]

ListPlot3D[,ImageSize->800]

Plot3D[Evaluate[Table[r[[k]],{k,Length[r]}]],{u,0,10},{ρ,-5,5},ImageSize->500,AspectRatio->{1,1}]

Minimize[(Subscript[f, 1,1,1/2])[g],g]//N

NSolve[(Subscript[P, 1,1/2])[g]==0]
NSolve[D[(Subscript[P, 1,1/2])[g],g]==0]
NSolve[D[(Subscript[P, 1,1/2])[g],g,g]==0]
NSolve[D[(Subscript[P, 1,1/2])[g],g,g,g]==0]

CoefficientList[(Subscript[P, u,ρ])[g]-((g^2)+(x^2))(a (g^4)+b (g^3)+c (g^2)+d g+e)//Expand,g]

Solve[CoefficientList[(Subscript[P, u,ρ])[g]-((Subscript[a, 1]) (g^3)+(Subscript[a, 2]) (g^2)+(Subscript[a, 3]) g +(Subscript[a, 3]))((Subscript[b, 1]) (g^3)+(Subscript[b, 2]) (g^2)+ (Subscript[b, 3])g+ (Subscript[b, 4]))//Expand,g]==0,{(Subscript[b, 1]),a,b,c,d,e}]

FullSimplify[(Subscript[P, u,ρ])[g]==(1+((((u/2)))^2))(g^6)+((ρ^2)(3+(u^2))-1)(g^4)+2u (ρ^2)(g^3)+(ρ^4)(3+(u^2))(g^2)-(ρ^2)(2+(u^2))(g^2)-2u (ρ^2)(1-(ρ^2))g-(ρ^4)(1-(ρ^2))]

Solve[{(r^2)==(g^2)+(ρ^2),R==g ρ},{ρ,g}]

(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2))))/.{ρ->r Cos[θ],g->r Sin[θ]}//FullSimplify

(Subscript[h, ρ_])[g_]:=(((g^2)(1-(Sqrt[1-(g^2)-(ρ^2)]))+2(ρ^2))/(2((g^2)+(ρ^2))));(Subscript[H, θ_])[r_]:=(1/4) (3-(Sqrt[1-(r^2)])+(1+(Sqrt[1-(r^2)])) Cos[2 θ]);

Solve[D[(Subscript[H, θ])[r],r]-(1/u)==0/.{u->1,ρ->1/2},g]//N

Solve[D[(Subscript[h, ρ])[g],g]-(1/u)==0/.{u->1,ρ->1/2},g]//N

Solve[D[(Subscript[h, ρ])[g],g]-(1/u)==0/.{u->1,ρ->1/2},g]//N

Collect[(Subscript[P, u,ρ])[g],u]
Collect[(Subscript[P, u,ρ])[g],ρ]

𝒫𝒯-symmetric Hubbard dimer

(({{2(Subscript[v, 1])+U, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, 2 (Subscript[v, 2])+U}}))==((Subscript[v, 1])+(Subscript[v, 2]))(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))+(({{-Δv+U, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, Δv+U}}))/.Δv->(Subscript[v, 2])-(Subscript[v, 1])/.(Subscript[v, 2])->-(Subscript[v, 1])//Simplify

|1↑1↓⟩, |1↑2↓⟩,|1↓2↑⟩,|2↑2↓⟩
			↓ 𝒫 = site flip
|2↑2↓⟩, |2↑1↓⟩,|2↓1↑⟩,|1↑1↓⟩
			↓ 𝒯 = spin flip
|2↓2↑⟩, |2↓1↑⟩,|2↑1↓⟩,|1↓1↑⟩

(Subscript[ϵ, 2])[t_,U_,Δv_]=Eigenvalues[(({{-I Δv+ U, -I t, I t, 0}, {-I t, 0, 0, -I t}, {I t, 0, 0, I t}, {0, -I t, I t, I Δv+U}}))]//FullSimplify

Show[{
With[{t=(1/2),U=1},Plot[Evaluate[(Subscript[ϵ, 2])[t,U,Δv]],{Δv,0,3},AspectRatio->1,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E(N=2)"},PlotLegends->Placed[Automatic,{Right,Center}]]]
,
With[{t=(1/2),U=1},Plot[Evaluate[Im[(Subscript[ϵ, 2])[t,U,Δv]]],{Δv,0,3},AspectRatio->1,PlotStyle->Dashed,PlotTheme->"Web",BaseStyle->16,FrameLabel->{"Δv","E(N=2)"}]]
}]

Asymmetric Hubbard dimer in the complex plane

(Subscript[ϵ, 2])[t_,U_,Δv_]=Eigenvalues[(({{-Δv+U, -t, t, 0}, {-t, 0, 0, -t}, {t, 0, 0, t}, {0, -t, t, Δv+U}}))]//FullSimplify

With[{t=1/2,Δv=1},
Show[{
Graphics3D[{Specularity[White,30],Gray,Opacity[0.2],Sphere[{0,0,0},1]}],
Plot3D[Re[(Subscript[ϵ, 2])[t,Ux+I Uy,Δv]][[2]],{Ux,-3,3},{Uy,-3,3},Mesh->None,PlotStyle->Directive[Red,Opacity[0.2]],PlotTheme->"Web",AspectRatio->1,BaseStyle->16,PlotLegends->Placed[Automatic,{Right,Center}]],
Plot3D[Re[(Subscript[ϵ, 2])[t,Ux+I Uy,Δv]][[3]],{Ux,-3,3},{Uy,-3,3},Mesh->None,PlotStyle->Directive[Blue,Opacity[0.2]],PlotTheme->"Web",AspectRatio->1,BaseStyle->16,PlotLegends->Placed[Automatic,{Right,Center}]],
Plot3D[Re[(Subscript[ϵ, 2])[t,Ux+I Uy,Δv]][[4]],{Ux,-3,3},{Uy,-3,3},Mesh->None,PlotStyle->Directive[Green,Opacity[0.2]],PlotTheme->"Web",AspectRatio->1,BaseStyle->16,PlotLegends->Placed[Automatic,{Right,Center}]]

}]
]

With[{t=1/2,Δv=0},
Plot[Evaluate[Re[(Subscript[ϵ, 2])[t,U,Δv][[2;;4]]]],{U,-3,3},Mesh->None,PlotStyle->Opacity[0.2],PlotTheme->"Web",AspectRatio->1,BaseStyle->16,PlotLegends->Placed[Automatic,{Right,Center}]]
]