Implement 1D SHO lowering and raising operator functions, and use to compute x^4 |0>.  This is something that I can do by hand easily, but it initially looked like I had to compute x^8 |0> for ps8 p1, which I didn’t want to do by hand.  It turned out that we only required x^4|0> for the problem.

Methods.

ClearAll[raise, lower, raiseKet, lowerKet, k, l, k1, x0, xOp, normalize, sumByList , display, formatKet, formatKetList, isNonZero, hbar, omega, m, sopts]

(* Given a ket |k> = {v,n} representing a single basis state v|n> apply the lowering operator: 
a|k> = v(Sqrt[n]) |n-1>.  Handling of |0> is special cased. *) 
lower[{v_,n_Integer}] := {Sqrt[n]v, n-1} ;
lower[{v_,0}] := {0, 0} ;

(* Given a ket |k> = {v,n} representing a single basis state v|n> apply the raising operator: 
(SuperDagger[a])|k> = v(Sqrt[n + 1]) |n+1>. *) 
raise[{v_,n_Integer}] := {Sqrt[n+1]v, n+1} ;

(* Given a ket |k> = (∑_i)(Subscript[v, i]) |(Subscript[n, i])>, represented as a list of kets{{(Subscript[v, 1]), (Subscript[n, 1])},{(Subscript[v, 2]), (Subscript[n, 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_] := (raise[#] &/@ k) ;
lowerKet[k_] := (lower[#] &/@ k) ;

(*Given a list
l = {{v_1,k_1}, {v_2, k_2}, ...} 
find and sum the v_i values for any k_j == k
*)
sumByList[l_,k_] := #[[1]] &/@ Select[l, #[[2]] == k &] // Total ;

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

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

(*Grid[{{ "k = ", k}, 
{"scale = ", scale}, 
{"nonzero = ", nonzero},
{"basisElements = ", basisElements},
{"basisElementMultipliers = ", basisElementMultipliers},
{"result = ", result}
} ] // Print ;*)
result
] ;


(*A display method for |k> = {v,n} *)
display  = {x0 -> (Subscript[x, 0]), hbar -> ℏ, omega -> ω} ;
sopts = {14, "Text"} ;
formatKet[{v_,n_Integer}] := Row[{v /. display, "\[VerticalSeparator]", n,"⟩"}] ;
formatKetList[k_] := (formatKet[#] &/@ k)  // Total;


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

On[Assert] ;

(* 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, #[[2]]==k &]  ;
c2 = Select[k2, #[[2]]==k &]  ;
1 == Length[c1] // Assert ;
1 == Length[c2] // Assert ;

(*Grid[{
{"dotOne:",""},
{"k1 = ", k1}, {"k2 = ", k2},
{"c1 = ", c1}, {"c2 = ", c2},
{"k = ", k}
} ] // Print ;*)

(c1[[1]][[1]])  (c2[[1]][[1]])
] ;

(*compute the inner product of two kets *)
innerProd[k1_, k2_] := Module[{t1, t2, k, dotlist, r, tmp},
t1 = (#[[2]] &/@ k1) ;
t2 = (#[[2]] &/@ 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}
} ] // Print ;*)
r
] ;

Calculation of some (x^4) operators.

ClearAll[ketZero, xp, n, yp, ketOne]
ketZero = {{1,0}}  ;
n = 8 ;
xp = Range[n+1] ;
xp[[1]] = ketZero ;
For[i=1,i ≤ n,i++, xp[[i + 1]] = xp[[i ]] // xOp // Simplify]

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

ketOne = {{1,1}}  ;
yp = Range[n+1] ;
yp[[1]] = ketOne;
For[i=1,i ≤ n,i++, yp[[i + 1]] = yp[[i ]] // xOp // Simplify]

Style[Row[{(x^(#-1)), "|1> = ", yp[[#]] // formatKetList}] & /@ Range[n+1] // TableForm, sopts]
(*innerProd[ xp[[#]], xp[[#]] ] & /@  Range[5] /. display
innerProd[ xp[[#]], xp[[#+2]] ] & /@  Range[3] /. display*)