Initialisation

(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}"}];

Spherium: RHF, UHF and symmetry breaking

The general function for the Hartree-Fock energy is:

EHF[R_,λ_]=(∑_(L=0)^1)(Subscript[a, L]^2)((L(L+1))/(2(R^2)))+(∑_(L=0)^1)(Subscript[b, L]^2)((L(L+1))/(2(R^2)))+(∑_(i=0)^1)(∑_(j=0)^1)(∑_(k=0)^1)(∑_(l=0)^1)(Subscript[a, i])(Subscript[a, j])(Subscript[b, k])(Subscript[b, 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[a, 0])->Cos[x],(Subscript[a, 1])->Sin[x],(Subscript[b, 0])->Cos[y],(Subscript[b, 1])->Sin[y]};

We search stationnary solution

sol=Solve[{(∂_x)EHF[R,1]==0,(∂_y)EHF[R,1]==0},{x,y}]/.{C[1]->0,C[2]->0}//FullSimplify

Some solutions correspond to RHF solutions et some to UHF solutions.

solRHF={{x->0,y->0},{x->0,y->π},{x->-(π/2),y->-(π/2)},{x->-(π/2),y->(π/2)},{x->(π/2),y->-(π/2)},{x->(π/2),y->(π/2)},{x->π,y->0},{x->π,y->π},{x->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))]},{x->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))]},{x->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))]},{x->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),((5 (Sqrt[3+2 R]))/(Sqrt[R]))]}};

solUHF={{x->0,y->-(π/2)},{x->0,y->(π/2)},{x->-(π/2),y->0},{x->-(π/2),y->π},{x->(π/2),y->0},{x->(π/2),y->π},{x->π,y->-(π/2)},{x->π,y->(π/2)},{x->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]},{x->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]},{x->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]},{x->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[-((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]},{x->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]}};

The two electrons integrals are equal to

(Subscript[U, i_,j_,k_,l_]):=(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);

The Fock operator for the electron α is

Fα=Table[KroneckerDelta[i,j]((i(i+1))/(2(R^2)))+(λ/R)(∑_(k=0)^1)(∑_(l=0)^1)(Subscript[a, k])(Subscript[a, l])((Subscript[U, i,j,k,l])-(Subscript[U, i,l,k,j]))+(λ/R)(∑_(k=0)^1)(∑_(l=0)^1)(Subscript[b, k])(Subscript[b, l])((Subscript[U, i,j,k,l]))/.{(Subscript[a, 0])->Cos[x],(Subscript[a, 1])->Sin[x],(Subscript[b, 0])->Cos[y],(Subscript[b, 1])->Sin[y]}/.λ->1//FullSimplify//Expand,{i,0,1},{j,0,1}]

And the Fock operator β

Fβ=Table[KroneckerDelta[i,j]((i(i+1))/(2(R^2)))+(λ/R)(∑_(k=0)^1)(∑_(l=0)^1)(Subscript[b, k])(Subscript[b, l])((Subscript[U, i,j,k,l])-(Subscript[U, i,l,k,j]))+(λ/R)(∑_(k=0)^1)(∑_(l=0)^1)(Subscript[a, k])(Subscript[a, l])((Subscript[U, i,j,k,l]))/.{(Subscript[a, 0])->Cos[x],(Subscript[a, 1])->Sin[x],(Subscript[b, 0])->Cos[y],(Subscript[b, 1])->Sin[y]}/.λ->1//FullSimplify//Expand,{i,0,1},{j,0,1}]

RHF solutions

In the RHF case Fα = Fβ

Fα/.solRHF//FullSimplify//Expand

Fβ/.solRHF//FullSimplify//Expand

The 16 stationnary RHF solutions yields 3  differents Fockian which are all diagonal.

{{{(1/R),0},{0,(1/(R^2))+(5/(3 R))}}//MatrixForm,{{(5/(3 R)),0},{0,(1/(R^2))+(29/(25 R))}}//MatrixForm,{{(25/(44 (R^2)))+(91/(66 R)),0},{0,(25/(44 (R^2)))+(91/(66 R))}}//MatrixForm}

The first solution corresponds to a density matrix spatially symmetric where the two electrons are in (Subscript[Y, 0]) and the second to a density matrix spatially symmetric where the two electrons are in (Subscript[Y, 1]). The last solution corresponds to the symmetry broken RHF solution. In the sb-RHF solution the two orbitals are degenerate.

{(∂_(x,x))EHF[R,1],(∂_(y,y))EHF[R,1]}/.{{x->0,y->π},{x->(π/2),y->-(π/2)},{x->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[-75+38 R])/(Sqrt[R])),-((5 (Sqrt[3+2 R]))/(Sqrt[R]))]}}//FullSimplify

Solve[-((50+8 R)/(25 (R^2)))==0,R]

Plot[{(1/R),(5/(3R)),(1/(R^2))+(5/(3 R)),(1/(R^2))+(29/(25 R)),(25/(44 (R^2)))+(91/(66 R))},{R,-3,3},PlotLegends->{"RHF_GS_min","RHF_GS_max","RHF_ES_min","RHF_ES_max","sb_RHF"}]

Plot[{(1/R),(5/(3R)),(1/(R^2))+(5/(3 R)),(1/(R^2))+(29/(25 R)),(25/(44 (R^2)))+(91/(66 R))},{R,1.5,2.5},PlotLegends->{"RHF_GS_min","RHF_GS_max","RHF_ES_min","RHF_ES_max","sb_RHF"}]

Solve[(5/(3 R))==(1/(R^2))+(29/(25 R)),R]

Plot[{(1/R),(5/(3R)),(1/(R^2))+(5/(3 R)),(1/(R^2))+(29/(25 R)),(25/(44 (R^2)))+(91/(66 R))},{R,0,-3},PlotLegends->{"RHF_GS_min","RHF_GS_max","RHF_ES_min","RHF_ES_max","sb_RHF"}]

Solve[(1/R)==(1/(R^2))+(5/(3 R)),R]

The sb-RHF energy interacts with the minimized RHF for R<0 and the maximized RHF for R>0. For negative value of R there is a point at which the attraction of the sb-RHF will minimize the energy than the kinetical energy.  For positive value of lambda there is a point at which the repulsion of the sb RHF solution is bigger than the kinetical energy.

UHF solutions

Fα/.solUHF//FullSimplify//Expand

Fβ/.solUHF//FullSimplify//Expand

The 16 stationnary UHF solutions yields 3  differents Fockian, only two of them are already diagonal.

{{{(5/(3 R)),0},{0,(1/(R^2))+(1/R)}}//MatrixForm,{{(1/R),0},{0,(1/(R^2))+(137/(75 R))}}//MatrixForm,DiagonalMatrix[{{-(25/(56 (R^2)))+(109/(84 R)),((5 (Sqrt[-3+2 R]) (Sqrt[75+62 R]))/(84 (R^2)))},{((5 (Sqrt[-3+2 R]) (Sqrt[75+62 R]))/(84 (R^2))),(75/(56 (R^2)))+(121/(84 R))}}//Eigenvalues//Expand]//MatrixForm}

{(∂_(x,x))EHF[R,1],(∂_(y,y))EHF[R,1]}/.{{x->0,y->(π/2)},{x->-(π/2),y->0},{x->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),((5 (Sqrt[-3+2 R]))/(Sqrt[R]))],y->ArcTan[((Sqrt[75+62 R])/(Sqrt[R])),-((5 (Sqrt[-3+2 R]))/(Sqrt[R]))]}}//FullSimplify

Plot[{(1/R),(5/(3R)),(1/(R^2))+(1/R),(1/(R^2))+(137/(75 R)),(25/(56 (R^2)))+(59/(84 R)),(25/(56 (R^2)))+(57/(28 R))},{R,-3,3},PlotLegends->{"UHF_GS_min","UHF_GS_max","UHF_ES_min","UHF_ES_max","sb_UHF1","sb_UHF2"}]

Plot[{(1/R),(5/(3R)),(1/(R^2))+(1/R),(1/(R^2))+(137/(75 R)),(25/(56 (R^2)))+(59/(84 R)),(25/(56 (R^2)))+(57/(28 R))},{R,1,3},PlotLegends->{"UHF_GS_min","UHF_GS_max","UHF_ES_min","UHF_ES_max","sb_UHF1","sb_UHF2"}]

Plot[{(1/R),(5/(3R)),(1/(R^2))+(1/R),(1/(R^2))+(137/(75 R)),(25/(56 (R^2)))+(59/(84 R)),(25/(56 (R^2)))+(57/(28 R))},{R,0,-3},PlotLegends->{"UHF_GS_min","UHF_GS_max","UHF_ES_min","UHF_ES_max","sb_UHF1","sb_UHF2"}]

For R>0 the sb-UHF solution is a minimum because the symmetry breaking minimize the repulsion and for R<0 this is a maximum because the symmetry breaking minimize the attraction.
For R>0 the sb-UHF 1 is the “bonding” solution and the sb-UHF 2 is the antibonding solution. There is a point a which the bonding solution becomes better than the ground state of the symmetric solution, and an other point at which the anti bonding solution becomes higher than the symetric solution.
We can do the same reasoning for R<0.

sb-RHF and sb-UHF

Table[ EHF[R,1]/.sol[[k]],{k,Length[sol]}]//FullSimplify//Expand//DeleteDuplicates

Plot[{(1/R),(1/(R^2))+(1/R),(2/(R^2))+(29/(25 R)),(75/(88 (R^3)))+(25/(22 (R^2)))+(91/(66 R)),-(75/(112 (R^3)))+(25/(28 (R^2)))+(59/(84 R))},{R,1.5,3},PlotLegends->{"(s^2)_RHF","((Subscript[sp, z])_UHF) ???","((Subscript[p, z])^2)_RHF","sb_RHF/h_RHF","sb_UHF/h_UHF"}]

Plot[{(1/R),(1/(R^2))+(1/R),(2/(R^2))+(29/(25 R)),(75/(88 (R^3)))+(25/(22 (R^2)))+(91/(66 R)),-(75/(112 (R^3)))+(25/(28 (R^2)))+(59/(84 R))},{R,-1,-3},PlotLegends->{"(s^2)_RHF","((Subscript[sp, z])_UHF) ???","((Subscript[p, z])^2)_RHF","sb_RHF/h_RHF","sb_UHF/h_UHF"}]

For R > 0:
-the sb_RHF solution is a maximum and exist for R>75/38 (~1.97) before it is a holomorphic solution
-the sb_UHF solution is a minimum and exist for R>3/2  before it is a holomorphic solution
For R<0:
-the sb_RHF solution is a minimum and exist for R<-3/2 after it is a holomorphic solution
-the sb_UHF solution is a maximum and exist for R<-75/62(~1.2) after it is a holomorphic solution

Those points are the Coulson-Fischer points of the sb-RHF and sb-UHF solutions.