Hubbard Model in Complex Plane

Contributors: Hugh G. A. Burton
Last Updated: 14-Nov-2020 22:03:00

We develop the Hubbard dimer as a consistent example for the convergence behaviour of MP theories.

1) Singlet FCI on real axis and complex plane
2) RMP expansion
3) UMP expansion

Initialisation

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

(* Core Hamiltonian *)
Hc=(({{-δϵ, -t}, {-t, δϵ}}));
(* Two-electron integrals *)
G=Table[If[(μ==ν&&μ==λ&&μ==σ),U,0],{μ,2},{ν,2},{λ,2},{σ,2}];
c[θ_]={{Cos[θ],-Sin[θ]},{Sin[θ],Cos[θ]}};
P[θ_,w_,x_]:=Table[c[θ][[i,w]]c[θ][[j,x]],{i,2},{j,2}]
(* Core Hamiltonian and Fock matrices in MO basis *)
Fmo[θ1_,θ2_]= Transpose[c[θ1]].(Hc+Table[(∑_(λ=1)^2)(∑_(σ=1)^2)((P[θ1,1,1][[λ,σ]]+P[θ2,1,1][[λ,σ]])G[[μ,ν,λ,σ]]-P[θ1,1,1][[λ,σ]]G[[μ,σ,λ,ν]]),{μ,2},{ν,2}]).c[θ1];
Hmo[θ_]= Transpose[c[θ]].Hc.c[θ];
(* HF energy *)
Ehf[θα_,θβ_]:=(1/2)(Hmo[θα][[1,1]]+Hmo[θβ][[1,1]]+Fmo[θα,θβ][[1,1]]+Fmo[θβ,θα][[1,1]])

(*Equations to evaluate matrix elements *)
SlaterCondon[θα_,θβ_,iOcc_,jOcc_]:=Module[{iα,jα,iβ,jβ},
{iα,iβ}=iOcc;
{jα,jβ}=jOcc;
If[(iα==jα&&iβ==jβ),(∑_(μ=1)^2)(∑_(ν=1)^2)Hc[[μ,ν]](P[θα,iα,jα][[ν,μ]]+P[θβ,iβ,jβ][[ν,μ]])+(1/2)(∑_(μ=1)^2)(∑_(ν=1)^2)(∑_(σ=1)^2)(∑_(τ=1)^2)((P[θα,iα,jα][[ν,μ]]+P[θβ,iβ,jβ][[ν,μ]])(P[θα,iα,jα][[σ,τ]]+P[θβ,iβ,jβ][[σ,τ]])-(P[θα,iα,jα][[τ,μ]]P[θα,iα,jα][[ν,σ]]+P[θβ,iβ,jβ][[τ,μ]]P[θβ,iβ,jβ][[ν,σ]]))G[[μ,ν,σ,τ]],
If[(iα≠jα&&iβ==jβ),(∑_(μ=1)^2)(∑_(ν=1)^2)Hc[[μ,ν]]P[θα,iα,jα][[ν,μ]]+(∑_(μ=1)^2)(∑_(ν=1)^2)(∑_(σ=1)^2)(∑_(τ=1)^2)P[θα,iα,jα][[ν,μ]]P[θβ,iβ,jβ][[σ,τ]]G[[μ,ν,σ,τ]],
If[(iα== jα&&iβ≠jβ),(∑_(μ=1)^2)(∑_(ν=1)^2)Hc[[μ,ν]]P[θβ,iβ,jβ][[ν,μ]]+(∑_(μ=1)^2)(∑_(ν=1)^2)(∑_(σ=1)^2)(∑_(τ=1)^2)P[θα,iα,jα][[ν,μ]]P[θβ,iβ,jβ][[σ,τ]]G[[μ,ν,σ,τ]],
If[(iα≠jα&&iβ≠jβ),(∑_(μ=1)^2)(∑_(ν=1)^2)(∑_(σ=1)^2)(∑_(τ=1)^2)P[θα,iα,jα][[ν,μ]]P[θβ,iβ,jβ][[σ,τ]]G[[μ,ν,σ,τ]],
0]
]]]
]
(* Occupation numbers for Hilbert space *)
Occ={{1,1},{1,2},{2,1},{2,2}};
(Subscript[H, FCI])[θα_,θβ_]=Table[SlaterCondon[θα,θβ,Occ[[i]],Occ[[j]]],{i,1,4},{j,1,4}]//FullSimplify;
(Subscript[H, sFCI])[θα_,θβ_]=(({{1, 0, 0, 0}, {0, (1/(Sqrt[2])), (1/(Sqrt[2])), 0}, {0, 0, 0, 1}})).(Subscript[H, FCI])[θα,θβ].(({{1, 0, 0}, {0, (1/(Sqrt[2])), 0}, {0, (1/(Sqrt[2])), 0}, {0, 0, 1}}))//FullSimplify;

(* MP reference Hamiltonian *)
H0[θα_,θβ_]=DiagonalMatrix[Table[Fmo[θα,θβ][[Occ[[i]][[1]],Occ[[i]][[1]]]]+Fmo[θβ,θα][[Occ[[i]][[2]],Occ[[i]][[2]]]],{i,1,4}]]//FullSimplify;
Hmp[θα_,θβ_]=H0[θα,θβ]+λ((Subscript[H, FCI])[θα,θβ]-H0[θα,θβ])//FullSimplify;

Complex FCI

(Subscript[H, sFCI])[0,0]/.{δϵ->0}//MatrixForm

First construct symmetric dimer FCI Hamiltonian in the conventional RHF-based (θ = π/4) singlet CSFs.

(Subscript[H, FCI])[0,0]/.{δϵ->0}//MatrixForm

Here we see a block form, so we only get coupling between the ground-state and doubly-excited state. 
In this minimal basis, the open-shell singlet is given exactly by the corresponding CSF.

We can therefore apply the 2x2 Hamiltonian example to solve for the ground-state and doubly-excited state. 
This gives analytic energies as...

Eigenvalues[(({{U, -t, -t, 0}, {-t, 0, 0, -t}, {-t, 0, 0, -t}, {0, -t, -t, U}}))]

We can analytically find the exceptional point and find the corresponding energy (it’s pure imaginary!)

Solve[(Sqrt[16 (λ^2)(t^2)+(U^2)])==0]

(1/2) (U-(Sqrt[16 (t^2)+(U^2)]))/.{{U->-4 I t},{U->4 I t}}

Plotting these energies as a function of U/t shows the avoided crossing on the real axis, 
and the exceptional point structure in the complex plane.

With[{U=4,t=1},Plot[{(1/2) (U-(Sqrt[16 (λ^(2 ))(t^2)+(U^2)])),0,U,(1/2) (U+(Sqrt[16 (λ^2)(t^2)+(U^2)]))},{λ,-5,5},
PlotRange->{-7,12},PlotRangePadding->0,FrameStyle->Directive[Thick,18,Black],AspectRatio->0.7,
PlotStyle->{Orange,Darker[Green],Red,Blue},
PlotLegends->Placed[{MaTeX["E_{-}",Magnification->1.5],MaTeX["E_{\\text{T}}",Magnification->1.5],MaTeX["E_{\\text{S}}",Magnification->1.5],MaTeX["E_{\\text{+}}",Magnification->1.5]},{0.5,0.8}],
FrameLabel->{MaTeX["\\lambda",Magnification->1.5],MaTeX["E/t",Magnification->1.5]},
ImageSize->Large,PlotTheme->"Scientific"]]

With[{t=1,U=4},
Show[Plot3D[{Re[(1/2) (U-(Sqrt[16 (((x+I y))^(2 ))(t^2)+(U^2)]))],Re[(1/2) (U+(Sqrt[16(((x+I y))^2)(t^2)+(U^2)]))]},{x,-2,2},{y,0,3},
BoxRatios->{1, 1, 1},ImageSize->Large,PlotStyle->Opacity[0.7],Mesh->5,ViewPoint->{1.3`,-1.6,0.7},Boxed-> False,
AxesEdge->{{-1,-1},{1,-1},{-1,-1}},AxesStyle->Directive[Thick,18,Black],
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}(E)/t",Magnification->1.5],(π/2)]}
],
ListPointPlot3D[{{0,1,2}},PlotStyle->Black]]]

RHF and UHF

We will need to identify analytic forms for the  RHF and UHF solutions to use in the RMP and UMP expansions. 
We’ve done this to death before, so here it is in a relatively un-interesting form...

Ehf[((θα)/2),((θβ)/2)]/.{δϵ->0}//FullSimplify

Solve[{D[(1/2) (U+U Cos[θα] Cos[θβ]-2 t (Sin[θα]+Sin[θβ])),θα]==0,D[(1/2) (U+U Cos[θα] Cos[θβ]-2 t (Sin[θα]+Sin[θβ])),θβ]==0},{θα,θβ}]/.{ (C[1])->0, (C[2])->0}

{RHF1,UHF1,UHF2,RHF2}={{θα->-(π/2),θβ->-(π/2)},{θα->-(π/2),θβ->(π/2)},{θα->(π/2),θβ->-(π/2)},{θα->(π/2),θβ->(π/2)}};
{sbRHF1,sbUHF1,sbRHF2,sbUHF2}={{θα-> ArcTan[-((Sqrt[-4 (t^2)+(U^2)])/U),-((2 t)/U)],θβ-> ArcTan[-((Sqrt[-4 (t^2)+(U^2)])/U),-((2 t)/U)]},{θα-> ArcTan[-((Sqrt[-4 (t^2)+(U^2)])/U),((2 t)/U)],θβ-> ArcTan[((Sqrt[-4 (t^2)+(U^2)])/U),((2 t)/U)]},{θα->ArcTan[((Sqrt[-4 (t^2)+(U^2)])/U),-((2 t)/U)],θβ->ArcTan[((Sqrt[-4 (t^2)+(U^2)])/U),-((2 t)/U)]},{θα-> ArcTan[((Sqrt[-4 (t^2)+(U^2)])/U),((2 t)/U)],θβ->ArcTan[-((Sqrt[-4 (t^2)+(U^2)])/U),((2 t)/U)]}};

(1/2) (U+U Cos[θα] Cos[θβ]-2 t (Sin[θα]+Sin[θβ]))/.{RHF1,UHF1,UHF2,RHF2,sbRHF1,sbUHF1,sbUHF2,sbRHF2}//FullSimplify

With[{t=1},Plot[{(1/2) (4 t+U),(U/2),(U/2),(1/2) (-4 t+U),((2 (t^2))/U)+U,-((2 (t^2))/U),-((2 (t^2))/U),((2 (t^2))/U)+U},{U,-10,10}]]

Self-Consistency as a perturbation

We consider the introduction of the two-electron interaction as a perturbation to the one-electron core Hamiltonian.
For simplicity, we consider only the ground-state RHF and UHF solutions, with energies given by

Show[
Plot[{(1/2) (-4 +U ),If[U >2,-((2 )/(U ))]},{U,0,10},
PlotRangePadding->0,FrameStyle->Directive[Thick,18,Black],AspectRatio->0.6,
PlotStyle->{Orange,Darker[Green],Red,Blue},
PlotLegends->Placed[{
MaTeX["\\text{RHF}",Magnification->1.5],
MaTeX["\\text{UHF}",Magnification->1.5]},{0.9,0.6}],
FrameLabel->{
MaTeX["U / t",Magnification->1.5],
MaTeX["E_{\\text{HF}} /t",Magnification->1.5]},
ImageSize->Large,PlotTheme->"Scientific"],
ListPlot[{{2,-1}},PlotStyle->Black]
]

If we now consider the two-electron interaction as a perturbation to the one-electron core Hamiltonian,
we introduce the transformation U -> λ U.
In this case, the corresponding energies become

(Subscript[ϵ, RHF])[U_,t_,λ_]=(1/2) (-4 t+U λ);
(Subscript[ϵ, UHF])[U_,t_,λ_]=If[U λ>2(t^2),-((2 (t^2))/(U λ))];

We can then consider the UHF orbital parameters as a perturbation expansion of the one-electron orbital 
coefficients. If we consider the the complex λ plane, and introduce a complex-analytic continuation of 
the HF orbitals, we realise that the coalescence point corresponds to 
an exceptional point in both the orbital rotation angles and the HF energies.

With[{U=2,t=1},
Show[Plot3D[ Re[(1/π)ArcTan[-((Sqrt[-4 (t^2)+(U^2)(((rλ+I zλ))^2)])/(U(rλ+I zλ))),((2 t)/(U(rλ+I zλ)))]],{rλ,0.01,2},{zλ,-1,1},
PlotRange->{0,1},PlotStyle->Opacity[0.8],PlotPoints->50,MaxRecursion->6,Mesh->5,
ColorFunction->Function[{rλ,zλ},Hue[(1/(2π))Arg[(1/π)ArcTan[-((Sqrt[-4 (t^2)+(U^2)(((rλ+I zλ))^2)])/(U(rλ+I zλ))),((2 t)/(U(rλ+I zλ)))]-(1/2)],0.8,1.0]],ColorFunctionScaling->False,
BoxRatios->{1, 1, 1},ImageSize->Large,Boxed->False,
AxesEdge->{{-1,-1},{1,-1},{-1,-1}},AxesStyle->Directive[Thick,18,Black],
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}\\left(\\theta_{\\alpha}^{\\text{UHF}} \\right)/\\pi",Magnification->1.5],(π/2)]}],
Plot3D[ Re[(1/π)ArcTan[((Sqrt[-4 (t^2)+(U^2)(((rλ+I zλ))^2)])/(U(rλ+I zλ))),((2 t)/(U(rλ+I zλ)))]],{rλ,0.01,2},{zλ,-1,1},PlotRange->{0,1},
PlotRange->{0,1},PlotStyle->Opacity[0.8],PlotPoints->50,MaxRecursion->6,Mesh->5,
ColorFunction->Function[{rλ,zλ},Hue[(1/(2π))Arg[(1/π)ArcTan[((Sqrt[-4 (t^2)+(U^2)(((rλ+I zλ))^2)])/(U(rλ+I zλ))),((2 t)/(U(rλ+I zλ)))]-(1/2)],0.8,1.0]],ColorFunctionScaling->False,
BoxRatios->{1, 1, 1},ImageSize->Large],
Plot3D[0.5,{rλ,0.01,2},{zλ,-1,1},PlotRange->{0,1},Mesh->5,
PlotRange->{0,1},PlotStyle->{Opacity[0.75],RGBColor[0.9,0.9,0.99]},PlotPoints->50,MaxRecursion->6],
ListPointPlot3D[{{1,0,0.5}},PlotStyle->Black]

]]

With[{U=2,t=1},
Show[
Plot3D[ Re[(1/2) (-4 t+U (rλ+I zλ))],{rλ,0.01,2.7},{zλ,-1,1},
ColorFunction->Function[{rλ,zλ},Hue[-(1/(2π))Arg[ (1/2) (-4 t+U (rλ+I zλ))],0.75,1.0]],ColorFunctionScaling->False,
BoxRatios->{1, 1, 1},ImageSize->Large,PlotStyle->Opacity[0.7],Boxed->False,Mesh->5,PlotRange->{{0,3},{-1,1},{-2,1}},
AxesEdge->{{-1,-1},{1,-1},{-1,-1}},AxesStyle->Directive[Thick,18,Black],
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}(E_{\\text{HF}})/t",Magnification->1.5],(π/2)]}],
Plot3D[Re[ -((2 (t^2))/(U (rλ+I zλ)))],{rλ,0.01,3},{zλ,-1,1},PlotRange->{-2,2},Mesh->5,
ColorFunction->Function[{rλ,zλ},Hue[-(1/(2π))Arg[ -((2 (t^2))/(U (rλ+I zλ)))],0.75,1]],ColorFunctionScaling->False,
BoxRatios->{1, 1, 1},ImageSize->Large,PlotStyle->Opacity[0.7]],
ListPointPlot3D[{{1,0,-0.98}},PlotStyle->Black]
]]

With[{U=3,t=1},
Plot[{(1/2) (-4 t+U λ),If[U λ>2(t^2),-((2 (t^2))/(U λ))]},{λ,0,2}]
]

If we look at the optimal values of the orbital rotation angles in the complex plane, we find that
the Coulson-Fischer point represents an EP.

{Re[Exp[2π I t]],Im[Exp[2π I t]],0.2}/.{t->0.5}

RGBColor[Abs[Re[Exp[2π I t]]],Im[Exp[2π I t]],0.2]/.{t->0.5}

RGBColor[0.5,0.7,0.5]

RMP

We can now consider the RMP series convergence. To do this, we get the Fock operator in RMP as our reference and

(Subscript[H, RMP])=Hmp[(π/4),(π/4)]/.{δϵ->0}//FullSimplify//MatrixForm

Eigenvalues[(({{-2 t+U-((U λ)/2), 0, 0, ((U λ)/2)}, {0, U-((U λ)/2), ((U λ)/2), 0}, {0, ((U λ)/2), U-((U λ)/2), 0}, {((U λ)/2), 0, 0, 2 t+U-((U λ)/2)}}))]

First, we identify that the radius of convergence is defined by the exceptional points at λ = ± 4t / U. 
We illustrate this by plotting the energy expansions at orders of λ just below and above the radius of convergence

SeriesCoefficient[(1/2) (2 U-U λ-(Sqrt[16 (t^2)+(U^2) (λ^2)]))/.{U->3},{λ,0,0}]

With[{t=1},ListPlot[
{Table[{i,SeriesCoefficient[(1/2) (2 U-U λ-(Sqrt[16 (t^2)+(U^2) (λ^2)]))/.{U->3.5},{λ,0,i}]},{i,0,50}],
Table[{i,SeriesCoefficient[(1/2) (2 U-U λ-(Sqrt[16 (t^2)+(U^2) (λ^2)]))/.{U->4.5},{λ,0,i}]},{i,0,50}]},
PlotStyle->{Darker[Green],Red},PlotRange->{-5,5},PlotTheme->"Scientific",PlotRangePadding->0,Joined->True,PlotMarkers->Automatic,
AspectRatio->0.7,FrameStyle->Directive[Thick,18,Black],
FrameLabel->{MaTeX["n",Magnification->1.5],MaTeX["E^{(n)}_\\text{RMP}/t",Magnification->1.5]},
ImageSize->Large,
PlotLegends->Placed[{MaTeX["U = 3.5 t",Magnification->1.5],MaTeX["U = 4.5 t",Magnification->1.5]},{0.2,0.85}]
]]

We can also plot the corresponding Riemann surface that connects the reference and exact values at the same two U/t ratios

With[{t=1,U=3.5},
Show[{Plot3D[{Re[(1/2) (2 U-U (rλ+I zλ)-(Sqrt[16 (t^2)+(U^2) (((rλ+I zλ))^2)]))],Re[ U-U (rλ+I zλ)],Re[(1/2) (2 U-U (rλ+I zλ)+(Sqrt[16 (t^2)+(U^2) (((rλ+I zλ))^2)]))]},{rλ,-1.5,1.5},{zλ,0,2.5},
BoxRatios->{1, 1, 1},ImageSize->Large,PlotStyle->Opacity[0.7],Mesh->5, ViewPoint->{1.3,-1.6,0.7},
PlotRange-> {-5,10},Boxed->False,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}(E_\\text{FCI})/t",Magnification->1.5],(π/2)]},
AxesStyle->Directive[Thick,18,Black]],
ListPointPlot3D[{{0,4t/U, U }},PlotStyle->Black],
Graphics3D[{Opacity[0.5],Cylinder[{{0,0,-10},{0,0,20}},1]}]},
PlotRange->{{-1.2,1.2},{0,2.4},{-2,10}}]]

With[{t=1,U=4.5},
Show[{Plot3D[{Re[(1/2) (2 U-U (rλ+I zλ)-(Sqrt[16 (t^2)+(U^2) (((rλ+I zλ))^2)]))],Re[ U-U (rλ+I zλ)],Re[(1/2) (2 U-U (rλ+I zλ)+(Sqrt[16 (t^2)+(U^2) (((rλ+I zλ))^2)]))]},{rλ,-1.5,1.5},{zλ,0,2.4},
BoxRatios->{1, 1, 1},ImageSize->Large,PlotStyle->Opacity[0.7],Mesh->5,ViewPoint->{1.3,-1.6,0.7},PlotRange-> {-5,12},Boxed->False,
AxesEdge->{{-1,-1},{1,-1},{-1,-1}},
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}(E_\\text{RMP})/t",Magnification->1.5],(π/2)]},
AxesStyle->Directive[Thick,18,Black]],
ListPointPlot3D[{{0,4t/U, U }},PlotStyle->Black],
Graphics3D[{Opacity[0.5],Cylinder[{{0,0,-10},{0,0,20}},1]}]},
PlotRange->{{-1.2,1.2},{0,2.4},{-5,12}}
]]

An interesting connection...

Let’s take a quick look at the exceptional points we have found so far...

FCI:             1   = ± i 4t / U
RMP:          (Subscript[λ, c])= ± i 4t / U
CFP:           1   = 2t / U

So we see that the magnitudes of each exceptional point are all closely related.
Is it possible that each of these exceptional point is dictated by the fundamental exceptional point in the FCI energy?
Or is this an artefact of our incredible simple model?

UMP

We can now consider the MP expansion from the symmetry-broken UHF reference state.
The optimal rotation angles are computed in the HF section and are given as...

sbUHF1

We use these to construct the UMP Hamiltonian problem ...

Hump[U_,t_,λ_]=Hmp[(1/2) ArcTan[-((Sqrt[-4 (t^2) +(U^2)])/U),((2 t)/U)],(1/2) ArcTan[((Sqrt[-4 (t^2)+(U^2)])/U),((2 t)/U)]]/.{δϵ->0}//FullSimplify;
Hump[U,t,λ]//MatrixForm

As expected, the spatial and spin contamination makes the exact eigenvalues more difficult to work with and introduces many more
coupling elements into the expression. We get closed forms, but they are quite unpleasant!

ToRadicals[Eigenvalues[(({{-((2 (t^2) λ)/U), 0, 0, ((2 (t^2) λ)/U)}, {0, U-((2 (t^2) λ)/U), ((2 (t^2) λ)/U), ((2 t (Sqrt[-4 (t^2)+(U^2)]) λ)/U)}, {0, ((2 (t^2) λ)/U), U-((2 (t^2) λ)/U), -((2 t (Sqrt[-4 (t^2)+(U^2)]) λ)/U)}, {((2 (t^2) λ)/U), ((2 t (Sqrt[-4 (t^2)+(U^2)]) λ)/U), -((2 t (Sqrt[-4 (t^2)+(U^2)]) λ)/U), -2 U (-1+λ)+((6 (t^2) λ)/U)}}))]]//FullSimplify

ϵSingletClosed[U_,t_,λ_]={
(1/(3 U))(3 (U^2)-2 (U^2) λ+((-12 (t^2) (U^2) (-2+λ) λ+(U^4) (3-6 λ+4 (λ^2)))/((((U^6) λ (-9+2 (9-4 λ) λ)+18 (t^2) (U^4) λ (3+2 (-2+λ) λ)+(1/2) (Sqrt[4 (U^8) (λ^2) ((((U^2) (-9+2 (9-4 λ) λ)+18 (t^2) (3+2 (-2+λ) λ)))^2)-4 (((-12 (t^2) (U^2) (-2+λ) λ+(U^4) (3-6 λ+4 (λ^2))))^3)])))^(1/3)))+((((U^6) λ (-9+2 (9-4 λ) λ)+18 (t^2) (U^4) λ (3+2 (-2+λ) λ)+(1/2) (Sqrt[4 (U^8) (λ^2) ((((U^2) (-9+2 (9-4 λ) λ)+18 (t^2) (3+2 (-2+λ) λ)))^2)-4 (((-12 (t^2) (U^2) (-2+λ) λ+(U^4) (3-6 λ+4 (λ^2))))^3)])))^(1/3))),U-((2 U λ)/3)+(((1+I (Sqrt[3])) U (12 (t^2) (-2+λ) λ+(U^2) (-3+6 λ-4 (λ^2))))/(6 ((((U^6) λ (-9+2 (9-4 λ) λ)+18 (t^2) (U^4) λ (3+2 (-2+λ) λ)+(1/2) (Sqrt[4 (U^8) (λ^2) ((((U^2) (-9+2 (9-4 λ) λ)+18 (t^2) (3+2 (-2+λ) λ)))^2)-4 (((-12 (t^2) (U^2) (-2+λ) λ+(U^4) (3-6 λ+4 (λ^2))))^3)])))^(1/3))))+((I (I+(Sqrt[3])) ((((U^6) λ (-9+2 (9-4 λ) λ)+18 (t^2) (U^4) λ (3+2 (-2+λ) λ)+(1/2) (Sqrt[4 (U^8) (λ^2) ((((U^2) (-9+2 (9-4 λ) λ)+18 (t^2) (3+2 (-2+λ) λ)))^2)-4 (((-12 (t^2) (U^2) (-2+λ) λ+(U^4) (3-6 λ+4 (λ^2))))^3)])))^(1/3)))/(6 U)),U-((2 U λ)/3)+((I (I+(Sqrt[3])) U (-12 (t^2) (-2+λ) λ+(U^2) (3-6 λ+4 (λ^2))))/(6 ((((U^6) λ (-9+2 (9-4 λ) λ)+18 (t^2) (U^4) λ (3+2 (-2+λ) λ)+(1/2) (Sqrt[4 (U^8) (λ^2) ((((U^2) (-9+2 (9-4 λ) λ)+18 (t^2) (3+2 (-2+λ) λ)))^2)-4 (((-12 (t^2) (U^2) (-2+λ) λ+(U^4) (3-6 λ+4 (λ^2))))^3)])))^(1/3))))-(((1+I (Sqrt[3])) ((((U^6) λ (-9+2 (9-4 λ) λ)+18 (t^2) (U^4) λ (3+2 (-2+λ) λ)+(1/2) (Sqrt[4 (U^8) (λ^2) ((((U^2) (-9+2 (9-4 λ) λ)+18 (t^2) (3+2 (-2+λ) λ)))^2)-4 (((-12 (t^2) (U^2) (-2+λ) λ+(U^4) (3-6 λ+4 (λ^2))))^3)])))^(1/3)))/(6 U))};

(We can probably get analytic forms for the exceptional points using these expressions).

In the subsequent analysis, it is easier to work with the Root forms for the solutions that are better manipulated by the SeriesCoefficient and Plot3D routines.

ϵSinglet[U_,t_,λ_]={
((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,1])/U),
((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,2])/U),
((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,3])/U)};

UMP Convergence Behaviour

Here we compare the convergence behaviour for the MP ground- and excited-state expansions using UMP ground-state orbitals.

Ground State Convergence

Here we should the UMP contributions for the ground state at two different U/t ratios. The first (U = 3t) is well within the RMP convergence region, while the second (U = 7t) falls outside.
Using the UHF orbitals, we see convergence at both U/t ratios, although the rate of convergence becomes increasingly slow for larger U. 
This behaviour is reflected in the position of the exceptional point (see below).

ListLogPlot[
{Table[{i,Abs[SeriesCoefficient[((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,1])/U)/.{U->3,t->1},{λ,0,i}]]},{i,0,200}],
Table[{i,Abs[SeriesCoefficient[((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,1])/U)/.{U->7,t->1},{λ,0,i}]]},{i,0,200}]},
PlotStyle->{Darker[Green],Red},PlotRange->{(10^(-8)),1},
PlotTheme->"Scientific",PlotRangePadding->0,Joined->True,PlotMarkers->Automatic,
AspectRatio->0.7,FrameStyle->Directive[Thick,18,Black],
FrameLabel->{MaTeX["n",Magnification->1.5],MaTeX["E^{(n)}_\\text{UMP} /t",Magnification->1.5]},
ImageSize->Large,PlotLabel->MaTeX["\\text{Ground State}",Magnification-> 1.5],
PlotLegends->Placed[{MaTeX["U = 3 t",Magnification->1.5],MaTeX["U = 7 t",Magnification->1.5]},{0.2,0.85}]
]

Open-Shell Singlet Convergence

For the open-shell singlet (the first excited state), we see catastrophic divergence at both U/t ratios. 
In fact we will get convergence for ALL U/t ratios using the UMP orbitals (ie beyond the CFP).

ListLogPlot[
{Table[{i,Abs[SeriesCoefficient[((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,2])/U)/.{U->3,t->1},{λ,0,i}]]},{i,0,50}],
Table[{i,Abs[SeriesCoefficient[((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,2])/U)/.{U->7,t->1},{λ,0,i}]]},{i,0,50}]},
PlotStyle->{Darker[Green],Red},PlotRange->{(10^(-1)),(10^12)},
PlotTheme->"Scientific",PlotRangePadding->0,Joined->True,PlotMarkers->Automatic,
AspectRatio->0.7,FrameStyle->Directive[Thick,18,Black],
FrameLabel->{MaTeX["n",Magnification->1.5],MaTeX["E^{(n)}_\\text{UMP} /t",Magnification->1.5]},
ImageSize->Large,PlotLabel->MaTeX["\\text{Open-Shell Singlet}",Magnification-> 1.5],
PlotLegends->Placed[{MaTeX["U = 3 t",Magnification->1.5],MaTeX["U = 7 t",Magnification->1.5]},{0.2,0.85}]
]

Double Excitation Convergence

The same catastrophic convergence is seen for the double excitation, in fact with an almost identical pattern and oscillation period (spikes due to zero-crossing).
These two states share the same exceptional point, suggesting that the position of this point controls the oscillations in the radius of convergence. 
I suspect this is covered in the idea of alpha/beta and oscillatory convergence elsewhere?

ListLogPlot[
{Table[{i,Abs[SeriesCoefficient[((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,3])/U)/.{U->3,t->1},{λ,0,i}]]},{i,0,50}],
Table[{i,Abs[SeriesCoefficient[((Root[4 (t^2) (U^4) λ-4 (t^2) (U^4) (λ^2)+(2 (U^4)-8 (t^2) (U^2) λ-2 (U^4) λ+4 (t^2) (U^2) (λ^2)) #1+(-3 (U^2)+2 (U^2) λ) ((#1)^2)+((#1)^3)&,3])/U)/.{U->7,t->1},{λ,0,i}]]},{i,0,50}]},
PlotStyle->{Darker[Green],Red},PlotRange->{0.1,(10^12)},
PlotTheme->"Scientific",PlotRangePadding->0,Joined->True,PlotMarkers->Automatic,
AspectRatio->0.7,FrameStyle->Directive[Thick,18,Black],
FrameLabel->{MaTeX["n",Magnification->1.5],MaTeX["E^{(n)}_\\text{UMP} /t",Magnification->1.5]},
ImageSize->Large,PlotLabel->MaTeX["\\text{Double Excitation}",Magnification-> 1.5],
PlotLegends->Placed[{MaTeX["U = 3 t",Magnification->1.5],MaTeX["U = 7 t",Magnification->1.5]},{0.2,0.85}]
]

UMP Riemann Surfaces

In this section, we plot the Riemann surfaces for the UMP expansions at the two U/t ratios considered above.

U = 7t

Here we are well towards the strong correlation regime, where it is known that UMP will be slowly convergent but RMP diverges.
We see a single exceptional point on the ground-state surface which falls just outside (maybe on?) the radius of convergence. 
An exceptional point this close to the radius of convergence gives an increasingly slow convergence of the UMP series, but it will converge evenutally!

Note that the ground-state energy is remarkably flat since the UHF energy is a pretty good estimate. 
Most of the UMP expansion is actually correcting the spin-contamination in the wave function.

On the other hand, there is an exceptional point on the excited energy surface that is well within the radius of convergence.
We can therefore say that the use of a symmetry-broken UHF wave function can retain a convergent ground-state perturbation series
at the expense of a divergent excited-state perturbation series. (Note: the orbitals are not optimised for excited-state here).
In contrast, the RMP expansion was always convergent for the open-shell excited state (which was a single CSF) while
the radius of convergence for the doubly-excited state was identical to the ground-state as this was the only exceptional point.

BoxRatios->{1, 1, 1},ImageSize->Large,PlotStyle->Opacity[0.7],ViewPoint->{1.3,-1.6,0.7},PlotRange-> {-5,12},Boxed->False,
AxesEdge->{{-1,-1},{1,-1},{-1,-1}},
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}(E_\\text{RMP})/t",Magnification->1.5],(π/2)]},
AxesStyle->Directive[Thick,18,Black]]

Show[
{Plot3D[Evaluate[Re[ϵSinglet[7,1,rλ+I zλ]]],{rλ,-0,1.5},{zλ,0,1},PlotStyle->Opacity[0.7],BoxRatios->{1, 1, 1},
ImageSize->Large,MaxRecursion->6,Mesh->5,Boxed->False,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}(E_\\text{UMP})/t",Magnification->1.5],(π/2)]},
AxesStyle->Directive[Thick,18,Black]
],
Graphics3D[{Opacity[0.5],Cylinder[{{0,0,-4},{0,0,20}},1]}]},PlotRange->{{0,1.2},{0,1.2},{-2.4,17}},Boxed->False
]

U = 3t

If we now consider a case within the radius of convergence for RMP, we see that the ground-state exceptional point has moved
slightly further beyond the radius of convergence. However, this is still much closer than the RMP exceptional point, so the convergence
is now slower than the RMP convergence.

On the other hand, the excited-state exceptional point is still within the radius of convergence and its expansion will therefore diverge.

Show[
{Plot3D[Evaluate[Re[ϵSinglet[3,1,rλ+I zλ]]],{rλ,-00,1.2},{zλ,0,1},PlotStyle->Opacity[0.7],BoxRatios->{1, 1, 1},ImageSize->Large,MaxRecursion->6,
ImageSize->Large,MaxRecursion->6,Mesh->{5,5},Boxed->False,AxesEdge->{{-1,-1},{1,-1},{-1,-1}},
AxesLabel->{
MaTeX["\\text{Re}(\\lambda)",Magnification->1.5],
MaTeX["\\text{Im}(\\lambda)",Magnification->1.5],
Rotate[MaTeX["\\text{Re}(E_\\text{UMP})/t",Magnification->1.5],(π/2)]},
AxesStyle->Directive[Thick,18,Black]],
Graphics3D[{Opacity[0.5],Cylinder[{{0,0,-4},{0,0,20}},1]}]},PlotRange->{{0,1.2},{0,1.2},{-2,6.}},Boxed->False
]

Will the UMP ground-state always converge?

This important question can be paraphrased as “does the ground-state EP ever fall within the radius of convergence?”. 
To answer this, we can consider the fixed FCI point on the corresponding Riemann surface

Firstly the nature of the symmetry-breaking in UHF will always give an energy that is an increasingly accurate estimate of the exact energy. 
This is because UHF can capture the strong correlation effects by spatially separating the two electrons.
At the fixed point λ=1, the energies of the states will always be the FCI energy, and this is independent of the reference orbitals.

If the ground-state is a good FCI estimate, it will cross from λ=0 to 1 in a “flat” fashion. This makes it far less likely to interact with 
a higher energy state, so it is less likely to encounter an exceptional point. Equally, the first excited state will tend towards the 
exact first-excited singlet FCI energy. It is possible for these states to cross on the real-λ line, but any crossing must occur
beyond λ=1 because otherwise they would not connect to the correct fixed point.

So what happens to the RMP exceptional point?

In the RMP case, the convergence is controlled by an exceptional point between the ground- and doubly-excited RHF configurations.
These two states interact strongly in the λ expansion in the large U/t limit because the ground RHF state is heavily 
contaminated by the 3rd FCI singlet  state. It is this strong interaction that causes the exceptional point, and ultimately
the divergence. 

In the UMP case, this “energetic” contamination from the 3rd FCI state is reduced (removed?) by replacing it with triplet contamination.
The strong interaction is then shifted from the ground configuration to the “open-shell” configurations in the UHF framework, 
since the 3rd FCI singlet state contributions are rotated into the virtual orbitals.
This explains why we see the severe excited-state exceptional point in the UMP Riemann surface.
In a way, we can think of this as the “open-shell” configurations “shielding” the ground-state from the strong interaction (and
exceptional point), hence we have more robust convergence in ground-state, but less robust convergence in excited-state.

Summary

SeriesCoefficient[((-2(t^2))/(U λ)),{λ,0,5}]