Implement 2D SHO lowering and raising operator functions, and use them to calculate (x^4)+ (y^4), xy, and (x^2)(y^2) operator actions on some two particle states, and then some specific calculations for ps8 p3.

Methods.

ClearAll[raiseOrLower,raise, lower, raiseKet, lowerKet, k, l, k1, x0, xOp, yOp,xyOp,normalize, sumByList , display, formatKet, formatKetList, isNonZero, hbar, omega, m, flattenKet, debugMode, debugPrint, sopts, projection, rejection, sumKets, diffKets, lambda]

On[Assert] ;
debugMode = False ;
(*debugMode = True ;*)

(*3 // debugPrint ;
debugMode = True ;
4 // debugPrint ;

debugMode = False*)
debugPrint := If[ debugMode == True, # // Print,,] & ;

(* Given a ket |k> = {v,(Subscript[n, 1]),(Subscript[n, 2])} representing a single basis state v|(Subscript[n, 1]),(Subscript[n, 2])> apply:
(a) sign=-1, delta=0: the lowering operator: (Subscript[a, p]), as
(Subscript[a, 1])|k> = v (Sqrt[(Subscript[n, 1])])|(Subscript[n, 1])-1,(Subscript[n, 2])>, or (Subscript[a, 2])|k> = v (Sqrt[(Subscript[n, 2])])|(Subscript[n, 1]),(Subscript[n, 2])-2>.  Handling of |0,0>, |n,0>, and |0,n> are special cased. 
(b) sign=+1, delta=+1: Apply (SuperDagger[(Subscript[a, p])])|k> to |k> = {v,(Subscript[n, 1]),(Subscript[n, 2])} representation.
*) 
raiseOrLower[{v_,n1_Integer, n2_Integer}, p_Integer, sign_Integer, delta_Integer] := Module[{nn,n},
nn = {n1, n2} ;
n = nn[[p]] ;
((p == 1) || (p == 2)) // Assert ;

{Sqrt[n + delta]v, n1 + sign KroneckerDelta[p,1],n2 + sign KroneckerDelta[p,2] } 
] ;

lower[{v_,n1_Integer, n2_Integer}, p_Integer] :=  raiseOrLower[{v,n1, n2}, p, -1, 0];
raise[{v_,n1_Integer, n2_Integer}, p_Integer] :=  raiseOrLower[{v,n1, n2}, p, +1, 1];
lower[{v_,0, 0}, n_Integer] := {0, 0, 0} ;
lower[{v_,0, n_Integer}, 1] := {0, 0, 0} ;
lower[{v_,n_Integer, 0}, 2] := {0, 0, 0} ;


(* Given a ket |k> = (∑_i)(Subscript[v, i]) |(Subscript[n, 1,i]),(Subscript[n, 2,i])>, represented as a list of kets{{(Subscript[v, 1]), (Subscript[n, 1,1]), (Subscript[n, 2,1])},{(Subscript[v, 2]), , (Subscript[n, 1,2]), (Subscript[n, 2,2])}, ...}, apply the raising ( lowering) operators to each element.  Zero elements that result from the lowering operator are not weeded out by these functions. *)
raiseKet[k_, p_Integer] := (raise[#, p] &/@ k) ;
(*raiseKet[{{1,0,0}}, 1]*)
lowerKet[k_, p_Integer] := (lower[#, p] &/@ k) ;

(*Given a list
l = {{(Subscript[v, 1]),(Subscript[n, 1,1]),(Subscript[n, 2,1])}, {(Subscript[v, 1]),(Subscript[n, 1,2]),(Subscript[n, 2,2])}, ...} 
find and sum the (Subscript[v, i]) values for any (Subscript[k, j]) == k = {(Subscript[k, 1]),(Subscript[k, 2])}
*) 
sumByList[l_,k_] := #[[1]] &/@ Select[l, Drop[#,1] == k &] // Total ;

(* Test if a value is a symbol or a non-zero number *)
isNonZero = ((!NumberQ[#]) || (# ≠ 0)) & ;

(* {{(x0/(Sqrt[2])),{1,0}},{(x0/(Sqrt[2])),{1,1}}} => {{(x0/(Sqrt[2])),1,0},{(x0/(Sqrt[2])),1,1}} *)
flattenKet[k_] := {#[[1]],#[[2]][[1]], #[[2]][[2]]} &/@ k ;
(*flattenKet[{{(x0/(Sqrt[2])),{1,0}},{(x0/(Sqrt[2])),{1,1}}}]*)

(* Weed out any zeros and merge (sum) the coefficients for any kets that are duplicated.  Each of the resulting ket coefficients are multiplied by the optional scale.
 *)
normalize[k_List, scale_ : 1] := Module[{nonzero,basisElements, basisElementMultipliers, result},
nonzero = Select[k, isNonZero[#[[1]]] & ] ;
basisElements = Part[ nonzero , All,2;;3]// DeleteDuplicates ;
basisElementMultipliers = sumByList[nonzero, #] &/@ basisElements ;
result = {scale basisElementMultipliers, basisElements}  ;

Grid[{{ "k = ", k}, 
{"scale = ", scale}, 
{"nonzero = ", nonzero},
{"basisElements = ", basisElements},
{"basisElementMultipliers = ", basisElementMultipliers},
{"result = ", result}
} ] // debugPrint ;
result // Transpose // flattenKet
] ;

(*A display method for |k> = {v,n} *)
display  = {x0 -> (Subscript[x, 0]), hbar -> ℏ, omega -> ω, g1 -> "(Subscript[g, 1])", g2 -> "(Subscript[g, 2])", lambda -> λ} ;
sopts = {14, "Text"} ;

(* formatKet[{2, 3, 4}]  => 2|3,4> *)
formatKet[{v_,n1_Integer, n2_Integer}] := Row[{v /. display, "\[VerticalSeparator]", n1,", ",n2,"⟩"}] ;

(* formatKetList[{{2, 3, 4},{4,5,6}}] => 2|3,4> + 4|5,6> *)
formatKetList[k_] := (formatKet[#] &/@ k)  // Total;

sumKets[k1_, k2_, scale_ : 1] := Module[{both},
both = Flatten[{k1, k2}, 1] ;
both // debugPrint ;
normalize[ both, scale ]
] ;

diffKets[k1_, k2_, scale_ : 1] := Module[{both},
both = Flatten[{k1, normalize[k2,-1]}, 1] ;
both // debugPrint ;
normalize[ both, scale ]
] ;

(* Apply the X operator to a ket. *)
(*xOp[{{1,0,0},{1,0,1}}]*)
xyOp[k_, p_Integer ] := sumKets[ raiseKet[k, p], lowerKet[k, p], x0/Sqrt[2] ] ;
(*Module[{r, l, both},
r = raiseKet[k, p] ;
l = lowerKet[k, p] ;
both = Flatten[{r, l}, 1] ;
both // debugPrint ;
normalize[ both, x0/Sqrt[2]]
] ;*)

(* Apply the X,Y operator to a ket. *)
xOp[k_ ] := xyOp[k, 1] ;
yOp[k_ ] := xyOp[k, 2] ;

(* for two kets each with a nonzero projection onto state k,take the dot products of those projections (assuming real coefficients) *)
dotOne[k1_, k2_, k_] := Module[{c1, c2},

c1 = Select[k1, Drop[#,1]==k &]  ;
c2 = Select[k2, Drop[#,1]==k &]  ;
1 == Length[c1] // Assert ;
1 == Length[c2] // Assert ;

Grid[{
{"dotOne:",""},
{"k1 = ", k1}, {"k2 = ", k2},
{"c1 = ", c1}, {"c2 = ", c2},
{"k = ", k}
} ] // debugPrint ;

(c1[[1]][[1]])  (c2[[1]][[1]])
] ;
(*dotOne[{{11,2,3},{4,5,6}}, {{7,2,3},{8,7,6}},{2,3}]*)

(*compute the inner product of two kets *)
innerProd[k1_, k2_] := Module[{t1, t2, k, dotlist, r, tmp},
t1 = (Drop[#,1] &/@ k1) ;
t2 = (Drop[#,1] &/@ k2) ;

k = Intersection[ t1, t2] ;
tmp = {} ;
r= If [ k ≠ {},
tmp = dotOne[k1, k2, #] &/@ k ;
Total[ tmp ],
0 ] ;

Grid[{
{"innerProd:",""},
{"k1 = ", k1}, {"k2 = ", k2},
{"t1 = ", t1}, {"t2 = ", t2},
{"k = ", k},
{"tmp = ", tmp},
{"r = ", r}
} ] // debugPrint ;
r
] ;
(*innerProd[{{11,2,3},{4,5,6}}, {{7,2,3},{8,7,6},{3,5,6}}]*)

(* signed projection of b on a: (sign) ((|a><a|)/(<a|a>)) |b> *)
projection[a_, b_, sign_] := Module[{norma, scale},
norma = innerProd[a,a] ;
scale = innerProd[a,b] ;
normalize[a, sign scale/norma]
] ;

(*subtract any component of |a> from |b>*)
rejection[a_, b_] := normalize[sumKets[ b, projection[a,b,-1] ],1] ;

Calculation of  (x^4)+(y^4), xy, (x^n), (y^n) operator actions.

ClearAll[ketZero, xp, n, yp, ketOne, i, xFourthPlusyFourth, xTimesY, basis, xxTimesYy]
ketZero = {{1,0, 0}}  ;
debugMode = False ;
(*x0 = 1 ;*)
xFourthPlusyFourth[k_] :=
normalize[Flatten[{k // xOp // xOp // xOp// xOp,
k // yOp // yOp // yOp// yOp},1]] // Simplify ;
xTimesY[k_] := (k // xOp // yOp // Simplify) ;
xxTimesYy[k_] := (k // xOp //xOp// yOp //yOp // Simplify) ;

basis= {{{1,0,0}},{{1,1,0}},{{1,0,1}}} ;
Style[Grid[{
{"((x^4)+ (y^4))|0,0> = ", {{1,0,0}} // xFourthPlusyFourth // formatKetList},
{"((x^4)+ (y^4))|1,0> = ", {{1,1,0}} // xFourthPlusyFourth // formatKetList},
{"((x^4)+ (y^4))|0,1> = ", {{1,0,1}} // xFourthPlusyFourth // formatKetList},
{"x y|0,0> = ", {{1,0,0}} // xTimesY // formatKetList},
{"x y|1,0> = ", {{1,1,0}} // xTimesY // formatKetList},
{"x y|0,1> = ", {{1,0,1}} // xTimesY // formatKetList}(*,
{"(x^2)(y^2)|0,0> = ", {{1,0,0}} // xxTimesYy // formatKetList},
{"(x^2)(y^2)|1,0> = ", {{1,1,0}} // xxTimesYy // formatKetList},
{"(x^2)(y^2)|0,1> = ", {{1,0,1}} // xxTimesYy // formatKetList}*)
}] ,sopts] 
Style[Grid[{ 
{"basis for n ≤ 1 subspace:", formatKetList[#] &/@ basis},
{"<(Subscript[b, i])|((x^4)+ (y^4))|0,0> = ", (innerProd[{{1,0,0}} // xFourthPlusyFourth,#] &/@ basis )/. display},
{"<(Subscript[b, i])|((x^4)+ (y^4))|1,0> = ", (innerProd[{{1,1,0}} // xFourthPlusyFourth,#] &/@ basis)/. display},
{"<(Subscript[b, i])|((x^4)+ (y^4))|0,1> = ", (innerProd[{{1,0,1}} // xFourthPlusyFourth,#] &/@ basis)/. display},
{"<(Subscript[b, i])|x y|0,0> = ", (innerProd[{{1,0,0}} // xTimesY,#] &/@ basis)/. display},
{"<(Subscript[b, i])|x y|1,0> = ", (innerProd[{{1,1,0}} // xTimesY,#] &/@ basis)/. display},
{"<(Subscript[b, i])|x y|0,1> = ", (innerProd[{{1,0,1}} // xTimesY,#] &/@ basis)/. display}(*,
{"<(Subscript[b, i])|(x^2)(y^2)|0,0> = ", (innerProd[{{1,0,0}} // xxTimesYy,#] &/@ basis)/. display},
{"<(Subscript[b, i])|(x^2)(y^2)|1,0> = ", (innerProd[{{1,1,0}} // xxTimesYy,#] &/@ basis)/. display},
{"<(Subscript[b, i])|(x^2)(y^2)|0,1> = ", (innerProd[{{1,0,1}} // xxTimesYy,#] &/@ basis)/. display}*)
}] ,sopts]

n = 4 ;
xp = Range[n+1] ;
xp[[1]] = ketZero ;
For[i=1,i ≤ n,i++, 
xp[[i + 1]] = (xp[[i]] // xOp) ;
]

Style[Row[{(x^(#-1)), "|0,0> = ", xp[[#]] // formatKetList}] & /@ Range[n+1] // TableForm, sopts] 

(*yp = Range[n+1] ;
yp[[1]] = ketZero;
For[i=1,i ≤ n,i++, yp[[i + 1]] = yp[[i ]] // yOp // Simplify]

Style[Row[{(y^(#-1)), "|1> = ", yp[[#]] // formatKetList}] & /@ Range[n+1] // TableForm, sopts] *)

Calculation of  (x^4)+(y^4), xy pertubations.

ClearAll[ v1, v2, g1, g2, v1zero, rv1zero, scaleN1, energy,pv1zero, v2zero, vOne00, vTwo00 ]

v1[k_] :=  sumKets[xOp[xOp[xOp[xOp[k]]]],yOp[yOp[yOp[yOp[k]]]],hbar omega g1/x0^4 ] ;
v2[k_] := normalize[xOp[yOp[k]], hbar omega g2/x0^2] ;

(* energy of one component of a state *)
energy[k_] := hbar omega ((Drop[k,1] // Total) + 1) ;

(* rescale each component of the state |(Subscript[n, 1]),(Subscript[n, 2])> -> (1/(-ℏ ω ((Subscript[n, 1])+ (Subscript[n, 2])))) |(Subscript[n, 1]),(Subscript[n, 2])> *)
scaleN1[k_] := ({#[[1]]/(hbar omega - energy[#]), Drop[#,1]} & /@ k) // flattenKet ;

v1zero = v1[ketZero] ;
v2zero = v2[ketZero] ;
rv1zero = rejection[ketZero, v1zero];
pv1zero = scaleN1[rv1zero] ;
vOne00 = innerProd[ketZero,v1zero] ;
vTwo00 = innerProd[ketZero,v2zero] ;
vOneNprime = v1[pv1zero] ;
vOne01 = innerProd[ketZero,vOneNprime] ;

Style[{#1," = ", #2 /. display} &@@@{
{"(Subscript[V, 1])|0,0>", v1zero // formatKetList},
{"(Subscript[V, 2])|0,0>", v2zero // formatKetList},
{"(Δ^1) = <0,0|(Subscript[V, 1])|0,0>", vOne00},
{"<0,0|(Subscript[V, 2])|0,0>", vTwo00},
{"(Subscript[(OverBar[P]), 0, 0]) (Subscript[V, 1])|0,0>", rv1zero // formatKetList},
{"|0,0((>)^1)= \!\(\*FractionBox[SubscriptBox[OverscriptBox[\(P\), \(_\)], \(0, 0\)], \(\*SuperscriptBox[\(E\), \(0\)] - \\ \*SubscriptBox[\(H\), \(0\)]\)]\) (Subscript[V, 1])|0,0>", pv1zero // formatKetList},
{"(Subscript[V, 1])|0,0((>)^1)", vOneNprime // formatKetList},
{"<0,0|(Subscript[V, 1])|0,0((>)^1)", vOne01},
{"(Δ^2) = <0,0|(Subscript[V, 2])|0,0> + <0,0|(Subscript[V, 1])|0,0((>)^1)",vTwo00 + vOne01}
} // Grid ,sopts]

Degenerate calculation (also done by hand in the problem set report.)

Degenerate computation, forming the matrix element for the {|1,0>, |0,1>} subspace basis.

ClearAll[h]
h = (2 hbar omega + 9 g1 lambda hbar omega/2)IdentityMatrix[2] + (g2 lambda^2 hbar omega/2) PauliMatrix[1] ;


Style[{#1," = ", #2 /. display} &@@@{
{"H", h // MatrixForm },
{"Eigenvalues",Eigenvalues[h]//Column}
} // Grid ,sopts]