Initialization

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

J3[i_,a_,l_]:=ThreeJSymbol[{i,0},{a,0},{l,0}];

J6[{a_,b_,c_},{e_,f_,g_}]:=SixJSymbol[{a,b,c},{e,f,g}];

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

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

(Subscript[P, ℓ_])[θ_]=LegendreP[ℓ,θ];

(Subscript[𝒴, ℓ_,m_])[θ_,ϕ_]:=FullSimplify[Which[m<0,(I/(Sqrt[2]))((Subscript[Y, ℓ,m])[θ,ϕ]-(((-1))^m)(Subscript[Y, ℓ,-m])[θ,ϕ]),m==0,(Subscript[Y, ℓ,0])[θ,ϕ],m>0,(1/(Sqrt[2]))((Subscript[Y, ℓ,m])[θ,ϕ]+(((-1))^m)(Subscript[Y, ℓ,-m])[θ,ϕ])],0<θ<π&&0<ϕ<2π];

(Subscript[δ, a_,b_])=KroneckerDelta[a,b];

(Subscript[δ, x_])=KroneckerDelta[x];

Table[(Subscript[ψ, ℓ(ℓ+1)+m+1])[θ_,ϕ_]=(Subscript[𝒴, ℓ,m])[θ,ϕ],{ℓ,0,5},{m,-ℓ,ℓ}];

SortEigensystem[eigsys_]:=Transpose[Sort[Transpose[eigsys]]]

YY[{ℓ1_,m1_},{ℓ2_,m2_}]:=(∑_(ℓ3=Abs[ℓ1-ℓ1])^(ℓ1+ℓ2))(∑_(m3=-ℓ3)^(ℓ3))(Sqrt[(((2ℓ1+1)(2ℓ2+1)(2ℓ3+1))/(4π))])J3[ℓ1,ℓ2,ℓ3]ThreeJSymbol[{ℓ1,m1},{ℓ2,m2},{ℓ3,m3}](((-1))^m3)(Subscript[Y, ℓ3,-m3])

Transient uniform electron gases

Useful relations

With[{L=1},Table[(∫_0^π)(∫_0^(2π))(Subscript[Y, (Subscript[ℓ, 1]),(Subscript[m, 1])])[θ,ϕ](Subscript[Y, (Subscript[ℓ, 2]),(Subscript[m, 2])])[θ,ϕ](Subscript[Y, (Subscript[ℓ, 3]),(Subscript[m, 3])])[θ,ϕ]Sin[θ]\[DifferentialD]ϕ\[DifferentialD]θ==(Sqrt[(((2(Subscript[ℓ, 1])+1)(2(Subscript[ℓ, 2])+1)(2(Subscript[ℓ, 3])+1))/(4π))])J3[(Subscript[ℓ, 1]),(Subscript[ℓ, 2]),(Subscript[ℓ, 3])]ThreeJSymbol[{(Subscript[ℓ, 1]),(Subscript[m, 1])},{(Subscript[ℓ, 2]),(Subscript[m, 2])},{(Subscript[ℓ, 3]),(Subscript[m, 3])}]//Simplify,{(Subscript[ℓ, 1]),0,L},{(Subscript[ℓ, 2]),0,L},{(Subscript[ℓ, 3]),0,L},{(Subscript[m, 1]),-(Subscript[ℓ, 1]),(Subscript[ℓ, 1])},{(Subscript[m, 2]),-(Subscript[ℓ, 2]),(Subscript[ℓ, 2])},{(Subscript[m, 3]),-(Subscript[ℓ, 3]),(Subscript[ℓ, 3])}]]

With[{L=2},Table[(Subscript[Y, (Subscript[ℓ, 1]),(Subscript[m, 1])])[θ,ϕ](Subscript[Y, (Subscript[ℓ, 2]),(Subscript[m, 2])])[θ,ϕ]==(∑_((Subscript[ℓ, 3])=Abs[(Subscript[ℓ, 1])-(Subscript[ℓ, 2])])^((Subscript[ℓ, 1])+(Subscript[ℓ, 2])))(∑_((Subscript[m, 3])=-(Subscript[ℓ, 3]))^(Subscript[ℓ, 3]))(Sqrt[(((2(Subscript[ℓ, 1])+1)(2(Subscript[ℓ, 2])+1)(2(Subscript[ℓ, 3])+1))/(4π))])J3[(Subscript[ℓ, 1]),(Subscript[ℓ, 2]),(Subscript[ℓ, 3])]ThreeJSymbol[{(Subscript[ℓ, 1]),(Subscript[m, 1])},{(Subscript[ℓ, 2]),(Subscript[m, 2])},{(Subscript[ℓ, 3]),(Subscript[m, 3])}](((-1))^(Subscript[m, 3]))(Subscript[Y, (Subscript[ℓ, 3]),-(Subscript[m, 3])])[θ,ϕ]//Simplify,{(Subscript[ℓ, 1]),0,L},{(Subscript[ℓ, 2]),0,L},{(Subscript[m, 1]),-(Subscript[ℓ, 1]),(Subscript[ℓ, 1])},{(Subscript[m, 2]),-(Subscript[ℓ, 2]),(Subscript[ℓ, 2])}]]

With[{L=5},Table[(Subscript[Y, (Subscript[ℓ, 1])])[θ](Subscript[Y, (Subscript[ℓ, 2])])[θ]==(∑_((Subscript[ℓ, 3])=Abs[(Subscript[ℓ, 1])-(Subscript[ℓ, 2])])^((Subscript[ℓ, 1])+(Subscript[ℓ, 2])))(Sqrt[(((2(Subscript[ℓ, 1])+1)(2(Subscript[ℓ, 2])+1)(2(Subscript[ℓ, 3])+1))/(4π))])((J3[(Subscript[ℓ, 1]),(Subscript[ℓ, 2]),(Subscript[ℓ, 3])])^2)(Subscript[Y, (Subscript[ℓ, 3])])[θ]//Simplify,{(Subscript[ℓ, 1]),0,L},{(Subscript[ℓ, 2]),0,L}]]

Table[(2ℓ+1)(∫_0^π)(Subscript[P, ℓ])[Cos[θ]] (((2R Sin[θ/2]))^n)Sin[θ]\[DifferentialD]θ==(((2R))^n)(((2ℓ+1)(((-1))^ℓ)2 ((n/2))!((n/2))!)/( ((n/2)+ℓ+1)!((n/2)-ℓ)!)),{ℓ,0,2},{n,0,10}]//MatrixForm

Jacobian on a d-sphere

For two particles on a d-sphere, we need to convert the Jacobian in the principal angle θ to a Jacobian in u

u=2R Sin[(θ/2)]	⟺	θ=2ArcSin[(u/(2R))]

((((Sin[θ])^(d-1))\[DifferentialD]θ)/((∫_0^π)((Sin[θ])^(d-1))\[DifferentialD]θ))=((Γ[((d+1)/2)])/((Sqrt[π]) Γ[(d/2)]))((Sin[θ])^(d-1))((\[DifferentialD]θ)/(\[DifferentialD]u))\[DifferentialD]u

=((Γ[((d+1)/2)])/((Sqrt[π]) Γ[(d/2)]))((Sin[2ArcSin[(u/(2R))]])^(d-1))(\[DifferentialD]/(\[DifferentialD]u))(2ArcSin[(u/(2R))])\[DifferentialD]u

=((Γ[((d+1)/2)])/((Sqrt[π]) Γ[(d/2)]R))((((u/R)))^(d-1))(((1-((((u/(2R))))^2)))^((d-2)/2))\[DifferentialD]u

Spherium

Legendre basis set

-(1/(R^2))(((\[DifferentialD]^2)/(\[DifferentialD](θ^2)))+(2Cot[θ]+(1/(Sin[θ])))(\[DifferentialD]/(\[DifferentialD]θ)))+(1/(2R Sin[θ/2]))+(1/(R^2)) 	(* Hamiltonian for the triplet P state *)

(Subscript[ψ, ℓ_])[θ_]=((Sqrt[2ℓ+1])/(4π))(Subscript[P, ℓ])[Cos[θ]];	With[{M=3},Table[8(π^2)(∫_0^π)(Subscript[ψ, i])[θ](Subscript[ψ, j])[θ]Sin[θ]\[DifferentialD]θ,{i,0,M},{j,0,M}]==IdentityMatrix[M+1]]

With[{M=2},H[α_]=Table[8(π^2)(∫_0^π)(Subscript[ψ, i])[θ](-((Subscript[ψ, j]^″)[θ]+(2Cot[θ]+(1/(Sin[θ])))(Subscript[ψ, j]^′)[θ])+((α/(2 Sin[θ/2]))+1)(Subscript[ψ, j])[θ])Sin[θ]\[DifferentialD]θ,{i,0,M},{j,0,M}]];

With[{M=2},
Table[
{ϵ,c}=SortEigensystem[Eigensystem[N[H[α],15]]];
f=Table[(Subscript[ψ, n])[θ],{n,0,M}];
ρ=Table[8(π^2)(2ℓ+1)(∫_0^π)(Subscript[P, ℓ])[Cos[θ]](((c[[1]].f))^2)Sin[θ]\[DifferentialD]θ,{ℓ,0,2}];
Print["α: ",N[α]];
Print["ρ: ",ρ];
Print[{ρ[[1]]-(2/3)ρ[[2]]+(1/5)ρ[[3]]}];
,{α,-10,0,1}]
]

Explicitly-correlated calculations for ground state

With[{M=10},
tab=Table[
{ϵ,c}=SortEigensystem[Eigensystem[{N[H[R],100],N[S[R],100]}]];
Print[ϵ];
f=Table[(Subscript[ψ, n])[ R (Sqrt[2-2Cos[θ]])],{n,M}];
ρ=Table[(ℓ+(1/2))NIntegrate[(Subscript[P, ℓ])[Cos[θ]](((c[[1]].f))^2)Sin[θ],{θ,0,π},WorkingPrecision->15],{ℓ,0,2}];
(*Print["R: ",N[R]];
Print["ρ: ",ρ];*)
{N[R],Chop[ρ[[1]]-(2/3)ρ[[2]]+(1/5)ρ[[3]]]}
,{R,{-Rationalize[5.3253]}}]
]
ListPlot[tab,Joined->True,PlotRange->All]

(* R = 5.3253 *)

(Subscript[c, 0])-(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2])==0

(Subscript[c, ℓ])∝(2ℓ+1)(∫_0^π)(Subscript[P, ℓ])[Cos[θ]]((ψ[u])^2)Sin[θ]\[DifferentialD]θ

((ψ[u])^2)=(∑_(n=0)^∞)(Subscript[a, n])(u^n)

(Subscript[c, 0])-(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2])==0

a=CoefficientList[1.` +0.6666467666189491` u+0.233711497250696` (u^2)+0.057892269462279114` (u^3)+0.011416060718727025` (u^4)+0.0019001512003665785` (u^5)+0.00027454051463773074` (u^6)+0.00003485368552570877` (u^7)+3.909550758995941`*^-6 (u^8)+3.886613919337475`*^-7 (u^9)+3.423932142765105`*^-8 (u^10)+2.6560194575496133`*^-9 (u^11)+1.7871244353862165`*^-10 (u^12)+1.017581390449739`*^-11 (u^13)+4.734306593585661`*^-13 (u^14)+1.714745687814212`*^-14 (u^15)+4.501929931069031`*^-16 (u^16)+7.574161874717616`*^-18 (u^17)+6.100436598427603`*^-20 (u^18),u];

NSolve[(∑_(n=0)^(Length[a]-1))a[[n+1]]((((2^(1+n)) (R^n) (n/2)!)/((1+(n/2))!))-(2/3)(-((3 (2^(1+n)) (R^n) ((((n/2)!))^2))/((-1+(n/2))! (2+(n/2))!)))+(1/5)((5 (2^(1+n)) (R^n) ((((n/2)!))^2))/((-2+(n/2))! (3+(n/2))!)))==0,R,Reals]

With[{M=2},
tab=Table[
{ϵ,c}=SortEigensystem[Eigensystem[{N[H[R],100],N[S[R],100]}]];
f=Table[(Subscript[ψ, n])[2 R Sin[θ/2]],{n,M}];
ρ=Table[(ℓ+(1/2))NIntegrate[(Subscript[P, ℓ])[Cos[θ]](((c[[1]].f))^2)Sin[θ],{θ,0,π},WorkingPrecision->15],{ℓ,0,2}];
ρ=Normalize[ρ];
(*Print["R: ",N[R]];
Print["ρ: ",ρ];*)
{N[R],Chop[ρ[[1]]-(2/3)ρ[[2]]+(1/5)ρ[[3]]]}
,{R,-8,-1}]
]
ListPlot[tab,Joined->True]

Table[
ρ=Table[(ℓ+(1/2))NIntegrate[(Subscript[P, ℓ])[Cos[θ]](((1+2R Sin[θ/2]/3))^2)Sin[θ],{θ,0,π},WorkingPrecision->15],{ℓ,0,2}];
ρ=Normalize[ρ];
(*Print["R: ",N[R]];
Print["ρ: ",ρ];*)
{N[R],Chop[ρ[[1]]-(2/3)ρ[[2]]+(1/5)ρ[[3]]]}
,{R,-2,-1,1/10}]

c=RecurrenceTable[{a[0]==1,a[1]==γ,a[k+2]==(γ/((k+2)((k+1)γ+1)))(a[k+1]+(((k(k+δ-1))/(4(R^2)))-Ε)a[k])}/.{γ->(1/3),δ->5},a,{k,0,10}];
c2=CoefficientList[((∑_(n=0)^(Length[c]-1))c[[n+1]](u^n))((∑_(m=0)^(Length[c]-1))c[[m+1]](u^m))/.Ε->-0.157412987479/.R->-5.3253,u]//N
NSolve[(∑_(n=0)^(Length[c2]-1))c2[[n+1]]((((2^(1+n)) (R^n) (n/2)!)/((1+(n/2))!))-(2/3)(-((3 (2^(1+n)) (R^n) ((((n/2)!))^2))/((-1+(n/2))! (2+(n/2))!)))+(1/5)((5 (2^(1+n)) (R^n) ((((n/2)!))^2))/((-2+(n/2))! (3+(n/2))!)))==0,R,Reals]

a=CoefficientList[1.` +0.6666467666189491` u+0.233711497250696` (u^2)+0.057892269462279114` (u^3)+0.011416060718727025` (u^4)+0.0019001512003665785` (u^5)+0.00027454051463773074` (u^6)+0.00003485368552570877` (u^7)+3.909550758995941`*^-6 (u^8)+3.886613919337475`*^-7 (u^9)+3.423932142765105`*^-8 (u^10)+2.6560194575496133`*^-9 (u^11)+1.7871244353862165`*^-10 (u^12)+1.017581390449739`*^-11 (u^13)+4.734306593585661`*^-13 (u^14)+1.714745687814212`*^-14 (u^15)+4.501929931069031`*^-16 (u^16)+7.574161874717616`*^-18 (u^17)+6.100436598427603`*^-20 (u^18),u]//N

Triplet P state

((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]-((8 π)/3)(Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]

(Subscript[ρ, Λ])[(Subscript[Ω, 1])]=2∫((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])
	=2∫(((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]-((8 π)/3)(Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])])\[DifferentialD](Subscript[Ω, 2])
	=2(((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]((32 (π^(3/2)))/(3 (Sqrt[5])))+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]((32 (π^(3/2)))/3))
	=((512 (π^(5/2)))/9)((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+(1/5)(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])

(Subscript[ρ, ψ])[(Subscript[Ω, 1])]=2∫((ψ[γ])^2)\[DifferentialD](Subscript[Ω, 2])=2(∑_(ℓ=0)^∞)(Subscript[c, ℓ])((4π)/(2ℓ+1))(∑_(m=-ℓ)^ℓ)(Subscript[Y, ℓ,m]^( * ))[(Subscript[Ω, 1])]∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]\[DifferentialD](Subscript[Ω, 2])

Φ[(Subscript[Ω, 1]),(Subscript[Ω, 2])]=Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])]ψ[γ]	((ψ[γ])^2)=(∑_(ℓ=0)^∞)(Subscript[c, ℓ])(Subscript[P, ℓ])[Cos[γ]]	(Subscript[c, ℓ])=((2ℓ+1)/2)(∫_0^π)(Subscript[P, ℓ])[Cos[γ]]((ψ[γ])^2)Sin[γ]\[DifferentialD]γ

ρ[(Subscript[Ω, 1])]=2∫((ψ[γ])^2)((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])=2(∑_(ℓ=0)^∞)(Subscript[c, ℓ])((4π)/(2ℓ+1))(∑_(m=-ℓ)^ℓ)(Subscript[Y, ℓ,m]^( * ))[(Subscript[Ω, 1])]∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])

∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])=∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])](((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[Ω, 1])](Subscript[Y, 0,0])[(Subscript[Ω, 2])]+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 0,0])[(Subscript[Ω, 2])]-((8 π)/3)(Subscript[Y, 1,0])[(Subscript[Ω, 1])](Subscript[Y, 1,0])[(Subscript[Ω, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 2,0])[(Subscript[Ω, 2])]+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 0,0])[(Subscript[Ω, 2])])\[DifferentialD](Subscript[Ω, 2])
					     =((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,0])(Subscript[δ, m,0])+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,0])(Subscript[δ, m,0])-((8 π)/3)(Subscript[Y, 1,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,1])(Subscript[δ, m,0])+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,2])(Subscript[δ, m,0])+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,0])(Subscript[δ, m,0])

Simplify[(((Cos[(Subscript[θ, 1])]-Cos[(Subscript[θ, 2])]))^2)==((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]-((8 π)/3)(Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]]
(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//Simplify

ρ[(Subscript[Ω, 1])]=((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 2,0])[(Subscript[Ω, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)-((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 2,0])[(Subscript[Ω, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)
	   =((32 (π^2))/3) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)+((Subscript[c, 0])-(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2))

Solve[(∫_0^π)(∫_0^(2π)) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Subscript[c, 0])-(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))Sin[(Subscript[θ, 1])]\[DifferentialD](Subscript[ϕ, 1])\[DifferentialD](Subscript[θ, 1])==2//Simplify,(Subscript[c, 0])]//Expand

{(Subscript[c, 0])-(1/5)(Subscript[c, 2]),(Subscript[c, 0])-(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2])}/.(Subscript[c, 0])->1+((Subscript[c, 1])/3)//Expand

∫ρ[(Subscript[Ω, 1])]\[DifferentialD](Subscript[Ω, 1])=2	\[Implies]	(Subscript[c, 0])=1+((Subscript[c, 1])/3)	\[Implies]	ρ[(Subscript[Ω, 1])]=((32 (π^2))/3) ((1+((Subscript[c, 1])/3)-((Subscript[c, 2])/5)) (((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)+(1-((Subscript[c, 1])/3)+((Subscript[c, 2])/5))(((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2))

2-(Subscript[c, 0])+(1/5)(Subscript[c, 2])=0	\[Implies]	(Subscript[c, 2])==5(Subscript[c, 0])-10

Solve[(Subscript[c, 0])-(1/5)(Subscript[c, 2])==0,(Subscript[c, 2])]

Solve[(Subscript[c, 0])-(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2])==0/.(Subscript[c, 2])->5 (Subscript[c, 0]),(Subscript[c, 1])]

(Subscript[ψ, n_])[u_]=(u^(n-1));
With[{M=10,d=2},
H[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u]((((u^2)/(4(R^2)))-1)(Subscript[ψ, j]^″)[u]+((((2d+1)u)/(4(R^2)))-((d+1)/u))  (Subscript[ψ, j]^′)[u]+((1/u)+(d/(2(R^2))))(Subscript[ψ, j])[u])((((u/R)))^(d-1))(((1-((((u/(2R))))^2)))^((d-2)/2))\[DifferentialD]u,{i,M},{j,M}];
S[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u](Subscript[ψ, j])[u]((((u/R)))^(d-1))(((1-((((u/(2R))))^2)))^((d-2)/2))\[DifferentialD]u,{i,M},{j,M}];
];

With[{M=10},
tab3P3=Table[
{ϵ,c}=SortEigensystem[Eigensystem[{N[H[R],100],N[S[R],100]}]];
f=Table[(Subscript[ψ, n])[2 R Sin[θ/2]],{n,M}];
ρ=Table[(ℓ+(1/2))NIntegrate[(Subscript[P, ℓ])[Cos[θ]](((c[[3]].f))^2)Sin[θ],{θ,0,π},WorkingPrecision->30],{ℓ,0,2}];
(*n[θ_]=(ρ[[1]]-((ρ[[3]])/5))(((Subscript[Y, 0,0])[θ,0])^2)+(ρ[[1]]-2((ρ[[2]])/3)+((ρ[[3]])/5))(((Subscript[Y, 1,0])[θ,0])^2);
Print[NIntegrate[4π n[θ]Sin[θ],{θ,0,π}]];
Print[NIntegrate[((((Subscript[Y, 1,0])[(Subscript[θ, 1]),0]-(Subscript[Y, 1,0])[(Subscript[θ, 2]),0]))^2)Sin[(Subscript[θ, 1])]Sin[(Subscript[θ, 2])],{(Subscript[θ, 1]),0,π},{(Subscript[θ, 2]),0,π}]];
Print[Plot[{((Subscript[Y, 0,0])[θ,0]+(1/5)(Subscript[Y, 2,0])[θ,0]),((Subscript[Y, 0,0])[θ,0])},{θ,0,π}]];*)
{R,ρ[[1]]-(2/3)ρ[[2]]+(1/5)ρ[[3]]},
{R,
-40,-5,(5/10)}]
];

Show[{
ListPlot[{tab3P1,tab3P2,tab3P3},ImageSize->500,Joined->True,InterpolationOrder->2,PlotRange->{All,{-0.00135,0.001}},PlotTheme->"Scientific",AxesStyle->Thick,BaseStyle->18,GridLines->Automatic,AxesStyle->Directive[Black,Thick],FrameStyle->Directive[Black,Thick],FrameLabel->{Style["μ",FontSize->18],Style["(Subscript[c, 0])-((2 (Subscript[c, 1]))/3)+((Subscript[c, 2])/5)",FontSize->18]},PlotLegends->Placed[{Style["(^3)P ground state",FontSize->18],Style["(^3)P first excited state",FontSize->18],Style["(^3)P second excited state",FontSize->18]},{Left,Bottom}]]
,
ListPlot[{{{-5.3257,0}},{{-15.9,0}},{{-30.9,0}}},PlotTheme->"Scientific",PlotMarkers->{Automatic,Medium}]
}]

ψc[[1]].f

((Subscript[c, 0])-((Subscript[c, 2])/5))(((Subscript[Y, 0,0])[θ,ϕ])^2)+((Subscript[c, 0])-2((Subscript[c, 1])/3)+((Subscript[c, 2])/5))(((Subscript[Y, 1,0])[θ,ϕ])^2)

Γ[1,2]=ρ[1]ρ[2]+(Subscript[Γ, xc])[1,2]

(ρ^hole)[1,2]=((Γ[1,2])/(ρ[1]))

Triplet P state at the HF level

With[{M=5},
H=Table[((Subscript[δ, μ,ν])/(R^2))((μ(μ+1))/2),{μ,0,M-1},{ν,0,M-1}];
ERI=Table[(1/R)(Sqrt[(2μ+1)(2ν+1)(2λ+1)(2σ+1)])(∑_(L=Max[Abs[μ-ν],Abs[λ-σ]])^(Min[μ+ν,λ+σ]))((J3[μ,ν,L])^2)((J3[λ,σ,L])^2),{μ,0,M-1},{ν,0,M-1},{λ,0,M-1},{σ,0,M-1}];
P=Table[(Subscript[c, μ])(Subscript[c, ν])+(Subscript[d, μ])(Subscript[d, ν]),{μ,0,M-1},{ν,0,M-1}];
EHF=(∑_(μ=1)^M)(∑_(ν=1)^M)P[[μ,ν]] H[[μ,ν]]+(1/2)(∑_(μ=1)^M)(∑_(ν=1)^M)(∑_(λ=1)^M)(∑_(σ=1)^M)P[[μ,ν]](ERI[[μ,ν,λ,σ]]-ERI[[μ,σ,λ,ν]])P[[λ,σ]];
]

With[{M=5},
t1=Table[{R,(R^2)((1/(R^2))+(2/(3R)))//Chop},{R,0.5,5,0.5}];
t2=Table[{R,(R^2)FindMinimum[{EHF,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[d, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ])(Subscript[d, ℓ])==0},Sequence[Flatten[Table[{(Subscript[c, ℓ]),(Subscript[d, ℓ])},{ℓ,0,M-1}]]],MaxIterations->1000][[1]]//Chop},{R,0.5,5,0.5}];
ListPlot[{t1,t2}]
]

With[{M=3},
Table[Integrate
[Evaluate[(∑_((Subscript[ℓ, 1])=0)^(M-1))(∑_((Subscript[ℓ, 2])=0)^(M-1))(Subscript[Y, (Subscript[ℓ, 1])])[θ]P[[(Subscript[ℓ, 1])+1,(Subscript[ℓ, 2])+1]](Subscript[Y, (Subscript[ℓ, 2])])[θ]2π Sin[θ]/.FindMinimum[{EHF,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[d, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ])(Subscript[d, ℓ])==0},Sequence[Flatten[Table[{(Subscript[c, ℓ]),(Subscript[d, ℓ])},{ℓ,0,M-1}]]],MaxIterations->1000][[2]]//Chop],{θ,0,π}],{R,-5,-3,1}]
]

With[{M=5},
Table[Plot[Evaluate[(∑_((Subscript[ℓ, 1])=0)^(M-1))(∑_((Subscript[ℓ, 2])=0)^(M-1))(Subscript[Y, (Subscript[ℓ, 1])])[θ]P[[(Subscript[ℓ, 1])+1,(Subscript[ℓ, 2])+1]](Subscript[Y, (Subscript[ℓ, 2])])[θ]2π Sin[θ]/.NMinimize[{EHF,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[d, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ])(Subscript[d, ℓ])==0},Sequence[Flatten[Table[{(Subscript[c, ℓ]),(Subscript[d, ℓ])},{ℓ,0,M-1}]]],MaxIterations->1000][[2]]//Chop],{θ,0,π}],{R,-10,-1,1}]
]

ρ=(∑_(μν))(Subscript[ϕ, μ])(Subscript[P, μν])(Subscript[ϕ, ν])
=(Subscript[P, 00])(Subscript[Y, 0]^2)+(Subscript[P, 11])(Subscript[Y, 1]^2)+(Subscript[P, 22])(Subscript[Y, 2]^2)+((Subscript[P, 01])+(Subscript[P, 10]))(Subscript[Y, 0])(Subscript[Y, 1])+((Subscript[P, 02])+(Subscript[P, 20]))(Subscript[Y, 0])(Subscript[Y, 2])+((Subscript[P, 12])+(Subscript[P, 21]))(Subscript[Y, 1])(Subscript[Y, 2])
=(Subscript[P, 00])(Subscript[Y, 0]^2)+(Subscript[P, 11])(Subscript[Y, 1]^2)+(Subscript[P, 22])(Subscript[Y, 2]^2)+((Subscript[P, 01])+(Subscript[P, 10]))(Subscript[Y, 0])(Subscript[Y, 1])-((Sqrt[5])/2)((Subscript[P, 02])+(Subscript[P, 20])) (((Subscript[Y, 0])^2)-((Subscript[Y, 1])^2))-2(Sqrt[5])π((Subscript[P, 12])+(Subscript[P, 21]))(Subscript[Y, 1])(Subscript[Y, 0])(((Subscript[Y, 0])^2)-((Subscript[Y, 1])^2))
=((Subscript[P, 00])-((Sqrt[5])/2)((Subscript[P, 02])+(Subscript[P, 20])))(Subscript[Y, 0]^2)+((Subscript[P, 11])+((Sqrt[5])/2)((Subscript[P, 02])+(Subscript[P, 20])) )(Subscript[Y, 1]^2)+(Subscript[P, 22])(Subscript[Y, 2]^2)+((Subscript[P, 01])+(Subscript[P, 10]))(Subscript[Y, 0])(Subscript[Y, 1])-2(Sqrt[5])π((Subscript[P, 12])+(Subscript[P, 21]))(Subscript[Y, 1])(Subscript[Y, 0])(((Subscript[Y, 0])^2)-((Subscript[Y, 1])^2))

With[{M=5},
tabHF=Table[(∑_((Subscript[ℓ, 1])=0)^(M-1))(∑_((Subscript[ℓ, 2])=0)^(M-1))(Subscript[Y, (Subscript[ℓ, 1])])[θ]P[[(Subscript[ℓ, 1])+1,(Subscript[ℓ, 2])+1]](Subscript[Y, (Subscript[ℓ, 2])])[θ]2π Sin[θ]/.FindMinimum[{EHF,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[d, ℓ]^2)==1,(∑_(ℓ=0)^(M-1))(Subscript[c, ℓ])(Subscript[d, ℓ])==0},Sequence[Flatten[Table[{(Subscript[c, ℓ]),(Subscript[d, ℓ])},{ℓ,0,M-1}]]],MaxIterations->1000][[2]]//Chop,{R,{-8,-7,-6}}];
]

Plot[tabHF,{θ,0,π},ImageSize->400,PlotRange->{All,All},PlotTheme->"Scientific",AxesStyle->Thick,BaseStyle->18,GridLines->Automatic,AxesStyle->Directive[Black,Thick],FrameStyle->Directive[Black,Thick],FrameTicks->{{Automatic,None},{{0,π/2,π},None}},PlotStyle->Thick,FrameLabel->{Style["θ",FontSize->18],Style["((OverTilde[ρ])^HF)(θ)",FontSize->16]},PlotLegends->Placed[{Style["μ = -8",FontSize->18],Style["μ = -7",FontSize->18],Style["μ = -6",FontSize->18]},{Center,Bottom}]]

Singlet P state

Φ[(Subscript[Ω, 1]),(Subscript[Ω, 2])]=Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])]ψ[γ]	((ψ[γ])^2)=(∑_(ℓ=0)^∞)(Subscript[c, ℓ])(Subscript[P, ℓ])[Cos[γ]]	(Subscript[c, ℓ])=((2ℓ+1)/2)(∫_0^π)(Subscript[P, ℓ])[Cos[γ]]((ψ[γ])^2)Sin[γ]\[DifferentialD]γ

ρ[(Subscript[Ω, 1])]=2∫((ψ[γ])^2)((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])=2(∑_(ℓ=0)^∞)(Subscript[c, ℓ])((4π)/(2ℓ+1))(∑_(m=-ℓ)^ℓ)(Subscript[Y, ℓ,m]^( * ))[(Subscript[Ω, 1])]∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])

∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])=∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])](((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[Ω, 1])](Subscript[Y, 0,0])[(Subscript[Ω, 2])]+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 0,0])[(Subscript[Ω, 2])]+((8 π)/3)(Subscript[Y, 1,0])[(Subscript[Ω, 1])](Subscript[Y, 1,0])[(Subscript[Ω, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 2,0])[(Subscript[Ω, 2])]+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 0,0])[(Subscript[Ω, 2])])\[DifferentialD](Subscript[Ω, 2])
					     =((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,0])(Subscript[δ, m,0])+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,0])(Subscript[δ, m,0])+((8 π)/3)(Subscript[Y, 1,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,1])(Subscript[δ, m,0])+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,2])(Subscript[δ, m,0])+((4 π)/3)(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[δ, ℓ,0])(Subscript[δ, m,0])

Simplify[(((Cos[(Subscript[θ, 1])]+Cos[(Subscript[θ, 2])]))^2)==((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]]
Simplify[(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))]

((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==((32 (π^2))/3) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//Simplify

((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==((32 (π^2))/3) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//Simplify

((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]==((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)//FullSimplify

((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)//FullSimplify

((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(-((Sqrt[5])/2) ((((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)//FullSimplify

8 π(((Cos[(Subscript[θ, 1])]+Cos[(Subscript[θ, 2])]))^2)==((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)//FullSimplify

((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==((32 (π^2))/3) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//FullSimplify

((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==((32 (π^2))/3) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//FullSimplify

ρ[(Subscript[Ω, 1])]=((64 (π^2))/(3 (Sqrt[5])))(Subscript[c, 0])(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 2,0])[(Subscript[Ω, 1])]+((32 (π^2))/3)(Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)+((64 (π^2))/9)(Subscript[c, 1])(((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2)+((64 (π^2))/(15 (Sqrt[5])))(Subscript[c, 2])(Subscript[Y, 0,0])[(Subscript[Ω, 1])](Subscript[Y, 2,0])[(Subscript[Ω, 1])]+((32 (π^2))/3)(Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)
	   =((32 (π^2))/3) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)+((Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2))

Solve[(∫_0^π)(∫_0^(2π)) (((Subscript[c, 0])-(1/5)(Subscript[c, 2])) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2]))(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))Sin[(Subscript[θ, 1])]\[DifferentialD](Subscript[ϕ, 1])\[DifferentialD](Subscript[θ, 1])==2,(Subscript[c, 0])]//Expand

{(Subscript[c, 0])-(1/5)(Subscript[c, 2]),(Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2])}/.(Subscript[c, 0])->1-((Subscript[c, 1])/3)

(Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2])/.(Subscript[c, 2])->-10+5 (Subscript[c, 0])//Simplify

Plot[(Subscript[P, 0])[Cos[θ]]Sin[θ],{θ,0,π}]

Solve[(2/3) (3 (Subscript[c, 0])+(Subscript[c, 1]))==2,(Subscript[c, 0])]//Expand

(Subscript[c, 0])+(2/3)(Subscript[c, 1])+(1/5)(Subscript[c, 2])/.(Subscript[c, 0])->1-((Subscript[c, 1])/3)

1+((Subscript[c, 1])/3)+((Subscript[c, 2])/5)/.(Subscript[c, 0])->1-((Subscript[c, 1])/3)

(Subscript[c, 0])-(1/5)(Subscript[c, 2])

∫ρ[(Subscript[Ω, 1])]\[DifferentialD](Subscript[Ω, 1])=2	\[Implies]	(Subscript[c, 0])=1-((Subscript[c, 1])/3)	\[Implies]	ρ[(Subscript[Ω, 1])]=((32 (π^2))/3) ((1+((Subscript[c, 1])/3)+((Subscript[c, 2])/5)) (((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)+(1+((Subscript[c, 1])/3)+((Subscript[c, 2])/5))(((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2))

2-(Subscript[c, 0])+(1/5)(Subscript[c, 2])=0	\[Implies]	(Subscript[c, 2])==5(Subscript[c, 0])-10

(Subscript[ψ, n_])[u_]=(u^(n-1));
With[{M=10,d=2},
H[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u]((((u^2)/(4(R^2)))-1)(Subscript[ψ, j]^″)[u]+((((2d+1)u)/(4(R^2)))-((d-1)/u))  (Subscript[ψ, j]^′)[u]+((1/u)+(d/(2(R^2))))(Subscript[ψ, j])[u])((((u/R)))^(d-1))(((1-((((u/(2R))))^2)))^((d-2)/2))\[DifferentialD]u,{i,M},{j,M}];
S[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u](Subscript[ψ, j])[u]((((u/R)))^(d-1))(((1-((((u/(2R))))^2)))^((d-2)/2))\[DifferentialD]u,{i,M},{j,M}];
];

With[{M=10},
tab=Table[
{ϵ,c}=SortEigensystem[Eigensystem[{N[H[R],200],N[S[R],200]}]];
f=Table[(Subscript[ψ, n])[2 R Sin[θ/2]],{n,M}];
ρ=Table[(ℓ+(1/2))NIntegrate[(Subscript[P, ℓ])[Cos[θ]](((c[[2]].f))^2)Sin[θ],{θ,0,π},WorkingPrecision->50],{ℓ,0,2}];
(*Print["R: ",N[R]];
Print["ρ: ",ρ];*)
{N[R],Chop[ρ[[1]]+(2/3)ρ[[2]]-(1/5)ρ[[3]]]}
,{R,-20+1/10,-5+10/10,1}]
]
ListPlot[tab,Joined->True]

Triplet (p^2)

(((Sin[(Subscript[θ, 1])]Sin[(Subscript[θ, 2])]Cos[(Subscript[ϕ, 1])-(Subscript[ϕ, 2])]))^2)==(((((4 π)/3)))^2)((( (Subscript[Y, 1,1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,-1])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(Subscript[Y, 1,-1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 1,1])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]))^2)//FullSimplify
(Subscript[Y, 1,1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,-1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==-(1/(2 (Sqrt[π])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+(1/(2 (Sqrt[5 π])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]//Simplify
(Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==(1/96) (160 (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-80 (((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+16 (((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//Simplify
{(((Subscript[Y, 1,-1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==(Sqrt[(3/(10 π))]) (Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])],(((Subscript[Y, 1,1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==(Sqrt[(3/(10 π))]) (Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]}//Simplify
(((Sin[(Subscript[θ, 1])]Sin[(Subscript[θ, 2])]Cos[(Subscript[ϕ, 1])-(Subscript[ϕ, 2])]))^2)==((8π)/9) ((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]-(1/( (Sqrt[5])))  ((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])])+(9/45)  (Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(9/15) ((Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,-2])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]))//FullSimplify

((8π)/9) ((Subscript[c, 0])(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(1/( (Sqrt[5])))  ((Subscript[c, 0])(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((Subscript[c, 2])/5)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])+(9/45)((Subscript[c, 2])/5)(((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(9/15)((Subscript[c, 2])/5)(((-1))^2) ((((Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(((Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))==((8π)/9) ( (3/2)((Subscript[c, 0])+((Subscript[c, 2])/15))(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)- (1/2)((Subscript[c, 0])+((Subscript[c, 2])/5))(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(9/45)((Subscript[c, 2])/5)(((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(9/15)((Subscript[c, 2])/5)((((Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(((Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))//Simplify

(^3)(D^o)

((((Cos[(Subscript[θ, 1])]-Cos[(Subscript[θ, 2])])Sin[(Subscript[θ, 1])]Sin[(Subscript[θ, 2])]Cos[(Subscript[ϕ, 1])-(Subscript[ϕ, 2])]))^2)==(((((4 π)/3)))^2)(((Cos[(Subscript[θ, 1])]-Cos[(Subscript[θ, 2])]))^2)((( (Subscript[Y, 1,1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,-1])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(Subscript[Y, 1,-1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 1,1])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]))^2)//FullSimplify
(Subscript[Y, 1,1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,-1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==-(1/(2 (Sqrt[π])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+(1/(2 (Sqrt[5 π])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]//Simplify
(Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==(1/96) (160 (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-80 (((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+16 (((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//Simplify
{(((Subscript[Y, 1,-1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==(Sqrt[(3/(10 π))]) (Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])],(((Subscript[Y, 1,1])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==(Sqrt[(3/(10 π))]) (Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]}//Simplify
((((Cos[(Subscript[θ, 1])]-Cos[(Subscript[θ, 2])])Sin[(Subscript[θ, 1])]Sin[(Subscript[θ, 2])]Cos[(Subscript[ϕ, 1])-(Subscript[ϕ, 2])]))^2)==(((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]-((8 π)/3)(Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])])((8/9) π (Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]-((8 π (Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])] (Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])/(9 (Sqrt[5])))-((8 π (Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])])/(9 (Sqrt[5])))+(8/45) π (Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(8/15) π (Subscript[Y, 2,-2])[(Subscript[θ, 2]),(Subscript[ϕ, 2])] (Subscript[Y, 2,2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+(8/15) π (Subscript[Y, 2,-2])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 2]),(Subscript[ϕ, 2])])//FullSimplify

(((8 π)/(3 (Sqrt[5])))(Subscript[Y, 2,0])[(Subscript[θ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2])]+((8 π)/3)(Subscript[Y, 0,0])[(Subscript[θ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2])]-((8 π)/3)(Subscript[Y, 1,0])[(Subscript[θ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2])]+((8 π)/(3 (Sqrt[5])))(Subscript[Y, 0,0])[(Subscript[θ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 2])])((8/9) π (Subscript[Y, 0,0])[(Subscript[θ, 1])] (Subscript[Y, 0,0])[(Subscript[θ, 2])]-((8 π (Subscript[Y, 0,0])[(Subscript[θ, 2])] (Subscript[Y, 2,0])[(Subscript[θ, 1])])/(9 (Sqrt[5])))-((8 π (Subscript[Y, 0,0])[(Subscript[θ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2])])/(9 (Sqrt[5])))+(8/45) π (Subscript[Y, 2,0])[(Subscript[θ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2])]+(8/15) π (Subscript[Y, 2,-2])[(Subscript[θ, 2])] (Subscript[Y, 2,2])[(Subscript[θ, 1])]+(8/15) π (Subscript[Y, 2,-2])[(Subscript[θ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 2])])==(64/27) (π^2) (1/(Sqrt[4π]))(Subscript[Y, 0,0])[(Subscript[θ, 1])] (((Subscript[Y, 0,0])[(Subscript[θ, 2])])^2)-(64/27) (π^2) (1/(Sqrt[4π])) (1/(Sqrt[4π]))(Subscript[Y, 1,0])[(Subscript[θ, 1])] (Subscript[Y, 1,0])[(Subscript[θ, 2])]+((64 (π^2)(1/(Sqrt[4π])) (Subscript[Y, 1,0])[(Subscript[θ, 1])]  (Subscript[Y, 2,0])[(Subscript[θ, 1])](Subscript[Y, 1,0])[(Subscript[θ, 2])])/(27 (Sqrt[5])))-(64/135) (π^2) (1/(Sqrt[4π]))(Subscript[Y, 0,0])[(Subscript[θ, 2])] (((Subscript[Y, 2,0])[(Subscript[θ, 1])])^2)+((64 (π^2) (1/(Sqrt[4π])) (Subscript[Y, 1,0])[(Subscript[θ, 1])] (Subscript[Y, 1,0])[(Subscript[θ, 2])] (Subscript[Y, 2,0])[(Subscript[θ, 2])])/(27 (Sqrt[5])))-(64/135) (π^2) (1/(Sqrt[4π]))(Subscript[Y, 0,0])[(Subscript[θ, 2])] (Subscript[Y, 2,0])[(Subscript[θ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2])]-(64/135) (π^2) (Subscript[Y, 1,0])[(Subscript[θ, 1])] (Subscript[Y, 1,0])[(Subscript[θ, 2])] (Subscript[Y, 2,0])[(Subscript[θ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2])]+((64 (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 2])] (((Subscript[Y, 2,0])[(Subscript[θ, 1])])^2) (Subscript[Y, 2,0])[(Subscript[θ, 2])])/(135 (Sqrt[5])))-(64/135) (π^2) (((Subscript[Y, 0,0])[(Subscript[θ, 1])])^2) (((Subscript[Y, 2,0])[(Subscript[θ, 2])])^2)+((64 (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 1])] (((Subscript[Y, 2,0])[(Subscript[θ, 2])])^2))/(135 (Sqrt[5])))+(64/45) (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 1])] (Subscript[Y, 0,0])[(Subscript[θ, 2])] (Subscript[Y, 2,-2])[(Subscript[θ, 2])] (Subscript[Y, 2,2])[(Subscript[θ, 1])]-(64/45) (π^2) (Subscript[Y, 1,0])[(Subscript[θ, 1])] (Subscript[Y, 1,0])[(Subscript[θ, 2])] (Subscript[Y, 2,-2])[(Subscript[θ, 2])] (Subscript[Y, 2,2])[(Subscript[θ, 1])]+((64 (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 2])] (Subscript[Y, 2,-2])[(Subscript[θ, 2])] (Subscript[Y, 2,0])[(Subscript[θ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 1])])/(45 (Sqrt[5])))+((64 (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 1])] (Subscript[Y, 2,-2])[(Subscript[θ, 2])] (Subscript[Y, 2,0])[(Subscript[θ, 2])] (Subscript[Y, 2,2])[(Subscript[θ, 1])])/(45 (Sqrt[5])))+(64/45) (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 1])] (Subscript[Y, 0,0])[(Subscript[θ, 2])] (Subscript[Y, 2,-2])[(Subscript[θ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 2])]-(64/45) (π^2) (Subscript[Y, 1,0])[(Subscript[θ, 1])] (Subscript[Y, 1,0])[(Subscript[θ, 2])] (Subscript[Y, 2,-2])[(Subscript[θ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 2])]+((64 (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 2])] (Subscript[Y, 2,-2])[(Subscript[θ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 1])] (Subscript[Y, 2,2])[(Subscript[θ, 2])])/(45 (Sqrt[5])))+((64 (π^2) (Subscript[Y, 0,0])[(Subscript[θ, 1])] (Subscript[Y, 2,-2])[(Subscript[θ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2])] (Subscript[Y, 2,2])[(Subscript[θ, 2])])/(45 (Sqrt[5])))//FullSimplify

Triplet D state

Φ[(Subscript[Ω, 1]),(Subscript[Ω, 2])]=Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])]ψ[γ]	((ψ[γ])^2)=(∑_(ℓ=0)^∞)(Subscript[c, ℓ])(Subscript[P, ℓ])[Cos[γ]]	(Subscript[c, ℓ])=((2ℓ+1)/2)(∫_0^π)(Subscript[P, ℓ])[Cos[γ]]((ψ[γ])^2)Sin[γ]\[DifferentialD]γ

ρ[(Subscript[Ω, 1])]=2∫((ψ[γ])^2)((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])=8π(∑_(ℓ=0)^∞)((Subscript[c, ℓ])/(2ℓ+1))(∑_(m=-ℓ)^ℓ)(Subscript[Y, ℓ,m]^( * ))[(Subscript[Ω, 1])]∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])

∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]((Λ[(Subscript[Ω, 1]),(Subscript[Ω, 2])])^2)\[DifferentialD](Subscript[Ω, 2])=∫(Subscript[Y, ℓ,m])[(Subscript[Ω, 2])]((Sqrt[4π])(((Subscript[Y, 2,0])[(Subscript[Ω, 1])])^2)(Subscript[Y, 0,0])[(Subscript[Ω, 2])]-2 (Subscript[Y, 2,0])[(Subscript[Ω, 1])] (Subscript[Y, 2,0])[(Subscript[Ω, 2])]+(((Subscript[Y, 2,0])[(Subscript[Ω, 2])])^2))\[DifferentialD](Subscript[Ω, 2])
					     =(Sqrt[4π])(Subscript[c, 0])(((Subscript[Y, 2,0])[(Subscript[Ω, 1])])^2)+2 (Sqrt[5 π]) ((Subscript[c, 2])/5)((((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)-(((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2))+ 4π ((1/(Sqrt[4 π]))(Subscript[c, 0])+(1/7) (Sqrt[(5/π)]) ((Subscript[c, 2])/5)+((3 )/(7 (Sqrt[π])))((Subscript[c, 4])/9))(((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)
					      =(2 (Sqrt[π]) (Subscript[c, 0])+(18/7) (Sqrt[(π/5)]) (Subscript[c, 2])+((4 (Sqrt[π]) (Subscript[c, 4]))/21)) (((Subscript[Y, 0,0])[(Subscript[Ω, 1])])^2)-2 (Sqrt[(π/5)]) (Subscript[c, 2]) (((Subscript[Y, 1,0])[(Subscript[Ω, 1])])^2)+2 (Sqrt[π]) (Subscript[c, 0]) (((Subscript[Y, 2,0])[(Subscript[Ω, 1])])^2)

Simplify[(((((Cos[(Subscript[θ, 1])])^2)-((Cos[(Subscript[θ, 2])])^2)))^2)== ((16π)/45)((((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]-(Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 1])]))^2)]
(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==((Sqrt[5])/2)((((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//Simplify
(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 4,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]==(1/(Sqrt[4π]))(((3 )/(16 (Sqrt[π]))))(12π (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-30 ((4π)/3)(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(35/9)(((16 π)/5)(((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+8(Sqrt[(π/5)])(Sqrt[4π])((Sqrt[5])/2)((((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)-(((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))+4π (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)))//Simplify
(((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==(1/(Sqrt[4π]))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+(1/(Sqrt[5 π]))(Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]//Simplify
(((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)==(1/(Sqrt[4 π]))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+(1/7) (Sqrt[(5/π)]) (Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]+((3 )/(7 (Sqrt[π])))(Subscript[Y, 4,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]//Simplify
(((((Cos[(Subscript[θ, 1])])^2)-((Cos[(Subscript[θ, 2])])^2)))^2)==((32 π)/45)( (Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+((Sqrt[5])/7) ((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])](Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])])- (Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+(3/7)((Subscript[Y, 4,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 0,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]+ (Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 4,0])[(Subscript[θ, 2]),(Subscript[ϕ, 2])]))//Simplify

((32 π)/45)((Subscript[c, 0]) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Sqrt[5])/7) ((Subscript[c, 0])+((Subscript[c, 2])/5))(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])]-((Subscript[c, 2])/5) (((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(3/7)((Subscript[c, 0])+((Subscript[c, 4])/9) )(Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])] (Subscript[Y, 4,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])==((16 π)/45)(((Subscript[c, 0])-((Subscript[c, 2])/7)-((2 (Subscript[c, 4]))/63)) (((Subscript[Y, 0,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+(1/7)((Subscript[c, 2])-((5 (Subscript[c, 4]))/9)) (((Subscript[Y, 1,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2)+((Subscript[c, 0])-((2 (Subscript[c, 2]))/5)+((Subscript[c, 4])/9)) (((Subscript[Y, 2,0])[(Subscript[θ, 1]),(Subscript[ϕ, 1])])^2))//Simplify

(Subscript[ψ, n_])[u_]=(u^(n-1));
With[{M=10},
H[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u]((((u^2)/(4(R^2)))-1)(Subscript[ψ, j]^″)[u]+(((7u)/(4(R^2)))-(3/u))  (Subscript[ψ, j]^′)[u]+((1/u)+(3/(R^2)))(Subscript[ψ, j])[u])u\[DifferentialD]u,{i,M},{j,M}];
S[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u](Subscript[ψ, j])[u]u\[DifferentialD]u,{i,M},{j,M}];
];

With[{M=10},
tab=Table[
{ϵ,c}=SortEigensystem[Eigensystem[{N[H[R],100],N[S[R],100]}]];
f=Table[(Subscript[ψ, n])[2 R Sin[θ/2]],{n,M}];
ρ=Table[(ℓ+(1/2))NIntegrate[(Subscript[P, ℓ])[Cos[θ]](((c[[1]].f))^2)Sin[θ],{θ,0,π},WorkingPrecision->15],{ℓ,0,4}];
ρ=Normalize[ρ];
Chop[{N[R],ρ[[3]]-(5/9) ρ[[5]],ρ[[1]]-(2/5)ρ[[3]]+(1/9)ρ[[1]]}]
,{R,-20,20}];
]
ListPlot[{tab[[;;,{1,2}]],tab[[;;,{1,3}]]},Joined->True]