Bootstrapping Quantum Mechanics

(*potential*)

V[x]=g (x^2)+(x^4)

(*deriving the recurrence relation*)

(x^a)V'[x]+2a (x^(a-1))V[x]-2a e (x^(a-1))-((a(a-1)(a-2))/2)(x^(a-3))/.{V->((g ((#)^2)+((#)^4))&)}

Collect[Expand[-(1/2) (-2+a) (-1+a) a (x^(-3+a))-2 a e (x^(-1+a))+(x^a) (2 g x+4 (x^3))+2 a (x^(-1+a)) (g (x^2)+(x^4)),x],(x^(_)),Factor]

(x^(3+a))=((-1)/(2(a+2)))(-(1/2) (-2+a) (-1+a) a (x^(-3+a))-2 a e (x^(-1+a))+2 (1+a) g (x^(1+a)))

Collect[Expand[((-1)/(2(a+2)))(-(1/2) (-2+a) (-1+a) a (x^(-3+a))-2 a e (x^(-1+a))+2 (1+a) g (x^(1+a)))/.{a->a-3},x],(x^(_)),Factor]

(*final form of the recurrence*)

(x^a)=(((-5+a) (-4+a) (-3+a) (x^(-6+a)))/(4 (-1+a)))+(((-3+a) e (x^(-4+a)))/(-1+a))-(((-2+a) g (x^(-2+a)))/(-1+a))

(*a=4 requires special treatment (otherwise the recurrence tries to evaluate <(x^(-2))>)*)

(((-5+a) (-4+a) (-3+a) (x^(-6+a)))/(4 (-1+a)))+(((-3+a) e (x^(-4+a)))/(-1+a))-(((-2+a) g (x^(-2+a)))/(-1+a))/.{a->4}

Setting up the recurrence

(*ev[g][a] = <(x^a)>, <(x^2)> = t*)

ev[g_][a_]:=0/;OddQ[a];
ev[g_][0]:=1;
ev[g_][2]:=t;
ev[g_][4]:=(e/3)-((2 g )/3)t;
ev[g_][a_]:=(ev[g][a]=Factor[(((-5+a) (-4+a) (-3+a) )/(4 (-1+a)))ev[g][a-6]+(((-3+a) e )/(-1+a))ev[g][a-4]-(((-2+a) g )/(-1+a))ev[g][a-2]]);

Imposing positivity

(*hankel Matrix: (Subscript[M, ab]) = <(x^(a+b))>*)

hankelMatrix[g_,aMax_][ee_,tt_]:=Table[ev[g][a+b]/.{e->ee,t->tt},{a,0,aMax},{b,0,aMax}];

isAllowed[g_,aMax_][ee_,tt_]:=(isAllowed[g,aMax][ee,tt]=(hankelMatrix[g,aMax][ee,tt]//Eigenvalues//Min)>=0);

Bootstrapping

(*initial scan over a broad range of the (e,t) parameter space*)

Module[{g=-5,aMax=11},
Monitor[Select[Flatten[Table[{ee,tt,isAllowed[g,aMax][SetPrecision[ee,30],SetPrecision[tt,30]]},{ee,-8,8,(1/10)},{tt,0,5,(1/10)}],{1,2}],(#[[3]]==True)&][[All,1;;2]]//ListPlot,{ee,tt}//N]
]

(*zooming-in on the ground state + (1^st) excited state island*)

Module[{g=-5,aMax=11},
Monitor[Select[Flatten[Table[{ee,tt,isAllowed[g,aMax][SetPrecision[ee,30],SetPrecision[tt,30]]},{ee,((-39)/10),((-30)/10),(1/100)},{tt,(15/10),(24/10),(1/100)}],{1,2}],(#[[3]]==True)&][[All,1;;2]],{ee,tt}//N]//ListPlot
]

(*increasing the number of constraints to (Subscript[a, max]) = 15: still not enough to tell the ground state and (1^st) excited state apart*)

Module[{g=-5,aMax=15},
Monitor[Select[Flatten[Table[{ee,tt,isAllowed[g,aMax][SetPrecision[ee,30],SetPrecision[tt,30]]},{ee,((-36)/10),((-32)/10),(1/100)},{tt,(18/10),(23/10),(1/100)}],{1,2}],(#[[3]]==True)&][[All,1;;2]],{ee,tt}//N]//ListPlot
]

(*for (Subscript[a, max])=19, we finally get two nice islands around the two lowest states*)

Module[{g=-5,aMax=19},
Monitor[Select[Flatten[Table[{ee,tt,isAllowed[g,aMax][SetPrecision[ee,30],SetPrecision[tt,30]]},{ee,-(3415/(10^3)),-(3245/(10^3)),(1/500)},{tt,(2020/(10^3)),(2225/(10^3)),(1/500)}],{1,2}],(#[[3]]==True)&][[All,1;;2]],{ee,tt}//N]//ListPlot[#,PlotRange->{{-4,-2.5},{1.5,2.5}}]&
]