Initialisation

(Subscript[Y, ℓ_,m_])[θ_,ϕ_]=SphericalHarmonicY[ℓ,m,θ,ϕ];

(Subscript[Y, ℓ_])[θ_]=SphericalHarmonicY[ℓ,0,θ,ϕ];

Needs["MaTeX`"]
SetOptions[MaTeX,"Preamble"->{"\\usepackage{amssymb,amsmath,latexsym,amsfonts,amsthm,mathpazo,xcolor,bm,mhchem}","\\definecolor{darkgreen}{RGB}{0, 180, 0}","\\definecolor{darkcyan}{RGB}{0, 180, 180}","\\definecolor{darkmagenta}{RGB}{180,0, 180}"}];

Useful calculs and results

Hyperspherical rotation matrix

mapping=CoordinateTransformData[{"Hyperspherical"->"Cartesian",3},"Mapping"]

MatrixForm[rotation=CoordinateTransformData[{"Hyperspherical"-> "Cartesian",4},"OrthonormalBasisRotation",{r,θ,ϕ,α}]]

MatrixForm[rotation=CoordinateTransformData["Spherical"-> "Cartesian","OrthonormalBasisRotation",{r,θ,ϕ}]]

MatrixForm[rotation=CoordinateTransformData["Polar"-> "Cartesian","OrthonormalBasisRotation",{r,θ}]]

Two electron integrals

(ij|kl)=((((Subscript[Y, i,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])/R)(((Subscript[Y, j,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])/R)|(1/(Subscript[r, 12]))|(((Subscript[Y, k,0])[π-(Subscript[θ, 2]),(Subscript[ϕ, 2])])/R)(((Subscript[Y, l,0])[π-(Subscript[θ, 2]),(Subscript[ϕ, 2])])/R))
                 =((((-1))^(k+l))/R)(∫_0^π)(∫_0^(2π))(∫_0^π)(∫_0^(2π))(Subscript[Y, i])[(Subscript[θ, 1])](Subscript[Y, j])[(Subscript[θ, 1])]((∑_(L=0)^∞)((4π)/(2L+1))(Subscript[Y, L])[(Subscript[θ, 1])](Subscript[Y, L])[(Subscript[θ, 2])])(Subscript[Y, k])[(Subscript[θ, 2])](Subscript[Y, l])[(Subscript[θ, 2])]Sin[(Subscript[θ, 1])]Sin[(Subscript[θ, 2])]\[DifferentialD](Subscript[ϕ, 2])\[DifferentialD](Subscript[θ, 2])\[DifferentialD](Subscript[ϕ, 1])\[DifferentialD](Subscript[θ, 1])
	        =((((-1))^(k+l))/R)(∑_(L=0)^∞)((4π)/(2L+1))((∫_0^π)(∫_0^(2π))(Subscript[Y, i])[θ](Subscript[Y, j])[θ](Subscript[Y, L])[θ]Sin[θ]\[DifferentialD]ϕ\[DifferentialD]θ)((∫_0^π)(∫_0^(2π))(Subscript[Y, k])[θ](Subscript[Y, l])[θ](Subscript[Y, L])[θ]Sin[θ]\[DifferentialD]ϕ\[DifferentialD]θ)
	        =((((-1))^(k+l))/R)(Sqrt[(2i+1)(2j+1)(2k+1)(2l+1)])(∑_(L=0)^∞)(((({{i, j, L}, {0, 0, 0}})))^2)(((({{k, l, L}, {0, 0, 0}})))^2)

(im,jn|ko,lp)=((((Subscript[Y, i,m])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])/R)(((Subscript[Y, j,n])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])/R)|(1/(Subscript[r, 12]))|(((Subscript[Y, k,o])[π-(Subscript[θ, 2]),(Subscript[ϕ, 2])])/R)(((Subscript[Y, l,p])[π-(Subscript[θ, 2]),(Subscript[ϕ, 2])])/R))
                 =((((-1))^(k+l))/R)(∑_(L=0)^∞)((4π)/(2L+1))(∑_(M=-L)^L)((∫_0^π)(∫_0^(2π))(Subscript[Y, i,m])[θ](Subscript[Y, j,n])[θ](Subscript[Y, L,M])[θ]Sin[θ]\[DifferentialD]ϕ\[DifferentialD]θ)((∫_0^π)(∫_0^(2π))(Subscript[Y, k,o])[θ](Subscript[Y, l,p])[θ](Subscript[Y, L,M])[θ]Sin[θ]\[DifferentialD]ϕ\[DifferentialD]θ)
	        =((((-1))^(k+l))/R)(Sqrt[(2i+1)(2j+1)(2k+1)(2l+1)])(∑_(L=0)^(Max[i+j,k+l]))(∑_(M=-L)^L)(((({{i, j, L}, {m, n, M}})))^2)(((({{k, l, L}, {o, p, M}})))^2)

TEIu[{i_,m_},{j_,n_},{k_,o_},{l_,p_}]:=((((-1))^(k+l))/R)(Sqrt[(2i+1)(2j+1)(2k+1)(2l+1)])(∑_(L=Max[Abs[i-j],Abs[k-l]])^(Min[i+j,k+l]))(∑_(M=-L)^L)((ThreeJSymbol[{i,m},{j,n},{L,M}])^2)((ThreeJSymbol[{k,o},{l,p},{L,M}])^2)

Without the spherium shortcuts (π - (Subscript[θ, 2])) there are no (((-1))^(k+l)) factor

TEIr[{i_,m_},{j_,n_},{k_,o_},{l_,p_},R_]:=(1/R)(Sqrt[(2i+1)(2j+1)(2k+1)(2l+1)])(∑_(L=Max[Abs[i-j],Abs[k-l]])^(Min[i+j,k+l]))(∑_(M=-L)^L)((ThreeJSymbol[{i,m},{j,n},{L,M}])^2)((ThreeJSymbol[{k,o},{l,p},{L,M}])^2)

Minimal basis verification

{(Subscript[ψ, 1])[θ_]=(Subscript[Y, 0])[θ],(Subscript[ψ, 2])[θ_]=(Subscript[Y, 1])[θ]}

Table[((((2π))^2)/R)(∫_0^π)(∫_0^π)(Subscript[ψ, μ])[θ1](Subscript[ψ, ν])[θ1](1+Cos[θ1] Cos[θ2]+(1/4) (-1+3 ((Cos[θ1])^2)) (-1+3 ((Cos[θ2])^2))) (Subscript[ψ, λ])[θ2](Subscript[ψ, σ])[θ2]Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2,{μ,2},{ν,2},{λ,2},{σ,2}]

Table[TEIr[{μ,0},{ν,0},{λ,0},{σ,0}],{μ,0,1},{ν,0,1},{λ,0,1},{σ,0,1}]

Table[((((2π))^2)/R)(∫_0^π)(∫_0^π)(Subscript[ψ, μ])[θ1](Subscript[ψ, ν])[θ1](1+Cos[θ1] Cos[θ2]+(1/4) (-1+3 ((Cos[θ1])^2)) (-1+3 ((Cos[θ2])^2))) (Subscript[ψ, λ])[π-θ2](Subscript[ψ, σ])[π-θ2]Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2,{μ,2},{ν,2},{λ,2},{σ,2}]

Table[TEIu[{μ,0},{ν,0},{λ,0},{σ,0}],{μ,0,1},{ν,0,1},{λ,0,1},{σ,0,1}]

UHF equations and SCF algorithm

(N^α),(N^β) number of electrons with spin up or down

Introduction of a basis set:

(Subscript[ψ, i]^α)=(∑_(μ=1)^K)(Subscript[c, μ i]^α)(Subscript[ϕ, μ])       i=1,...,K
(Subscript[ψ, i]^β)=(∑_(μ=1)^K)(Subscript[c, μ i]^β)(Subscript[ϕ, μ])       i=1,...,K

We define the unrestricted density matrix as

(Subscript[P, μ ν]^α)=(∑_a^(N^α))(Subscript[c, μ a]^α)(SuperStar[((Subscript[c, ν a]^α))])
(Subscript[P, μ ν]^β)=(∑_a^(N^β))(Subscript[c, μ a]^β)(SuperStar[((Subscript[c, ν a]^β))])

We use the following notation for 2 electrons integrals

(μν|σλ)=((Subscript[ϕ, μ])(Subscript[ϕ, ν])|(Subscript[ϕ, σ])(Subscript[ϕ, λ]))
                  = ∫∫(Subscript[ϕ, μ])(SuperStar[(1)])(Subscript[ϕ, ν])(1)(1/(Subscript[r, 12]))(Subscript[ϕ, σ])(SuperStar[(2)])(Subscript[ϕ, λ])(2)\[DifferentialD](Subscript[Ω, 1])\[DifferentialD](Subscript[Ω, 2])

The expressions of the unrestricted Fock operator are

(Subscript[F, μ ν]^α)=(Subscript[H, μ ν]^core)+(∑_λ)(∑_σ)(Subscript[P, λ σ]^T)(μν|σλ)-(Subscript[P, λ σ]^α)(μλ|σν)
(Subscript[F, μ ν]^β)=(Subscript[H, μ ν]^core)+(∑_λ)(∑_σ)(Subscript[P, λ σ]^T)(μν|σλ)-(Subscript[P, λ σ]^β)(μλ|σν)

If we choose a basis set of size K, the basis set will be of size (K^2)

SCF algorithm in an orthonormal basis:

1) Calculate all required molecular integrals (Subscript[S, μν]) , (Subscript[((H^core)), μν]) and (μν | λσ) 
2) Obtain a guess at the density matrix P
3) Calculate G from P and the two electrons integrals, add (H^core)  to get the Fock matrix F
4)  Solve (C^t)F C=ϵ , obtain the new C’
5) Calculate the new density matrix P’ and determine whether the procedure has converged

UHF minimal basis minimization

{ψ1α[θ_]=Cos[χα](Subscript[Y, 0])[θ]+Sin[χα](Subscript[Y, 1])[θ],
ψ1β[θ_]=Cos[χβ](Subscript[Y, 0])[θ]+Sin[χβ](Subscript[Y, 1])[θ],
ψ2α[θ_]=-Sin[χα](Subscript[Y, 0])[θ]+Cos[χα](Subscript[Y, 1])[θ],
ψ2β[θ_]=-Sin[χβ](Subscript[Y, 0])[θ]+Cos[χβ](Subscript[Y, 1])[θ]};

Ground state energy:

EUHF[R_,λ_]=2(∑_(L=0)^1)(Subscript[c, L]^2)((L(L+1))/(2(R^2)))+(∑_(i=0)^1)(∑_(j=0)^1)(∑_(k=0)^1)(∑_(l=0)^1)(Subscript[c, i])(Subscript[c, j])(Subscript[c, k])(Subscript[c, l])(((-1))^(k+l))(λ/R)(Sqrt[(2i+1)(2j+1)(2k+1)(2l+1)])(∑_(L=0)^2)((ThreeJSymbol[{i,0},{j,0},{L,0}])^2)((ThreeJSymbol[{k,0},{l,0},{L,0}])^2)/.{(Subscript[c, 0])->Cos[χ],(Subscript[c, 1])->Sin[χ]}

Solve[(∂_χ)EUHF[R,λ]==0,χ]

We choose:

χα->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]

ArcTan

And the symmetry gives us

χβ->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]

(∂_(χ,χ))EUHF[R,1]/.χ->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]//FullSimplify

Plot[(((-3+2 R) (75+62 R))/(42 (R^3))),{R,-5,5}]

Log

ERHF[R_,λ_]=2(∑_(L=0)^1)(Subscript[c, L]^2)((L(L+1))/(2(R^2)))+(∑_(i=0)^1)(∑_(j=0)^1)(∑_(k=0)^1)(∑_(l=0)^1)(Subscript[c, i])(Subscript[c, j])(Subscript[c, k])(Subscript[c, l])(λ/R)(Sqrt[(2i+1)(2j+1)(2k+1)(2l+1)])(∑_(L=0)^2)((ThreeJSymbol[{i,0},{j,0},{L,0}])^2)((ThreeJSymbol[{k,0},{l,0},{L,0}])^2)/.{(Subscript[c, 0])->Cos[χ],(Subscript[c, 1])->Sin[χ]}

Solve[(∂_χ)ERHF[R,λ]==0,χ]

RHF vs UHF: minimal basis, Fock matrix, Hamiltonian and spin contamination

One electron basis and Fock matrix

RHF

{(Subscript[ψ, 1])[θ_]=(Subscript[Y, 0])[θ],(Subscript[ψ, 2])[θ_]=(Subscript[Y, 1])[θ]}

c={{1,0},{0,1}};o={{1,0}}\[Transpose];P=2o.o\[Transpose];

Hc=Table[(((2π))/(2(R^2)))(∫_0^π)Expand[((∂_θ)(Subscript[ψ, μ])[θ])((∂_θ)(Subscript[ψ, ν])[θ])]Sin[θ]\[DifferentialD]θ,{μ,2},{ν,2}]

Hc[R_]:=(({{0, 0}, {0, (1/(R^2))}}))

G=Table[TEIr[{μ,0},{ν,0},{λ,0},{σ,0}],{μ,0,1},{ν,0,1},{λ,0,1},{σ,0,1}];

f= Hc+Table[(∑_(λ=1)^2)(∑_(σ=1)^2)P[[λ,σ]](G[[μ,ν,λ,σ]]-(1/2)G[[μ,σ,λ,ν]]),{μ,2},{ν,2}]

F=c\[Transpose].f.c

F[R_]:=(({{(1/R), 0}, {0, (1/(R^2))+(5/(3 R))}}));

UHF

{(Subscript[ψα, 1])[θ_]=c10α (Subscript[Y, 0])[θ]+c11α (Subscript[Y, 1])[θ],
(Subscript[ψβ, 1])[θ_]=c10β (Subscript[Y, 0])[θ]+c11β (Subscript[Y, 1])[θ],
(Subscript[ψα, 2])[θ_]=c20α (Subscript[Y, 0])[θ]+c21α (Subscript[Y, 1])[θ],
(Subscript[ψβ, 2])[θ_]=c20β (Subscript[Y, 0])[θ]+c21β (Subscript[Y, 1])[θ]};

The coeffcients are obtained in the useful calculs and results part

{(Subscript[ψα, 1])[θ_,R_]:=Cos[χα](Subscript[Y, 0])[θ]+Sin[χα](Subscript[Y, 1])[θ]/.χα->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],
(Subscript[ψβ, 1])[θ_,R_]:=Cos[χβ](Subscript[Y, 0])[θ]+Sin[χβ](Subscript[Y, 1])[θ]/.χβ->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],
(Subscript[ψα, 2])[θ_,R_]:=-Sin[χα](Subscript[Y, 0])[θ]+Cos[χα](Subscript[Y, 1])[θ]/.χα->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],
(Subscript[ψβ, 2])[θ_,R_]:=-Sin[χβ](Subscript[Y, 0])[θ]+Cos[χβ](Subscript[Y, 1])[θ]/.χβ->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]}//FullSimplify;

-Sin[χα]y0+Cos[χα]y1/.χα->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]//FullSimplify

(Subscript[ψα, 1])[θ,R]//FullSimplify

(Subscript[ψβ, 1])[θ,R]//FullSimplify

(Subscript[ψα, 2])[θ,R]//FullSimplify

cα={{Cos[χα],Sin[χα]},{-Sin[χα],Cos[χα]}}\[Transpose];	oα={{Cos[χα],Sin[χα]}}\[Transpose];	Pα=oα.oα\[Transpose];
cβ={{Cos[χβ],Sin[χβ]},{-Sin[χβ],Cos[χβ]}}\[Transpose];	oβ={{Cos[χβ],Sin[χβ]}}\[Transpose];	Pβ=oβ.oβ\[Transpose];

G=Table[TEIr[{μ,0},{ν,0},{λ,0},{σ,0},R],{μ,0,1},{ν,0,1},{λ,0,1},{σ,0,1}];

fα=Hc[R]+Table[(∑_(λ=1)^2)(∑_(σ=1)^2)Pα[[λ,σ]](G[[μ,ν,λ,σ]]-G[[μ,σ,λ,ν]]),{μ,2},{ν,2}]+Table[(∑_(λ=1)^2)(∑_(σ=1)^2)Pβ[[λ,σ]](G[[μ,ν,λ,σ]]),{μ,2},{ν,2}];
fβ=Hc[R]+Table[(∑_(λ=1)^2)(∑_(σ=1)^2)Pβ[[λ,σ]](G[[μ,ν,λ,σ]]-G[[μ,σ,λ,ν]]),{μ,2},{ν,2}]+Table[(∑_(λ=1)^2)(∑_(σ=1)^2)Pα [[λ,σ]](G[[μ,ν,λ,σ]]),{μ,2},{ν,2}];

fα/.{χα->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],χβ->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]}//FullSimplify

{{((-75+218 R)/(168 (R^2))),-((5 (Sqrt[-3+2 R]) (Sqrt[75+62 R]))/(84 (R^2)))},{-((5 (Sqrt[-3+2 R]) (Sqrt[75+62 R]))/(84 (R^2))),((225+242 R)/(168 (R^2)))}}//Expand

{Fα,Fβ}={cα\[Transpose].fα.cα,cβ\[Transpose].fβ.cβ}/.{χα->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],χβ->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]}//FullSimplify

{Fα,Fβ}={(({{((75+118 R)/(168 (R^2))), 0}, {0, ((25+114 R)/(56 (R^2)))}})),(({{((75+118 R)/(168 (R^2))), 0}, {0, ((25+114 R)/(56 (R^2)))}}))}//Expand

Plot[{((75+118 R)/(168 (R^2))),((25+114 R)/(56 (R^2)))},{R,-5,5}]

NRJ of the Fock orbitals

Plot[{(1/R),(1/(R^2))+(5/(3R)),((75+118 R)/(168 (R^2))),((25+114 R)/(56 (R^2)))},{R,0,5},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->{MaTeX["\\text{GS RHF}"],MaTeX["\\text{ES RHF}"],MaTeX["\\text{GS UHF}"],MaTeX["\\text{ES UHF}"]},FrameLabel->{MaTeX["R",Magnification->1.5],MaTeX["E",Magnification->1.5]},AspectRatio->1]

We observe the Coulson-Fischer point at R=(3/2), the symmetry is broken for R>(3/2).

Two electron basis, Hamiltonian and spin contamination

RHF

{Ψ1=(Subscript[ψ, 1])[θ1](Subscript[ψ, 1])[θ2],Ψ2=(Subscript[ψ, 1])[θ1](Subscript[ψ, 2])[θ2],Ψ3=(Subscript[ψ, 2])[θ1](Subscript[ψ, 1])[θ2],Ψ4=(Subscript[ψ, 2])[θ1](Subscript[ψ, 2])[θ2]};

S[Ψ1_,Ψ2_]:=(((2π))^2)(∫_0^π)(∫_0^π) Expand[Ψ1 Ψ2] Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2//FullSimplify
H[Ψ1_,Ψ2_]:=((((2π))^2)/(2(R^2)))(∫_0^π)(∫_0^π)Expand[((∂_(θ1))Ψ1)((∂_(θ1))Ψ2)+((∂_(θ2))Ψ1)((∂_(θ2))Ψ2)]Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2+((((2π))^2)/R)(∫_0^π)(∫_0^π) Expand[Ψ1 Ψ2(1+Cos[θ1] Cos[θ2]+(1/4) (-1+3 ((Cos[θ1])^2)) (-1+3 ((Cos[θ2])^2)))] Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2//FullSimplify

{rH,rS}={(({{H[Ψ1,Ψ1], H[Ψ2,Ψ1], H[Ψ3,Ψ1], H[Ψ4,Ψ1]}, {H[Ψ1,Ψ2], H[Ψ2,Ψ2], H[Ψ3,Ψ2], H[Ψ4,Ψ2]}, {H[Ψ1,Ψ3], H[Ψ2,Ψ3], H[Ψ3,Ψ3], H[Ψ4,Ψ3]}, {H[Ψ1,Ψ4], H[Ψ2,Ψ4], H[Ψ3,Ψ4], H[Ψ4,Ψ4]}})),(({{S[Ψ1,Ψ1], S[Ψ2,Ψ1], S[Ψ3,Ψ1], S[Ψ4,Ψ1]}, {S[Ψ1,Ψ2], S[Ψ2,Ψ2], S[Ψ3,Ψ2], S[Ψ4,Ψ2]}, {S[Ψ1,Ψ3], S[Ψ2,Ψ3], S[Ψ3,Ψ3], S[Ψ4,Ψ3]}, {S[Ψ1,Ψ4], S[Ψ2,Ψ4], S[Ψ3,Ψ4], S[Ψ4,Ψ4]}}))}//FullSimplify;

{rH[R_],rS}:={(({{(1/R), 0, 0, (1/(3 R))}, {0, ((1+R)/(R^2)), (1/(3 R)), 0}, {0, (1/(3 R)), ((1+R)/(R^2)), 0}, {(1/(3 R)), 0, 0, ((50+29 R)/(25 (R^2)))}})),(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))};

The two electron basis is orthonormal. The Hamiltonian is block diagonal: (s^2) interact with ((Subscript[p, z])^2) and (Subscript[sp, z]) interact with (Subscript[p, z])s.
The singlet symmetry block has the same expression as in the minimal basis of legendre polynomials: those polynomials are singlet CSF.

UHF

{Ψ1=(Subscript[ψα, 1])[θ1](Subscript[ψβ, 1])[θ2],Ψ2=(Subscript[ψα, 1])[θ1](Subscript[ψβ, 2])[θ2],Ψ3=(Subscript[ψα, 2])[θ1](Subscript[ψβ, 1])[θ2],Ψ4=(Subscript[ψα, 2])[θ1](Subscript[ψβ, 2])[θ2]}//FullSimplify;

{uH,uS}={(({{H[Ψ1,Ψ1], H[Ψ2,Ψ1], H[Ψ3,Ψ1], H[Ψ4,Ψ1]}, {H[Ψ1,Ψ2], H[Ψ2,Ψ2], H[Ψ3,Ψ2], H[Ψ4,Ψ2]}, {H[Ψ1,Ψ3], H[Ψ2,Ψ3], H[Ψ3,Ψ3], H[Ψ4,Ψ3]}, {H[Ψ1,Ψ4], H[Ψ2,Ψ4], H[Ψ3,Ψ4], H[Ψ4,Ψ4]}})),(({{S[Ψ1,Ψ1], S[Ψ2,Ψ1], S[Ψ3,Ψ1], S[Ψ4,Ψ1]}, {S[Ψ1,Ψ2], S[Ψ2,Ψ2], S[Ψ3,Ψ2], S[Ψ4,Ψ2]}, {S[Ψ1,Ψ3], S[Ψ2,Ψ3], S[Ψ3,Ψ3], S[Ψ4,Ψ3]}, {S[Ψ1,Ψ4], S[Ψ2,Ψ4], S[Ψ3,Ψ4], S[Ψ4,Ψ4]}}))}//FullSimplify;

{uH[R_],uS}:={(({{((-225+4 R (75+59 R))/(336 (R^3))), 0, 0, ((-75+4 (-3+R) R)/(112 (R^3)))}, {0, ((5 (-45+60 R+92 (R^2)))/(336 (R^3))), ((-75+4 (-3+R) R)/(112 (R^3))), (((Sqrt[-3+2 R]) (25+2 R) (Sqrt[75+62 R]))/(280 (R^3)))}, {0, ((-75+4 (-3+R) R)/(112 (R^3))), ((5 (-45+60 R+92 (R^2)))/(336 (R^3))), (((Sqrt[-3+2 R]) (25+2 R) (Sqrt[75+62 R]))/(280 (R^3)))}, {((-75+4 (-3+R) R)/(112 (R^3))), (((Sqrt[-3+2 R]) (25+2 R) (Sqrt[75+62 R]))/(280 (R^3))), (((Sqrt[-3+2 R]) (25+2 R) (Sqrt[75+62 R]))/(280 (R^3))), ((16875+4 R (2775+1511 R))/(8400 (R^3)))}})),(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))};

The two electron basis is orthonormal. We see that the singlet state ((Subscript[p, z])^2) is spin contaminated for R>(3/2) i.e. the Coulson-Fischer point by the triplet states. Coherent with the UHF description of (Subscript[H, 2]) Gill 1988.
See equation (3.6) of Gill 1988 for the expression of <(S^2) >.

Plot[(((Sqrt[-3+2 R]) (25+2 R) (Sqrt[75+62 R]))/(280 (R^3))),{R,0,10}]

(((-75+4 (-3+R) R)/(112 (R^3)))/(((-225+4 R (75+59 R))/(336 (R^3)))-((16875+4 R (2775+1511 R))/(8400 (R^3)))))//FullSimplify

Plot[{Abs[((25 (75-4 (-3+R) R))/(12 (((25+2 R))^2)))],0.25},{R,0,20}]

Electron density

(Subscript[ψα, 1])[θ]

((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))+((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))//FullSimplify

(((Sqrt[75+62 R])+5 (Sqrt[-9+6 R]) Cos[θ])/(8 (Sqrt[7 π]) (Sqrt[R])))/.R->10//N//Simplify

(Subscript[ψβ, 1])[θ]

With[{R=10},Plot3D[(((((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))+((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))))^2),{θ,0,π},{ϕ,0,2π}]]

With[{R=2},Plot3D[(((((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))-((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))))^2),{θ,0,π},{ϕ,0,2π}]]

With[{R=1},Plot3D[{Re[(((((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))-((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ1])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))))^2)],Re[(((((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))+((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ2])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))))^2)],Re[(((((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))-((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ1])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))))^2)(((((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))+((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ2])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))))^2)]},{θ1,0,π},{θ2,0,π},Mesh->None,PlotStyle->Opacity[0.7]]]

With[{R=1},Plot3D[{((Abs[(((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))-((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ1])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)]))))])^2),((Abs[(((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))+((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ2])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)]))))])^2),((Abs[(((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))-((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ1])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)]))))(((Sqrt[75+62 R])/(2 (Sqrt[π]) (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)])))+((5 (Sqrt[(3/π)]) (Sqrt[-3+2 R]) Cos[θ2])/(2 (Sqrt[R]) (Sqrt[((25 (-3+2 R))/R)+((75+62 R)/R)]))))])^2)},{θ1,0,π},{θ2,0,π},Mesh->None,PlotStyle->Opacity[0.7]]]

Exact energies

(Subscript[A, i_,j_]):=(((i+1))^(3/2))(((j+1))^(3/2))If[i>j,(j/(2 (((1+i))^2)))+((j^2)/(4 (((1+i))^2))),If[i<j,(i/(2 (((1+j))^2)))+((i^2)/(4 (((1+j))^2))),((i (((2+i))^2))/(4 (((1+i))^3)))]];
(Subscript[B, i_,j_]):=(((i+1))^(3/2))(((j+1))^(3/2))If[i>j,(1/(((i+1))^2)),(1/(((j+1))^2))];

(Subscript[Ξ, l_])[u_,R_]:=((Sqrt[l+1])/(4π (R^2)))JacobiP[l,1,0,1-(u/R)]

Hexa[R_]:=N[Table[((Subscript[A, i,j])/(R^2))+((Subscript[B, i,j])/R),{i,0,14},{j,0,14}],15];

ϵExa[R_]:=Eigenvalues[Hexa[R]]

Eigensystem[Hexa[100]][[2]][[15]]

{-0.82465661425964856288024366682528176888`15.,0.53727701978549676082504563075775565322`15.,-0.17452095311349187291869504119527529972`15.,0.02852169793776642744342094261106605818`15.,-0.00195498686012338440102052266792391892`15.,0.00001629228242546981128114106308795977`15.,5.2018856554108174480210646638262`15.*^-7,3.121239702824123259766471667307`15.*^-8,2.50554144291865954665130547175`15.*^-9,2.3891310367697317150336409614`15.*^-10,2.553050250326141645923770064`15.*^-11,2.95737510925779517642776513`15.*^-12,3.6369580904908151731488518`15.*^-13,4.694710201501790322291318`15.*^-14,7.0560748171665399151488049893920937170603743`15.*^-15}//Chop

ψ[u_,R_]:=Table[(Subscript[Ξ, n])[u,R],{n,0,14}].(Eigensystem[Hexa[R]][[2]][[15]]//Chop);

ψ[u,4]

ζ[θ1_,θ2_,ϕ1_,ϕ2_]:=π - ArcCos[Cos[θ1] Cos[θ2]+Cos[ϕ1-ϕ2] Sin[θ1] Sin[θ2]]

Sin[(π-ω)/2]

Sin[((π-ω)/2)]

Sin[π/2-ω]

With[{θ2=π/2,ϕ2=π/2,R=5},Plot3D[(((ψ[2R Cos[ζ[θ1,θ2,ϕ1,ϕ2]/2],R]))^2),{θ1,0,π},{ϕ1,0,2π}]]

With[{θ2=0,ϕ2=0,R=5},Plot3D[(((ψ[2R Cos[ζ[θ1,θ2,ϕ1,ϕ2]/2],R]))^2),{θ1,0,π},{ϕ1,0,2π}]]

With[{θ2=π/2,ϕ2=π/2,R=0.2},Plot3D[(((ψ[2R Cos[ζ[θ1,θ2,ϕ1,ϕ2]/2],R]))^2),{θ1,0,π},{ϕ1,0,2π}]]

With[{ϕ1=0,ϕ2=0,R=0.4},Plot3D[(((ψ[2R Cos[ζ[θ1,θ2,ϕ1,ϕ2]/2],R]))^2),{θ1,0,π},{θ2,0,π}]]

With[{ϕ1=0,ϕ2=0,R=20},Plot3D[(((ψ[2R Cos[ζ[θ1,θ2,ϕ1,ϕ2]/2],R]))^2),{θ1,0,π},{θ2,0,π}]]

The Weak Correlation Expansion

Energies on the real axis

{Ψ1=(Subscript[ψ, 1])[θ1](Subscript[ψ, 1])[θ2],Ψ2=(Subscript[ψ, 1])[θ1](Subscript[ψ, 2])[θ2],Ψ3=(Subscript[ψ, 2])[θ1](Subscript[ψ, 1])[θ2],Ψ4=(Subscript[ψ, 2])[θ1](Subscript[ψ, 2])[θ2]};

T[Ψ1_,Ψ2_]:=((((2π))^2)/(2(R^2)))(∫_0^π)(∫_0^π)Expand[((∂_(θ1))Ψ1)((∂_(θ1))Ψ2)+((∂_(θ2))Ψ1)((∂_(θ2))Ψ2)]Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2//FullSimplify
V[Ψ1_,Ψ2_]:=((((2π))^2)/R)(∫_0^π)(∫_0^π) Expand[Ψ1 Ψ2(1+Cos[θ1] Cos[θ2]+(1/4) (-1+3 ((Cos[θ1])^2)) (-1+3 ((Cos[θ2])^2)))] Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2//FullSimplify
S[Ψ1_,Ψ2_]:=(((2π))^2)(∫_0^π)(∫_0^π) Expand[Ψ1 Ψ2] Sin[θ1]Sin[θ2]\[DifferentialD]θ1\[DifferentialD]θ2//FullSimplify

{t,v,s}={(({{T[Ψ1,Ψ1], T[Ψ2,Ψ1], T[Ψ3,Ψ1], T[Ψ4,Ψ1]}, {T[Ψ1,Ψ2], T[Ψ2,Ψ2], T[Ψ3,Ψ2], T[Ψ4,Ψ2]}, {T[Ψ1,Ψ3], T[Ψ2,Ψ3], T[Ψ3,Ψ3], T[Ψ4,Ψ3]}, {T[Ψ1,Ψ4], T[Ψ2,Ψ4], T[Ψ3,Ψ4], T[Ψ4,Ψ4]}})),(({{V[Ψ1,Ψ1], V[Ψ2,Ψ1], V[Ψ3,Ψ1], V[Ψ4,Ψ1]}, {V[Ψ1,Ψ2], V[Ψ2,Ψ2], V[Ψ3,Ψ2], V[Ψ4,Ψ2]}, {V[Ψ1,Ψ3], V[Ψ2,Ψ3], V[Ψ3,Ψ3], V[Ψ4,Ψ3]}, {V[Ψ1,Ψ4], V[Ψ2,Ψ4], V[Ψ3,Ψ4], V[Ψ4,Ψ4]}})),(({{S[Ψ1,Ψ1], S[Ψ2,Ψ1], S[Ψ3,Ψ1], S[Ψ4,Ψ1]}, {S[Ψ1,Ψ2], S[Ψ2,Ψ2], S[Ψ3,Ψ2], S[Ψ4,Ψ2]}, {S[Ψ1,Ψ3], S[Ψ2,Ψ3], S[Ψ3,Ψ3], S[Ψ4,Ψ3]}, {S[Ψ1,Ψ4], S[Ψ2,Ψ4], S[Ψ3,Ψ4], S[Ψ4,Ψ4]}}))};

{t[R_],v[R_],s}:={(({{0, 0, 0, 0}, {0, (1/(R^2)), 0, 0}, {0, 0, (1/(R^2)), 0}, {0, 0, 0, (2/(R^2))}})),(({{(1/R), 0, 0, (1/(3 R))}, {0, (1/R), (1/(3 R)), 0}, {0, (1/(3 R)), (1/R), 0}, {(1/(3 R)), 0, 0, (29/(25 R))}})),(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))};

Hλ[λ_,R_]:=t[R]+λ v[R]
ϵWCE[λ_,R_]:=Eigenvalues[Hλ[λ,R]]

With[{R=2},Plot[Evaluate[ϵWCE[λ,R]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{WCE}",Magnification->1.5]},AspectRatio->1]]

Plot of the coefficients

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵWCE[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,20,1}],{R,{0.1,0.5,1}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.1,0.5,1},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵWCE[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,20,1}],{R,{0.5,1,2,3}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.5,1,2,3},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵWCE[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,20,1}],{R,{5,25,100,200}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{5,25,100,200},ImageSize->Large]

Singularity structure

Can we give a physical signification to all the singularities ?

P[λ_,ϵ_]=CharacteristicPolynomial[Hλ[λ,2],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30]
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

With[{R=2},
NSolve[ϵWCE[λ,R][[3]]==ϵWCE[λ,R][[4]],λ]]

There is a pair of complex conjugate branch point between (s^2)and ((Subscript[p, z])^2)

With[{R=2},
NSolve[ϵWCE[λ,R][[3]]==ϵWCE[λ,R][[2]],λ]]

There s a pair of complex conjugate branch point between (s^2)and the (Subscript[sp, z]) singlet (x2 degen)

With[{R=2},
NSolve[ϵWCE[λ,R][[3]]==ϵWCE[λ,R][[1]],λ]]

There is a CI between (s^2) and (Subscript[sp, z]) triplet

With[{R=2},
NSolve[ϵWCE[λ,R][[4]]==ϵWCE[λ,R][[1]],λ]]

There is a CI between  ((Subscript[p, z])^2) and (Subscript[sp, z]) triplet

With[{R=2},
NSolve[ϵWCE[λ,R][[4]]==ϵWCE[λ,R][[2]],λ]]

With[{R=2},
NSolve[ϵWCE[λ,R][[1]]==ϵWCE[λ,R][[2]],λ]]

R dependance

P[λ_,ϵ_]=CharacteristicPolynomial[Hλ[λ,2],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[Hλ[λ,3],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

Solve[ϵWCE[λ,R][[3]]==ϵWCE[λ,R][[4]],λ]

LogLinearPlot[{Abs[(((450/661)-((1875 I)/661))/R)],1},{R,0.1,1000},PlotRange->{0.96,25}]

RHF vs UHF: the Møller-Plesset expansion

Energies on the real axis

RHF

HRMP=DiagonalMatrix[{F[[1,1]]+F[[1,1]],F[[1,1]]+F[[2,2]],F[[2,2]]+F[[1,1]],F[[2,2]]+F[[2,2]]}]

HRMP[R_]:=(({{(2/R), 0, 0, 0}, {0, (1/(R^2))+(8/(3 R)), 0, 0}, {0, 0, (1/(R^2))+(8/(3 R)), 0}, {0, 0, 0, (2/(R^2))+(10/(3 R))}}));

We define the adiabatic connection Hamiltonian as

rHλ[λ_,R_]:=HRMP[R]+λ(rH[R]-HRMP[R]);
ϵRMP[λ_,R_]:=Eigenvalues[rHλ[λ,R]];
ϵRHF[R_]:=Eigenvalues[rH[R]]//FullSimplify;

With[{R=10},Plot[Evaluate[ϵRMP[λ,R]],{λ,-1.5,3},PlotStyle->{RGBColor["LimeGreen"],RGBColor["Red"],RGBColor["Blue"],Darker[Purple]},FrameStyle->Directive[Thick,18],PlotLegends->{Placed[{MaTeX["\\text{sp}_\\text{z} ~ ^{1}P",Magnification->1.5],MaTeX["\\text{sp}_\\text{z} ~ ^{3}P",Magnification->1.5],MaTeX["\\text{s}^2",Magnification->1.5],MaTeX["\\text{p}_\\text{z}^2",Magnification->1.5]},{Right,Top}]},AxesLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E",Magnification->1.5]},AxesStyle->Directive[Thick,16],AspectRatio->1,PlotRange->{{-1,2.8},{-0.3,0.5}}]]

(Legended[, (LineLegend[Singly excited state singlet (antisymmetric), Singly excited state triplet (antisymmetric), Ground state singlet (symmetric), Doubly excited state singlet (symmetric)])])

UHF

HUMP=DiagonalMatrix[{Fα[[1,1]]+Fβ[[1,1]],Fα[[1,1]]+Fβ[[2,2]],Fα[[2,2]]+Fβ[[1,1]],Fα[[2,2]]+Fβ[[2,2]]}]

HUMP[R_]:=(({{((75+118 R)/(84 (R^2))), 0, 0, 0}, {0, ((25+114 R)/(56 (R^2)))+((75+118 R)/(168 (R^2))), 0, 0}, {0, 0, ((25+114 R)/(56 (R^2)))+((75+118 R)/(168 (R^2))), 0}, {0, 0, 0, ((25+114 R)/(28 (R^2)))}}));

uHλ[λ_,R_]:=HUMP[R]+λ(uH[R]-HUMP[R]);
ϵUMP[λ_,R_]:=Eigenvalues[uHλ[λ,R]];

With[{R=1},Plot[Evaluate[Re[ϵUMP[λ,R]]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{UHF}",Magnification->1.5]},AspectRatio->1]]

With[{R=1},Show[ListPlot[Table[{λ,Sort[Re[ϵUMP[λ,R]]][[i]]},{i,1,4},{λ,-3,0,0.01}],PlotStyle->{RGBColor["Red"],RGBColor["Blue"],RGBColor["LimeGreen"],Darker[Purple]},PlotMarkers->Graphics@{Disk[{0,0},Scaled@0.001]},Joined->True],ListPlot[Table[{λ,Sort[Re[ϵUMP[λ,R]]][[i]]},{i,1,4},{λ,0,3,0.01}],PlotStyle->{RGBColor["Red"],RGBColor["LimeGreen"],RGBColor["Blue"],Darker[Purple]},PlotMarkers->Graphics@{Disk[{0,0},Scaled@0.001]},Joined->True,PlotLegends->{Placed[{MaTeX["\\text{s}^2",Magnification->1.5],MaTeX["\\text{sp}_\\text{z} ~ ^{3}P",Magnification->1.5],MaTeX["\\text{sp}_\\text{z} ~ ^{1}P",Magnification->1.5],MaTeX["\\text{p}_\\text{z}^2",Magnification->1.5]},{Right,Top}]}],PlotRange->All,AxesLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["\\text{Re}(E)",Magnification->1.5]},ImageSize->Large,AspectRatio->1,AxesStyle->Directive[Thick,16]]]

With[{R=1},Show[ListPlot[Table[{λ,Sort[Im[ϵUMP[λ,R]]][[i]]},{i,1,4},{λ,-3,0.6,0.01}],PlotStyle->{RGBColor["LimeGreen"],RGBColor["Red"],RGBColor["Blue"],Darker[Purple]},PlotMarkers->Graphics@{Disk[{0,0},Scaled@0.001]},Joined->True],ListPlot[Table[{λ,Sort[Im[ϵUMP[λ,R]]][[i]]},{i,1,4},{λ,0.8,3,0.01}],PlotStyle->{RGBColor["Red"],RGBColor["Blue"],Darker[Purple],RGBColor["LimeGreen"]},PlotMarkers->Graphics@{Disk[{0,0},Scaled@0.001]},Joined->True],ListPlot[Table[{λ,Sort[Im[ϵUMP[λ,1]]][[i]]},{i,1,4},{λ,0.6,0.8,0.01}],PlotStyle->{RGBColor["LimeGreen"],RGBColor["Red"],RGBColor["Blue"],Darker[Purple]},PlotMarkers->Graphics@{Disk[{0,0},Scaled@0.001]},Joined->True],PlotRange->All,AspectRatio->1,AxesLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["\\text{Im}(E)",Magnification->1.5]},ImageSize->Large,AxesStyle->Directive[Thick,16]]]

ListPlot[Table[{λ,Sort[Im[ϵUMP[λ,1]]][[i]]},{i,1,4},{λ,0.6,0.8,0.01}],PlotStyle->{RGBColor["LimeGreen"],RGBColor["Red"],RGBColor["Blue"],Darker[Purple]},PlotMarkers->Graphics@{Disk[{0,0},Scaled@0.001]},Joined->True]

With[{R=1},Plot[Evaluate[Im[ϵUMP[λ,R]]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{UHF}",Magnification->1.5]},AspectRatio->1]]

With[{R=1.5},Plot[Evaluate[Re[ϵUMP[λ,R]]],{λ,-1.5,3},PlotStyle->{RGBColor["Red"],RGBColor["LimeGreen"],RGBColor["Blue"],Darker[Purple]},FrameStyle->Directive[Thick,18],AxesLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E",Magnification->1.5]},AxesStyle->Directive[Thick,16],AspectRatio->1,ImageSize->Large,PlotRange->{{0,3},{-1,2}}]]

With[{R=1.51},Plot[Evaluate[Re[ϵUMP[λ,R]]],{λ,-1.5,3},PlotStyle->{RGBColor["Blue"],RGBColor["Red"],RGBColor["LimeGreen"],Darker[Purple]},FrameStyle->Directive[Thick,18],PlotLegends->{Placed[{MaTeX["\\text{s}^2",Magnification->1.5],MaTeX["\\text{sp}_\\text{z} ~ ^{3}P",Magnification->1.5],MaTeX["\\text{sp}_\\text{z} ~ ^{1}P",Magnification->1.5],MaTeX["\\text{p}_\\text{z}^2",Magnification->1.5]},{Right,Top}]},AxesLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E",Magnification->1.5]},AxesStyle->Directive[Thick,16],AspectRatio->1,ImageSize->Large,PlotRange->{{0,3},{-1,2}}]]

RHF -> one avoided crossing 
UHF -> more than one because of the spin contamination
The singularity structure will be really different.

uHλ[λ,R]//MatrixForm

With[{R=1.49},Plot[Evaluate[Re[ϵUMP[λ,R]]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{UHF}",Magnification->1.5]},AspectRatio->1]]

(Legended)(Legended)(Legended)

R=1.49                                                                                                          R=1.5                                                                                                           R=1.51

Eigensystem[uHλ[1,151]]//N//Chop

Plot of the coefficients

RHF

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,20,1}],{R,{0.1,0.5,1}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.1,0.5,1},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,50,1}],{R,{0.5,1,2,3}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.5,1,2,3},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,100,1}],{R,{5,25,100,200}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{5,25,100,200},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,100,1}],{R,{100,300,500}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{100,300,500},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,10,1}],{R,{100,300,500}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{100,300,500},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,300,1}],{R,{500}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{500},ImageSize->Large]

UHF

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUMP[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,20,1}],{R,{0.1,0.5,1}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.1,0.5,1,1.25},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUMP[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,100,1}],{R,{2,3,4}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{2,3,4},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUMP[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,100,1}],{R,{5,25,100,200}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{5,25,100,200},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUMP[λ,R][[2]],{λ,0.,n}]//Normal]},{n,20,300,1}],{R,{5}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{5},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUMP[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,100,1}],{R,{100,200,300}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{100,200,300},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUMP[λ,R][[2]],{λ,0,n}]//Normal]},{n,0,100,1}],{R,{500}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{500},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUMP[λ,R][[2]],{λ,0,n}]//Normal]},{n,0,250,1}],{R,{500}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{500},ImageSize->Large]

Total energy corrections

ListPlot[
Table[Evaluate[Table[{n,Series[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal},{n,0,20,1}]/.λ->1.],{R,{0.1,0.5,1}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.1,0.5,1},ImageSize->Large]

ListPlot[
Table[Evaluate[Table[{n,Series[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal},{n,0,100,1}]/.λ->1.],{R,{25,100,200}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{25,100,200},ImageSize->Large]

Singularity structure

Can we give a physical signification to all the singularities ?

RHF

P[λ_,ϵ_]=CharacteristicPolynomial[rHλ[λ,1],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30]
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

There are 12 singularities as expected, but we can’t give them a physical signification for the moment. Let’s try to do this for this simple case

With[{R=2},
NSolve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[4]],λ]]

There is a pair of complex conjugate branch point between (s^2)and ((Subscript[p, z])^2)

With[{R=2},
NSolve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[2]],λ]]

There are 2 CI between (s^2)and the (Subscript[sp, z]) triplet (x2 degen)

With[{R=1},
NSolve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[1]],λ]]

The (Subscript[sp, z]) singlet is never degenerate with the ground state

With[{R=2},
NSolve[ϵRMP[λ,R][[4]]==ϵRMP[λ,R][[1]],λ]]

There is a pair of complex conjugate branch point between ((Subscript[p, z])^2) and (Subscript[sp, z]) singlet (x2 degen)

With[{R=2},
NSolve[ϵRMP[λ,R][[4]]==ϵRMP[λ,R][[2]],λ]]

With[{R=2},
NSolve[ϵRMP[λ,R][[1]]==ϵRMP[λ,R][[2]],λ]]

UHF

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ/.R->2,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Sort[λEP]
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

With[{R=2},
NSolve[ϵUMP[λ,R][[1]]==ϵUMP[λ,R][[4]],λ]]

There is a pair of complex conjugate branch point between (s^2)and ((Subscript[p, z])^2)

With[{R=2},
NSolve[ϵUMP[λ,R][[1]]==ϵUMP[λ,R][[2]],λ]]

With[{R=2},
NSolve[ϵUMP[λ,R][[1]]==ϵUMP[λ,R][[3]],λ]]

With[{R=2},
NSolve[ϵUMP[λ,R][[2]]==ϵUMP[λ,R][[3]],λ]]

With[{R=2},
NSolve[ϵUMP[λ,R][[2]]==ϵUMP[λ,R][[4]],λ]]

With[{R=2},
NSolve[ϵUMP[λ,R][[3]]==ϵUMP[λ,R][[4]],λ]]

R dependance

RHF

P[λ_,ϵ_]=CharacteristicPolynomial[rHλ[λ,1],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[rHλ/.R->10,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[rHλ/.R->20,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

We see that this is the EP between (s^2)and ((Subscript[p, z])^2) that control the convergence

Solve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[4]],λ]

Show[LogLinearPlot[Abs[((25 (88+(132/R)+((25 (Sqrt[-9-12 R-4 (R^2)]))/R)))/2561)],{R,1,2000},PlotRange->{0.96,1.06},PlotTheme->"Scientific",FrameStyle->Directive[Thick,16],FrameLabel->{MaTeX["R",Magnification->1.5],MaTeX["R_{\\text{CV}}",Magnification->1.5]}],LogLinearPlot[1,{R,1,2000},PlotStyle->{Dotted,Black}],AspectRatio->1]

LogLinearPlot[{Abs[((25 (88+(132/R)+((25 (Sqrt[-9-12 R-4 (R^2)]))/R)))/2561)],1},{R,0.1,1000},PlotRange->{0.96,25}]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal]},{n,0,100,1}],{R,{100,200}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{100,200},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal]},{n,0,400,1}],{R,{100,200}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{100,200},ImageSize->Large]

UHF

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,-2],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
λEP
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-5,5},{-5,5}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

With[{R=10},Plot[Evaluate[Re[ϵUMP[λ,R]]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{UHF}",Magnification->1.5]},AspectRatio->1]]

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,1.5],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
λEP
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,-75/62-0.1],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Sort[λEP]
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,2.6],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,5],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Sort[λEP];
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-0.5,0.5}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,10],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Sort[λEP]
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-0.2,0.2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ/.R->100,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Sort[λEP];
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,500],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Sort[λEP]
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

{{1,0.2374891467337787324088047102385541022612789519696924509021`29.057249529911832-0.8028546019039053543517392040254700532679904648679610977402`29.397940008672037 I},{1.1,0.4103602476827904930974827867820610912316749461629579191583`29.397940008672034-0.9199967823451507732947378036217367894690103625221676596362`29.382045606050877 I},{1.2,0.6325132614340338064606724797757945618337008173994847866943`29.397940008672037-0.9969463513600492866805439325006681051550875709224992364534`29.37657862532185 I},{1.3,0.8983661464806765204328282865490136810518356955335865118411`29.397940008672037-1.0117986604394248391677330434612120547211143805339788696858`29.338150095970537 I},{1.4,1.2132194606214072560152425658109284555296301326112167278446`29.358978980210164-0.9398425629846078013761494783848553304053852385701496357056`29.397940008672037 I},{1.5,1.7180788754388699718668967211291306945962368709039930005602`29.397940008672037+0.9761811792275888495142615075427764654293363424569219777264`29.218913175720342 I},{1.6,1.2213127273900431528823334673558248119938342390482437008927`29.397940008672037-0.2341254632770232770262092034821885597194057838759222967847`28.902119336225294 I},{1.7,1.1462056976220943513227441662158693646748908075043656549443`29.397940008672037-0.2049914959325541227299672600454722499046824129731852508417`28.87426974111576 I},{1.8,1.1070228917699778321081079490164397513201721157339728599782`29.397940008672037-0.1797376523337093194544259341484853513277700915761132064295`28.842045645066868 I},{2,1.0655678648783681952455499596090965480045810739079991800999`29.397940008672037-0.1414990387982489895023924064789054248077506698461980555895`28.793316676554895 I}}//Abs

ListPlot[{{1,0.8372434571942719174709848054876339028996391379036395056783`29.35818958434467},{1.1`,1.0073676649583864253023301680161145856827842190835168673413`29.384643216109342},{1.2`,1.1806672068707732360145619199598521452833111841474963223604`29.38260262929695},{1.3`,1.35306993995489932677509118723192298330639713334041406073`29.36349874541293},{1.4`,1.5346678803011336999024836477931887900818194739419630573141`29.37318460654532},{1.5`,1.9760376304381093501790635490498365642526666691994111520203`29.346978451060966},{1.6`,1.2435511692888180617058862112801644953777287893846624120874`29.366300894915515},{1.7`,1.164392122382326297549203881433715740881463333262773216045`29.367539851218485},{1.8`,1.1215191958095043752326967028200482032456150494419962200788`29.369901027080203},{2,1.074921789081546487205347181983407558609538486715624845771`29.375761277750755}}]

With[{R=1.5},P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,R],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[Sort[λEP]];
Print[Plot[Evaluate[Re[ϵUMP[λ,R]]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{UHF}",Magnification->1.5]},AspectRatio->1]]];
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

Eigensystem[uHλ[1,1]]//N

Riemann surfaces

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,2],ϵ];
NRJ =ϵ/. {ToRules[Roots[P[λ,ϵ]==0,ϵ]]};
Plot3D[Evaluate@Abs[NRJ/.λ-> λx+I λy],{λx,-2,2},{λy,-2,2},PlotLegends->Automatic,AxesLabel->{"Re(λ)","Im(λ)"},AxesStyle->Directive[16],BoxRatios->{1, 1, 1},Exclusions->{Automatic,Offset-> 0.05},Mesh->None,PlotStyle->Opacity[0.7]]

P[λ_,ϵ_]=CharacteristicPolynomial[uHλ[λ,5],ϵ];
NRJ =ϵ/. {ToRules[Roots[P[λ,ϵ]==0,ϵ]]};
Plot3D[Evaluate@Re[NRJ/.λ-> λx+I λy],{λx,-1,3},{λy,-1,1},PlotLegends->Automatic,AxesLabel->{"Re(λ)","Im(λ)"},AxesStyle->Directive[16],BoxRatios->{1, 1, 1},Exclusions->{Automatic,Offset-> 0.05},Mesh->None,PlotStyle->Opacity[0.7],ImageSize->Large]

Archetypes

RHF

rHλ[1,R]//FullSimplify//MatrixForm

rHλ[λ,R]//FullSimplify//MatrixForm

(29/25)+(163/75)==(250/75)

γ=(163/(75R));
δ=(1/(3R));

{γ//N,δ//N}

Abs[γ] > Abs[δ] but not Abs[γ] >> Abs[δ]: this is a symetric ripple archetype, this agree with the Plot of the coefficients part.

UHF

Clear[γ,δ]

uHλ[1,R]//FullSimplify//MatrixForm

uHλ[λ,R]//FullSimplify//MatrixForm

β=((16875+4 R (2775+1511 R))/(8400 (R^3)))

((16875 λ+4 R (1875+8550 R+900 λ-7039 R λ))/(8400 (R^3)))//Expand

β+γ-λ γ=(25/(28 (R^2)))+(57/(14 R))+((225 λ)/(112 (R^3)))+((3 λ)/(7 (R^2)))-((7039 λ)/(2100 R))

(25/(28 (R^2)))+(57/(14 R))-((16875+4 R (2775+1511 R))/(8400 (R^3)))//FullSimplify//Expand

γ[R_]:=-(225/(112 (R^3)))-(3/(7 (R^2)))+(7039/(2100 R))

(((-75+4 (-3+R) R))/(112 (R^3)))//Expand

δ[R_]:=-(75/(112 (R^3)))-(3/(28 (R^2)))+(1/(28 R))

LogLinearPlot[{γ[R],δ[R]},{R,0.5,500},PlotLegends->"Expressions"]

RMP2 vs UMP2

With[{R=1},{Series[ϵRMP[λ,R][[3]],{λ,0,2}]//Normal,Series[ϵUMP[λ,R][[2]],{λ,0,2}]//Normal}/.λ->1//N]

Table[With[{n=2},{R,Series[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal,Series[ϵUMP[λ,R][[2]],{λ,0,n}]//Normal,ϵExa[R]//Sort//First}]/.λ->1//N,{R,{0.1,0.2,0.5,1,2,3,4,5,10,20,50,100,1000}}]//TableForm

Table[With[{n=2},{R,(Series[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal)-(ϵExa[R]//Sort//First),(Series[ϵUMP[λ,R][[2]],{λ,0,n}]//Normal)-(ϵExa[R]//Sort//First)}]/.λ->1//N,{R,{0.1,0.2,0.5,1,2,3,4,5,10,20,50,100,1000}}]//TableForm

Table[With[{n=2},{R,(Series[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal)/(ϵExa[R]//Sort//First),(Series[ϵUMP[λ,R][[2]],{λ,0,n}]//Normal)/(ϵExa[R]//Sort//First)}]/.λ->1//N,{R,{0.1,0.2,0.5,1,2,3,4,5,10,20,50,100,1000}}]//TableForm

Table[With[{n=2},{R,((Series[ϵRMP[λ,R][[3]],{λ,0,n}]//Normal)-(ϵExa[R]//Sort//First))/(ϵExa[R]//Sort//First),((Series[ϵUMP[λ,R][[2]],{λ,0,n}]//Normal)-(ϵExa[R]//Sort//First))/(ϵExa[R]//Sort//First)}]/.λ->1//N,{R,{0.1,0.2,0.5,1,2,3,4,5,10,20,50,100,1000}}]//TableForm

Complete basis

RHF

We can restrict the basis to the CSF/Legendre polynomia for the RHF case.

HRMPn[R_,n_]:=DiagonalMatrix[Table[((l(l+1))/(R^2))+2((2/R)-(1/((2 l+1)R))),{l,0,n-1}]]
rHn[R_,n_]:=Table[ KroneckerDelta[i,j]((i(i+1))/(R^2))+TEIr[{i,0},{j,0},{i,0},{j,0},R],{i,0,n-1},{j,0,n-1}]
rHnλ[λ_,R_,n_]:=HRMPn[R,n]+λ(rHn[R,n]-HRMPn[R,n]);
ϵRMPn[λ_,R_,n_]:=Eigenvalues[rHnλ[λ,R,n]];

With[{R=10,n=5},Plot[Evaluate[ϵRMPn[λ,R,n]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{RHF}",Magnification->1.5]},AspectRatio->1]]

P[λ_,ϵ_]=CharacteristicPolynomial[rHnλ[λ,20,5]//N,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

PlotRcv[n_]:=Module[{HMP,H},
HMP[R_]:= HRMPn[R,n];
H[R_]:=rHn[R,n];
Hλ[R_,λ_]:=HMP[R]+λ(H[R]-HMP[R]);
ListR=Table[{R,Rcv[Hλ,R]},{R,logspace[0,3,100]}];
ListLogLinearPlot[ListR]
]

logspace[a_,b_,n_]:=10.0^Range[a,b,(b-a)/(n-1)]

Rcv[H_,R_]:=Module[{λ,ϵ,λEP},
P[λ_,ϵ_]=CharacteristicPolynomial[H[R,λ]//N,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Return[λEP//Abs//Sort//First]
]

PlotRcv[2]

UHF

“Final” results

RHF

With[{R=2},Plot[Evaluate[ϵRMP[λ,R]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{RHF}",Magnification->1.5]},AspectRatio->1]]

Manipulate[With[{R=a},Plot[Evaluate[ϵRMP[λ,R]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{RHF}",Magnification->1.5]},AspectRatio->1]],{a,0.1,5}]

P[λ_,ϵ_]=CharacteristicPolynomial[rHλ[λ,2],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30]
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

There are 12 singularities as expected, but we can’t give them a physical signification for the moment. Let’s try to do this for this simple case

With[{R=2},
NSolve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[4]],λ]]

There is a pair of complex conjugate branch point between (s^2)and ((Subscript[p, z])^2)

With[{R=2},
NSolve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[2]],λ]]

There are 2 CI between (s^2)and the (Subscript[sp, z]) triplet (x2 degen)

With[{R=1},
NSolve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[1]],λ]]

The (Subscript[sp, z]) singlet is never degenerate with the ground state

With[{R=2},
NSolve[ϵRMP[λ,R][[4]]==ϵRMP[λ,R][[1]],λ]]

There is a pair of complex conjugate branch point between ((Subscript[p, z])^2) and (Subscript[sp, z]) singlet (x2 degen)

With[{R=2},
NSolve[ϵRMP[λ,R][[4]]==ϵRMP[λ,R][[2]],λ]]

With[{R=2},
NSolve[ϵRMP[λ,R][[1]]==ϵRMP[λ,R][[2]],λ]]

This is the avoided crossing between (s^2)and ((Subscript[p, z])^2) which dictates the convergence properties for the RHF reference.

Solve[ϵRMP[λ,R][[3]]==ϵRMP[λ,R][[4]],λ]

LogLinearPlot[{Abs[((25 (88+(132/R)+((25 (Sqrt[-9-12 R-4 (R^2)]))/R)))/2561)],1},{R,1,1000},PlotRange->{0.96,1.06},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["R",Magnification->1.5],MaTeX["R_\\text{CV}",Magnification->1.5]},AspectRatio->1]

UHF

For R>(3/2) a broken symmetry solution is possible. By allowing the singly excited triplet state to mix with the (s^2)and ((Subscript[p, z])^2) states one can get a solution with a lower energy than in the RHF formalism.

This is the sharp avoided crossing between (s^2) and (Subscript[sp, z]) that dictates the convergence properties. Just after R=(3/2) the CI become a sharp avoided crossing which try to model the symmetry breaking of the real E(λ) function.

RHF vs UHF: the Epstein-Nesbet expansion

Energies on the real axis

RHF

We define the adiabatic connection Hamiltonian as

HRENλ[λ_,R_]:=DiagonalMatrix[Diagonal[rH[R]]]+λ (rH[R]-DiagonalMatrix[Diagonal[rH[R]]]);
ϵREN[λ_,R_]:=Eigenvalues[HRENλ[λ,R]];

With[{R=5},Plot[Evaluate[ϵREN[λ,R]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{RHF}",Magnification->1.5]},AspectRatio->1]]

UHF

HUENλ[λ,R]

We define the adiabatic connection Hamiltonian as

HUENλ[λ_,R_]:=DiagonalMatrix[Diagonal[uH[R]]]+λ (uH[R]-DiagonalMatrix[Diagonal[uH[R]]]);
ϵUEN[λ_,R_]:=Eigenvalues[HUENλ[λ,R]];

With[{R=10},Plot[Evaluate[ϵUEN[λ,R]],{λ,-3,3},FrameStyle->Directive[Thick,18],PlotTheme->"Scientific",PlotLegends->Automatic,FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E_\\text{UHF}",Magnification->1.5]},AspectRatio->1]]

Plot of the coefficients

RHF

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵREN[λ,R][[3]],{λ,0.,n}]//Normal]},{n,0,5,1}],{R,{0.1,0.5,1}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.1,0.5,1},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵREN[λ,R][[3]],{λ,0,n}]//Normal]},{n,0,10,1}],{R,{2,3,4}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{2,3,4},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵREN[λ,R][[3]],{λ,0,n}]//Normal]},{n,0,20,1}],{R,{5,25,100}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{5,25,100},ImageSize->Large]

UHF

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUEN[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,20,1}],{R,{0.1,0.5,1}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{0.1,0.5,1},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUEN[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,40,1}],{R,{2,3,4}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{2,3,4},ImageSize->Large]

5,25,100

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUEN[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,100,1}],{R,{5,10,25}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{5,10,25},ImageSize->Large]

ListPlot[
Table[Table[{n,N[SeriesCoefficient[ϵUEN[λ,R][[2]],{λ,0.,n}]//Normal]},{n,0,500,1}],{R,{2}}]
,Joined->True,PlotMarkers->{Automatic,Small},PlotRange->Automatic,
PlotLegends->{2},ImageSize->Large]

Singularity structure

Can we give a physical signification to all the singularities ?

RHF

P[λ_,ϵ_]=CharacteristicPolynomial[HRENλ/.R->2,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30]
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

There are only 8 singularities, but we can’t give them a physical signification for the moment. Let’s try to do this for this simple case

With[{R=2},
NSolve[ϵREN[λ,R][[3]]==ϵREN[λ,R][[4]],λ]]

There is a pair of complex conjugate branch point between (s^2)and ((Subscript[p, z])^2)

With[{R=2},
NSolve[ϵREN[λ,R][[3]]==ϵREN[λ,R][[2]],λ]]

There is 1 CI between (s^2)and the (Subscript[sp, z]) triplet

With[{R=2},
NSolve[ϵREN[λ,R][[3]]==ϵREN[λ,R][[1]],λ]]

The (Subscript[sp, z]) singlet is never degenerate with the ground state

With[{R=2},
NSolve[ϵREN[λ,R][[4]]==ϵREN[λ,R][[2]],λ]]

Solve[ϵREN[λ,R][[4]]==ϵREN[λ,R][[2]],λ]

With[{R=2},
NSolve[ϵREN[λ,R][[4]]==ϵREN[λ,R][[1]],λ]]

With[{R=2},
NSolve[ϵREN[λ,R][[1]]==ϵREN[λ,R][[2]],λ]]

UHF

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ/.R->2,ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Sort[λEP]
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-4,4},{-4,4}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

With[{R=2},
NSolve[ϵUEN[λ,R][[1]]==ϵUEN[λ,R][[4]],λ]]

There is a pair of complex conjugate branch point between (s^2)and ((Subscript[p, z])^2)

With[{R=2},
NSolve[ϵUEN[λ,R][[1]]==ϵUEN[λ,R][[2]],λ]]

With[{R=2},
NSolve[ϵUEN[λ,R][[1]]==ϵUEN[λ,R][[3]],λ]]

With[{R=2},
NSolve[ϵUEN[λ,R][[2]]==ϵUEN[λ,R][[3]],λ]]

With[{R=2},
NSolve[ϵUEN[λ,R][[2]]==ϵUEN[λ,R][[4]],λ]]

With[{R=2},
NSolve[ϵUEN[λ,R][[3]]==ϵUEN[λ,R][[4]],λ]]

R dependance

RHF

P[λ_,ϵ_]=CharacteristicPolynomial[HRENλ[λ,1],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[HRENλ[λ,2],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[HRENλ[λ,5],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

UHF

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ[λ,1.8],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-5,5},{-5,5}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ[λ,2],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ[λ,5],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Abs[λEP]//Sort
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ[λ,8],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ[λ,10],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ[λ,20],ϵ];
dPdϵ[λ_,ϵ_]=D[P[λ,ϵ],ϵ];
EP=NSolve[{P[λ,ϵ]==0,dPdϵ[λ,ϵ]==0},{ϵ,λ},WorkingPrecision->30];
{λEP=λ/.EP,ϵEP=ϵ/.EP};
Print[ListPlot[Table[{Re[λEP[[k]]],Im[λEP[[k]]]},{k,Length[λEP]}],PlotRange->{{-2,2},{-2,2}},Epilog->{Green,Opacity[0.5],Disk[{0,0}]},PlotMarkers->{Automatic,Medium},PlotTheme->"Scientific",GridLines->Automatic,BaseStyle->Directive[Thick,18],FrameLabel->{"Re[λ]","Im[λ]"}]];

Riemann surfaces

P[λ_,ϵ_]=CharacteristicPolynomial[HUENλ[λ,2],ϵ];
NRJ =ϵ/. {ToRules[Roots[P[λ,ϵ]==0,ϵ]]};
Plot3D[Evaluate@Abs[NRJ/.λ-> λx+I λy],{λx,-2,2},{λy,-2,2},PlotLegends->Automatic,AxesLabel->{"Re(λ)","Im(λ)"},AxesStyle->Directive[16],BoxRatios->{1, 1, 1},Exclusions->{Automatic,Offset-> 0.05},Mesh->None,PlotStyle->Opacity[0.7]]

Archetypes

RHF

HRENλ[1,R]//FullSimplify//MatrixForm

HRENλ[λ,R]//FullSimplify//MatrixForm

γ=0;
δ=(1/(3R));

Abs[δ] >> Abs[γ] : this is a symetric zigzag archetype, this agree with the Plot of the coefficients part.

UHF

Clear[γ,δ]

HUENλ[1,R]//FullSimplify//MatrixForm

HUENλ[λ,R]//FullSimplify//MatrixForm

β=((16875+4 R (2775+1511 R))/(8400 (R^3)))

((16875 λ+4 R (1875+8550 R+900 λ-7039 R λ))/(8400 (R^3)))//Expand

β+γ-λ γ=(25/(28 (R^2)))+(57/(14 R))+((225 λ)/(112 (R^3)))+((3 λ)/(7 (R^2)))-((7039 λ)/(2100 R))

(25/(28 (R^2)))+(57/(14 R))-((16875+4 R (2775+1511 R))/(8400 (R^3)))//FullSimplify//Expand

γ[R_]:=-(225/(112 (R^3)))-(3/(7 (R^2)))+(7039/(2100 R))

(((-75+4 (-3+R) R))/(112 (R^3)))//Expand

δ[R_]:=-(75/(112 (R^3)))-(3/(28 (R^2)))+(1/(28 R))

LogLinearPlot[{γ[R],δ[R]},{R,0.5,500},PlotLegends->"Expressions"]