(* Convert complex variables to real variables so differentiation works correctly. *)
ConvertZtoX[f_]:=f/.{(Subscript[z, 1])->(Subscript[x, 1]),(Subscript[z, 2])->(Subscript[x, 2]),Conjugate[(Subscript[z, 1])]->(Subscript[x, 3]),Conjugate[(Subscript[z, 2])]->(Subscript[x, 4]) }
(* Evaluate (D^α)((OverBar[D])^β)f, where f is in 4 real variables as above. Used to find the basis for (Subscript[H, p,q])((ℂ^2)). *)
MyDiffOp[a_,b_,f_]:=If[Length[a]==Length[b]==2,D[f,{(Subscript[x, 1]),a[[1]]},{(Subscript[x, 2]),a[[2]]},{(Subscript[x, 3]),b[[1]]},{(Subscript[x, 4]),b[[2]]}],"a and b must be vectors in N^2"]
(* Convert real variables to complex variables so we can display the answer properly. *)
ConvertXtoZ[f_]:=f/. {(Subscript[x, 1])->(Subscript[z, 1]),(Subscript[x, 2])->(Subscript[z, 2]),(Subscript[x, 3])-> Conjugate[(Subscript[z, 1])],(Subscript[x, 4])->Conjugate[(Subscript[z, 2])]}
(* Evaluate (D^α)((OverBar[D])^β)f, where f is in 2 complex variables. Used to find the basis for (Subscript[H, p,q])((ℂ^2)). *)
MyDO[a_,b_,f_]:=ConvertXtoZ[MyDiffOp[a,b,ConvertZtoX[f]]]
(* Take the Kelvin transform of a function in 2 complex variables. Done using a simple substitution. *)
KelvinTrans[f_]:=((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar])^(-2))f/.{(Subscript[z, 1])->(Subscript[z, 1])/((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar])^2),(Subscript[z, 2])->(Subscript[z, 2])/((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar])^2)}
(* Given multi-indices a & b, compute K[((OverBar[D])^α)(D^β)f]. Used to find a basis element of (Subscript[H, p,q])((ℂ^2)). *)
BasisElement[a_,b_]:=KelvinTrans[MyDO[b,a,ModZExpand[((\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar])^(-2))]]]
(* Expand the symbol \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar], as it is unevaluated by Mathematica. *)
ModZExpand[f_]:=f/.\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]->((((Subscript[z, 1])Conjugate[(Subscript[z, 1])]+(Subscript[z, 2])Conjugate[(Subscript[z, 2])]))^(1/2))
(* Simplify the symbol \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar], as it is unevaluated by Mathematica. *)
ModZSimplify[f_]:=f/.Conjugate[(Subscript[z, 1])] (Subscript[z, 1])+Conjugate[(Subscript[z, 2])] (Subscript[z, 2])-> \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]^2 /. Conjugate[\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]]->\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]
(* Find a basis for (Subscript[H, p,q])((ℂ^2)). *)
Basis[p_,q_]:=Join[Table[BasisElement[{0,p},{i,q-i}],{i,0,q}],Table[BasisElement[{j,p-j},{0,q}],{j,1,p}]]
(* Find a basis for (Subscript[H, p,q])((S^3)). *)
SBasis[p_,q_]:=Basis[p,q]/.{\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]-> 1 ,Conjugate[(Subscript[z, 1])] (Subscript[z, 1])+Conjugate[(Subscript[z, 2])] (Subscript[z, 2])->1,Conjugate[(Subscript[z, 2])] (Subscript[z, 2])+Conjugate[(Subscript[z, 1])] (Subscript[z, 1])->1}
(* Simplify the basis provided by Basis[p, q] as it is quite messy by default. Most of the messiness is because Mathematica doesn't know how to simplify \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar], so we help it do that. *)
SimplifiedBasis[p_, q_] := Expand[ModZSimplify[Simplify[ModZSimplify[Basis[p,q]]]]]
(* Find a basis for (Subscript[H, k])((ℂ^2)) by joining the bases for (Subscript[H, p,q])((ℂ^2)) together. *)
HBasis[k_] := Flatten[Table[SimplifiedBasis[i, k-i], {i, 0, k}]]
(* Define the ℒ, (OverBar[ℒ]), and (Subscript[□, b]) operators. We need to use the ZtoX and XtoZ conversions above for the differentiation to evaluate properly. *)
LOperator[f_]:=ConvertXtoZ[(Subscript[x, 4])*D[ConvertZtoX[f],(Subscript[x, 1])]-(Subscript[x, 3])*D[ConvertZtoX[f],(Subscript[x, 2])]]
LBarOperator[f_]:= ConvertXtoZ[(Subscript[x, 2])*D[ConvertZtoX[f],(Subscript[x, 3])]-(Subscript[x, 1])*D[ConvertZtoX[f],(Subscript[x, 4])]]
BoxB[f_]:=-1*LOperator[LBarOperator[f]]
(* Define the normalized (Subsuperscript[□, b, t]) operator. We don't use it at all because the un-normalized one is easier to work with. *)
normalizer = (1 + t Conjugate[t])/(((1 - t Conjugate[t]))^2)
BoxBT[f_] := -normalizer*(LOperator[LBarOperator[f] + t LOperator[f]] + Conjugate[t] LBarOperator[LBarOperator[f] + t LOperator[f]])
(* Define the un-normalized (Subsuperscript[□, b, t]) operator. *)
BoxBTU[f_] := -1*(LOperator[LBarOperator[f] + t LOperator[f]] + Conjugate[t] LBarOperator[LBarOperator[f] + t LOperator[f]])
(* Define the inner product on (L^2)((S^3)) by first parametrizing the 3-sphere as below in 3 variables, and then integrating ∫p(OverBar[q]) dσ under this parametrization. *)
Parametrize[p_] := p /. {(Subscript[z, 1]) -> r*(E^(I*θ)), (Subscript[z, 2]) -> Sqrt[1-(r^2)]*(E^(I*ϕ)), Conjugate[(Subscript[z, 1])] -> r*(E^(-I*θ)), Conjugate[(Subscript[z, 2])] -> Sqrt[1-(r^2)]*(E^(-I*ϕ))}
InnerProd[p_, q_] := Integrate[Parametrize[p] Parametrize[Conjugate[q]]r(((1-(r^2)))^(1/2)),{r,0,1},{θ,0,2π},{ϕ,0,2π}]
(* Find the matrix representation of O under the basis b by taking (L^2) inner products. Not used anymore due to being extremely slow to compute. *)
BuildMatrix[O_, b_] := (lengths = InnerProd[#, #]& /@ b; ParallelTable[Table[InnerProd[O[ b[[i]] ], b[[j]] ]/lengths[[j]], {j, 1, Length[b]}], {i, 1, Length[b]}])
(* Shorthand. *)
(Subscript[H, k_]) := HBasis[k]
L2[f_] := LOperator[LOperator[f]]
LBar2[f_] := LBarOperator[LBarOperator[f]]
(* Find the matrix representation of (Subsuperscript[□, b, t]) more efficiently, since we know there is one entry per row/column if we expand the definition. *)
SingleEntryMatrix[O_, b_] := Transpose[ParallelTable[Table[(entry=Simplify[Expand[O[b[[i]]]/b[[j]]]]; If[NumberQ[entry],entry,0]),{j,1,Length[b]}], {i, 1, Length[b]}]]
BoxBTUMatrix[b_] := -SingleEntryMatrix[LOperator @* LBarOperator,b] - t*Conjugate[t]*SingleEntryMatrix[LBarOperator @* LOperator,b]-t*SingleEntryMatrix[L2, b]-Conjugate[t]*SingleEntryMatrix[LBar2, b]
(* Compute the basis of (Subscript[H, p,q])((ℂ^2)) and (Subscript[H, k])((ℂ^2)), but append {(Subscript[α, 1]), (Subscript[α, 2]), (Subscript[β, 1]), (Subscript[β, 2])} to each basis element. We used it to explore and eventually find the (Subscript[V, i]) and (Subscript[W, i]) spaces, along with their definitions. *)
SimplifiedBasisElement[a_, b_] := Expand[ModZSimplify[Simplify[ModZSimplify[BasisElement[a,b]]]]]
SpecialBasis[p_, q_] := Join[Table[{SimplifiedBasisElement[{0,p},{i,q-i}],{0,p,i,q-i}},{i,0,q}], Table[{SimplifiedBasisElement[{j,p-j},{0,q}], {j, p-j, 0, q}},{j,1,p}]]
SpecialHBasis[k_] := Flatten[Table[SpecialBasis[i, k-i], {i, 0, k}],1]
(* Do SingleEntryMatrix and BlockDiagonalize but with a special basis instead of a normal basis. *)
SpecialSingleEntryMatrix[O_, sb_] :=Transpose[ParallelTable[Table[(entry=Simplify[Expand[O[sb[[i,1]]]/sb[[j,1]]]]; If[NumberQ[entry],entry,0]),{j,1,Length[sb]}], {i, 1, Length[sb]}]]
SpecialBlockDiagonalize[mat_, sb_] := 
(result = mat;
newsb = sb;
Do[
(nonzerorow = Position[result[[i,i+1;;]], x_/;x=!=0,1];
If[ Length[nonzerorow] ≤ 1,Nothing,(result = MoveMatrix[result, nonzerorow[[2,1]]+i,i+1]; newsb = MoveElem[newsb,nonzerorow[[2,1]]+i,i+1])]),
{i,Length[mat]-1}];
{result, newsb})
(* Move the element of the list l_ at position from_ to position to_. *)
MoveElem[l_,from_,to_]:= Insert[Delete[l,from],l[[from]],to] 
(* Swap rows from_ and to_, as well as columns from_ and to_ in the matrix mat_. This is equivalent to changing the order of which the basis elements appear. *)
MoveMatrix[mat_, from_, to_] := MoveElem[MoveElem[#, from, to]& /@mat, from, to]
(* Turn the normal matrix reperesentation of (Subsuperscript[□, b, t]) into a block diagonal matrix, where each block is tridiagonal. *)
BlockDiagonalize[mat_] := 
(result = mat;
Do[
(nonzerorow = Position[result[[i,i+1;;]], x_/;x=!=0,1];
If[ Length[nonzerorow] ≤ 1,Nothing,result = MoveMatrix[result, nonzerorow[[2,1]]+i,i+1]]),
{i,Length[mat]-1}];
result)
(* Determine how an operator O_ permutes the elements of the special basis sb_. 
Each arrow in the output {(Subscript[α, 1]), (Subscript[α, 2]), (Subscript[β, 1]), (Subscript[β, 2])} -> {(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[b, 1]), (Subscript[b, 2])} means O[((OverBar[(Subscript[D, 1])])^(Subscript[α, 1]))((OverBar[(Subscript[D, 2])])^(Subscript[α, 2]))((Subscript[D, 1])^(Subscript[β, 1]))((Subscript[D, 2])^(Subscript[β, 2]))|z((|)^(-2))] = λ((OverBar[(Subscript[D, 1])])^(Subscript[a, 1]))((OverBar[(Subscript[D, 2])])^(Subscript[a, 2]))((Subscript[D, 1])^(Subscript[b, 1]))((Subscript[D, 2])^(Subscript[b, 2]))|z((|)^(-2)).
Used to investigate the behavior of (ℒ^2) and ((OverBar[ℒ])^2) in the expansion of (Subsuperscript[□, b, t]). *)
HMappings[O_, sb_] :=
(ofuncs = O /@ sb[[;;,1]];
Flatten[Table[
If[ofuncs[[i]] === 0,
 sb[[i,2]] -> 0,
(*Else*)
 Table[(entry=Simplify[Expand[ ofuncs[[i]]/sb[[j,1]] ]]; If[IntegerQ[entry],sb[[i,2]] -> sb[[j,2]],Nothing]),
{j,1,Length[sb]}]],
{i, 1, Length[sb]}], 1])

(* Find the continuant {(Subscript[f, 0]), (Subscript[f, 1]), ..., (Subscript[f, n])} of a tridiagonal matrix mat_. *)
Continuant[mat_] := (list = {1, mat[[1,1]]}; 
Do[list = Append[ list, Expand[list[[i]]*mat[[i,i]] - list[[i-1]]*mat[[i-1,i]]*mat[[i,i-1]]] ],{i, 2, Length[mat]}];
	list)

EigenvalPlot[f_] := Plot[f, {x, 0, 1}, AxesOrigin->{0,0}]

(* Find the (Subscript[V, i]) and (Subscript[W, i]) spaces, but they were not named this at the time. *)
BasisChain[p_, q_] := (basis = SimplifiedBasis[p,q];  Expand /@Table[NestList[Expand @* LBar2, basis[[i]], (p+q-1)/2], {i, Length[basis]}])
BasisChains[k_] := Join[BasisChain[0,k], BasisChain[1,k-1]]

(* Miscellaneous functions used to investigate why ℒ takes basis elements to basis elements. We were unable to prove this, and this helped illuminate why it was so difficult to prove. *)
LPart1[f_]:=ConvertXtoZ[(Subscript[x, 4])*D[ConvertZtoX[f],(Subscript[x, 1])]]
LPart2[f_] := ConvertXtoZ[-(Subscript[x, 3])*D[ConvertZtoX[f],(Subscript[x, 2])]]
Nonzero[x_] := If[x === 0, 0, 1]

(* Find the matrix representations for (Subsuperscript[□, b, t]) on the V and W spaces on (Subscript[H, 2k-1])((ℂ^2)) using the formulas in the paper. *)

OddVMatrix[k_] := If[k ≠ 1,DiagonalMatrix[Table[-(2k-2i)(2k+1-2i)2i(2i-1)t, {i, 1, k-1}],1] + DiagonalMatrix[Table[(2i-1)(2k+1-2i) + (2i-2)(2k+2-2i)t Conjugate[t], {i, 1, k}]] + DiagonalMatrix[Table[-Conjugate[t], {i, 1, k-1}],-1], DiagonalMatrix[{1, t Conjugate[t]}]]
OddWMatrix[k_]  := If[k≠ 1,DiagonalMatrix[Table[-(2k-2i)(2k-1-2i)2i(2i+1)t, {i, 1, k-1}],1] + DiagonalMatrix[Table[(2i)(2k-2i) + (2i-1)(2k+1-2i)t Conjugate[t], {i, 1, k}]] + DiagonalMatrix[Table[-Conjugate[t], {i, 1, k-1}],-1], DiagonalMatrix[{1, t Conjugate[t]}]]

(* Make a tridiagonal matrix a symmetric tridiagonal matrix with the similarity relationship discussed in the paper. *)

MakeSymmetric[mat_] := DiagonalMatrix[Diagonal[mat]]  + DiagonalMatrix[Table[Sqrt[mat[[i,i+1]]*mat[[i+1,i]]], {i, Length[mat]-1}],1] + DiagonalMatrix[Table[Sqrt[mat[[i,i+1]]*mat[[i+1,i]]], {i, 1, Length[mat]-1}],-1]

(* Find the eigenvalues of (Subsuperscript[□, b, t]) on V and W spaces. Since the matrices are much smaller and don't have repeated eigenvalues, this is much quicker to compute. *)

EigenvalOddW[k_] :=(eigenvals = DeleteDuplicates[Eigenvalues[OddWMatrix[k]]]; plottables = # /. {t*Conjugate[t] -> (x^2), (t^(p_))((Conjugate[t])^(p_))-> (x^(2*p))}& /@ eigenvals)
EigenvalOddV[k_] :=(eigenvals = DeleteDuplicates[Eigenvalues[OddVMatrix[k]]]; plottables = # /. {t*Conjugate[t] -> (x^2), (t^(p_))((Conjugate[t])^(p_))-> (x^(2*p))}& /@ eigenvals)

bounds = TtoXReplace[Table[Det[OddWMatrix[i]]/Det[OddWMatrix[i][[;;-2,;;-2]]], {i, 1, 15}]]

bounds /. x-> 1

Plot[bounds, {x, 0, 1}, PlotRange->Full, PlotLegends->Automatic]

TtoXReplace[list_] := list /. {Conjugate[t] -> x,t -> x}

maybebound[k_, t_] := (2k-1)(t^(2k))(E^(Sum[1/(2i), {i, 1, k-1}]))