Rank3CentralFunctions.nb

A notebook for calculating central functions of SL(2,C) character varieties of rank 3 free groups.

Sean Lawton

July 8, 2008

Getting Started

We assume that the reader has read or has access to Spin Networks and SL(2,C)-Character Varieties (available on the arXiv).  So given that, the reader then want to compute some central functions.  First, under the Evaluation tab above, select  ``Evaluate Initialization Cells``.

For rank 1 central functions use the command: Central[a], and for rank 2 central functions use Central[a,b,c]:

Central[1]
Central[2]
Central[3]

Central[1,0,1]
Central[1,1,0]

For rank 3 cenral function, with respect to the preprint, Computing with SL(2,C) Central Functions, use the command Central[a,b,c,d,i,j]:

Central[1,1,1,1,1,1]

At the bottom we compute many further examples of rank 3 central functions.

Here is the code to make these computations, with respect to tensiorial contraction (it is only partially documented):

Preamble

A sample matrix for computing examples:

Mm={{a,b},{c,d}};

Here are the heads for monomials in the symmetric algebra:

Protect[pp,PP,qq,QQ,rr,RR,cc,CC];
HHeadQ[x_] := MemberQ[SymHeads={pp,qq,QQ,PP,rr,RR,cc,CC},x]

Since we are working with a free group of rank 3 we need three heads for the three tensors: pp, qq, and rr.  We have three more for their duals: PP, QQ, and RR.  The heads cc and CC are intermediate heads that are used for the first part of a two part injection that will be described later.

$IterationLimit=Infinity;
$RecursionLimit=Infinity;

The Goldman slice

Here are the three elements of SL(2,C) which realize the character (t1,t2,t12,t3,t23,t13). Here t12=  w+1/w.

T1 = {{t1,-1},{1,0}}; T2 = {{0,w},{-1/w,t2}};
T3[s_]:={{t3,(((t13-t1*t3)*w+t23)*w)/(w^2-1)},{((t13-t1*t3)+t23*w)/(w^2-1),0}}+s*(T1.T2-T2.T1);

There is a complex 1 dimensional family of solutions for the value of s that satisfy the trace conditions over the Lie algebra gl(2,C).  However, over the Lie group SL(2,C) there are two solutions give below (to view then remove the semi-colon):

Solve[Det[T3[s]]==1,s];

We will find it useful later to use one of these solutions:

ssub := Solve[Det[T3[s]] == 1, s][[1]]

Checking that these do in fact give a slice we first define the three matrix generators and their products:

Map[MatrixForm,GeneratorList={T1, T2,T3[s], T12 = T1.T2,T13=T1.T3[s],T23=T2.T3[s]}]

These three matrices have the desired traces:

Map[Tr,GeneratorList]//.{w+1/w->t12}

Simplify[%]

However, although locally these six traces determine a representation, globally they do not since the character variety is a hypersurface in (C^7) that is monic quadratic in an additional variable over the other six.  The triple products give the desired extention; that is, are the roots.

Expand[FullSimplify[Tr[T1.T3[s].T2]/.ssub]]
Expand[FullSimplify[Tr[T1.T2.T3[s]]/.ssub]]

In both we notice a squareroot.  It will come up in calculations later on, so we define the following substitution:

quadsub:={(√(-4+(t3^2) ((t1^2) (t2^2) (w^2)+(((-1+(w^2)))^2)-2 t1 t2 w (1+(w^2)))+w ((t13^2) (-4+(t2^2)) w+(t1^2) (-4+(t23^2)) w-4 t13 t23 (1+(w^2))-4 w (-2+(t2^2)+(t23^2)+(w^2))+2 t1 t2 (2+t13 t23 w+2 (w^2)))+2 t3 w (t23 (t1-(-2+(t1^2)) t2 w+t1 (w^2))+t13 (t2-t1 (-2+(t2^2)) w+t2 (w^2)))))-> (-2*w)*(t123-(1/2)*(t2*t13+t1*t23-t1*t2*t3+t3*t12))}

The two roots are related by the following substitution:

finalsub:={t132->t12*t3+t13*t2+t23*t1-t1*t2*t3-t123}

Bilinear and linear mappings

One basic operation is tensor product (⊗).  Bilinearity is its defining property:

Make it bilinear:

TensorProduct[X_ + Y_,Z_ ] := TensorProduct[X,Z] + TensorProduct[Y,Z]
TensorProduct[X_,Y_ + Z_] := TensorProduct[X,Y] + TensorProduct[X,Z]
TensorProduct[s_Integer  X_,Y_] := s TensorProduct[X,Y]
TensorProduct[s_Rational  X_,Y_] := s TensorProduct[X,Y]
TensorProduct[X_,s_Integer Y_] := s TensorProduct[X,Y]
TensorProduct[X_,s_Rational Y_] := s TensorProduct[X,Y]

TensorProduct[X_ + Y_,Z_ ,W_] := TensorProduct[X,Z,W] + TensorProduct[Y,Z,W]
TensorProduct[X_,Y_ + Z_,W_] := TensorProduct[X,Y,W] + TensorProduct[X,Z,W]
TensorProduct[X_,Y_ , Z_+W_] := TensorProduct[X,Y,Z] + TensorProduct[X,Y,W]
TensorProduct[s_Integer  X_,Y_,Z_] := s TensorProduct[X,Y,Z]
TensorProduct[s_Rational  X_,Y_,Z_] := s TensorProduct[X,Y,Z]
TensorProduct[X_,s_Integer Y_,Z_] := s TensorProduct[X,Y,Z]
TensorProduct[X_,s_Rational Y_,Z_] := s TensorProduct[X,Y,Z]
TensorProduct[X_,Y_,s_Integer Z_] := s TensorProduct[X,Y,Z]
TensorProduct[X_,Y_,s_Rational Z_] := s TensorProduct[X,Y,Z]

Use this function to display tensor products using ⊗:

DisplayT[x_]:=x/.(TensorProduct->CircleTimes)

The bilinear map ScalarProduct implements duality:

ScalarProduct[X_ + Y_,Z_ ] := ScalarProduct[X,Z] + ScalarProduct[Y,Z]
ScalarProduct[X_,Y_ + Z_] := ScalarProduct[X,Y] + ScalarProduct[X,Z]
ScalarProduct[s_Rational  X_,Y_] := s ScalarProduct[X,Y]
ScalarProduct[X_,s_Rational Y_] := s ScalarProduct[X,Y]
ScalarProduct[s_Integer  X_,Y_] := s ScalarProduct[X,Y]
ScalarProduct[X_,s_Integer Y_] := s ScalarProduct[X,Y]

ScalarProduct[TensorProduct[a1_,b1_],TensorProduct[a2_,b2_]]:=ScalarProduct[a1,a2]ScalarProduct[b1,b2]

ScalarProduct[TensorProduct[a1_,b1_,c1_],TensorProduct[a2_,b2_,c2_]]:=ScalarProduct[a1,a2]ScalarProduct[b1,b2]*ScalarProduct[c1,c2]

ScalarProduct[s_  f_[a___],Y_] := s ScalarProduct[f[a],Y]/;HHeadQ[f]
ScalarProduct[X_,s_ f_[a___]] := s ScalarProduct[X,f[a]]/;HHeadQ[f]

ScalarProduct implements the pairing between the symmetric power of a vector space and the symmetric power
of its dual. The (i^th) basic monomial (Subscript[ϕ, i]) in the (m^th) symmetric power (Subscript[V, m]) is the symmetrization of the (i^th) basic monomial in the (m^th) tensor power (⊗^m)V. In Mathematica we denote this element by pp[i,m].

Similarly, the dual ((SuperStar[(⊗^m)V)]) canonically identifies with (⊗^m)((SuperStar[V])), we obtain a corresponding basis Φ0, ..., (Subscript[Φ, m ])of (⊗^m)((SuperStar[V]))
obtained by symmetrizing basic monomials in (⊗^m)((SuperStar[V])) (denoted in Mathematica by PP[0,m], ..., PP[m,m]). However these bases are not dual,  due to the symmetrizations: The scalar product of (Subscript[ϕ, i])with (Subscript[Φ, j ])equals  (({{m}, {i}}))(Subscript[δ, ij]).

ScalarProduct[pp[i_,m_],PP[j_,n_]]:=If[i==j&&m==n,1/Binomial[m,i],0]
ScalarProduct[qq[i_,m_],QQ[j_,n_]]:=If[i==j&&m==n,1/Binomial[m,i],0]
ScalarProduct[rr[i_,m_],RR[j_,n_]]:=If[i==j&&m==n,1/Binomial[m,i],0]

GL(2) induces linear maps on the symmetric powers:

InducedMap[Mat_,X_ + Y_ ] := InducedMap[Mat,X ]+ InducedMap[Mat, Y]
InducedMap[Mat_,s_Integer X_] := s InducedMap[Mat,X]
InducedMap[Mat_,s_Rational X_] := s InducedMap[Mat,X]

Internal tensor multiplication

The next function implement  internal multiplication of symmetric tensors:

TensorTimes[TensorProduct[a_[i1_,m1_],b_[j1_,n1_]],TensorProduct[a_[i2_,m2_],b_[j2_,n2_]]]:=
TensorProduct[a[i1+i2,m1+m2],b[j1+j2,n1+n2]]/;HHeadQ[a]&&HHeadQ[b]

TensorTimes[TensorProduct[a_[i1_,m1_],b_[j1_,n1_]],c_[k1_,o1_]]:=
TensorProduct[a[i1,m1],b[j1,n1],c[k1,o1]]/;HHeadQ[a]&&HHeadQ[b]&&HHeadQ[c]

Make it bilinear:

TensorTimes[X_+Y_,Z_]:=TensorTimes[X,Z]+TensorTimes[Y,Z]
TensorTimes[X_,Y_+Z_] :=TensorTimes[X,Y]+TensorTimes[X,Z]
TensorTimes[s_ X_, Y_] := s TensorTimes[X,Y]
TensorTimes[X_,s_  Y_] := s TensorTimes[X,Y]

Later on, we will want the following function to be linear:

EvalMatCoeff[X_+Y_,{pMatrix_,qMatrix_,rMatrix_}]:=EvalMatCoeff[X,{pMatrix,qMatrix,rMatrix}]+EvalMatCoeff[Y,{pMatrix,qMatrix,rMatrix}]
EvalMatCoeff[s_Integer  X_,{pMatrix_,qMatrix_,rMatrix_}]:= s EvalMatCoeff[X,{pMatrix,qMatrix,rMatrix}]
EvalMatCoeff[s_Rational X_,{pMatrix_,qMatrix_,rMatrix_}]:=s EvalMatCoeff[X,{pMatrix,qMatrix,rMatrix}]

Symmetric powers

The function SymmetricPower returns the action of an element of GL(2) on the (m^th) symmetric power (Subscript[V, m]). First we compute the induced map on the basis pp[i,m] of (Subscript[V, m]), and then we use the basis PP[i,m] of (Subscript[((SuperStar[V])), m]) to compute matrix coefficients.

The next function provides a list corresponding to the basis of the symmetric monomials:

LList[f_,m_]:=Map[f[#,m]&,Range[0,m]]

The function supplied by Mathematica for obtaining the coefficients of a polynomial
is inadequate if the polynomial has missing terms:

NewCoefficientList[f_,xxx_,n_]:= Map[Coefficient[xxx f,xxx^#]&,Range[1,n+1]]

The next function computes the map induced on a symmetric power by a 2x2 matrix Mat:

InducedMap[Mat_,f_[i_,m_]]:=Inner[Times,NewCoefficientList[Expand[((Mat.{1,xx})[[1]])^(m-i)((Mat.{1,xx})[[2]])^i],xx,m],Map[f[#,m]&,Range[0,m]]]

The next function gives the symmetric power in terms of the bases of (Subscript[V, m])  and (Subscript[((SuperStar[V])), m]).

SymmetricPower[Mat_,m_]:=Outer[ScalarProduct,Map[InducedMap[Mat,#]&,LList[pp,m]],LList[PP,m]]

Map[MatrixForm[SymmetricPower[Mm,#]]&,Range[0,3]]

Recall that the (((i,j))^th) entry (Subscript[M, i j ]) = M[[i,j]]  ((i^(th ))row, (j^(th )) column)  of a matrix M equals the (i^(th ))basic covector evaluated on the image of the (j^(th )) basis vector under M.

The Clebsch-Gordan Summand

Suppose integers m, n, and r satisfy  0 ≤ m, n, r. Then the representation (Subscript[V, r]) includes into (Subscript[V, m])⊗ (Subscript[V, n]) if and only if  |m-n| ≤  r ≤  m+n and  m+n-r even. In that case the inclusion is
unique, and the Clebsch-Gordan summand (Subscript[V, r]^(m,n)) is the image of the inclusion
of  the representation (Subscript[V, r]) into (Subscript[V, m])⊗ (Subscript[V, n]).

The following tests whether a summand exists.

ClebschGordanQ[{m_,n_,r_}]:=IntegerQ[m]&&IntegerQ[n]&&IntegerQ[r]&&m≥0&&n≥0&&r≥0&& IntegerQ[(m+n-r)/2]&& (r≤ m+n) && Abs[m-n]≤r

Rank3AdmQ[{a_,b_,c_,d_}]:=IntegerQ[a]&&IntegerQ[b]&&IntegerQ[c]&&IntegerQ[d]&&a≥0&&b≥0&&c≥0&&d≥0&& MemberQ[ListRank3Adm[a,b,c],d]

ListRank2Adm[a_,b_]:=ListRank2Adm[a,b]=Module[{dummy={}},For[i=0,i<Min[a,b]+1,i++,
dummy=Join[dummy,{a+b-2*i}]]; dummy];
ListRank3Adm[a_,b_,c_]:=ListRank3Adm[a,b,c]=Module[{dummy={}},For[ii=1,ii<Length[ListRank2Adm[a,b]]+1,ii++,dummy=Join[dummy,ListRank2Adm[ListRank2Adm[a,b][[ii]],c]]]; dummy];

Denote the basis of (Subscript[V, m]) by (Subscript[φ, 0]) , ..., (Subscript[φ, m]) , the basis of (Subscript[V, n]) by (Subscript[ψ, 0, ])...,(Subscript[ψ, n])   and the basis of (Subscript[V, r]) by 
(Subscript[τ, 0, ]) ..., (Subscript[τ, r]).   If r = m + n,  the inclusion (Subscript[V, r]) into (Subscript[V, m])⊗ (Subscript[V, n]) is given by

                                        (({{m+n}, {k}}))(Subscript[τ, k       ])↦  (∑_( 0≤ i≤m
 0≤j≤n
 i+j=k))    (({{m}, {i}}))(Subscript[φ, i])  ⊗   (({{n}, {j}}))(Subscript[ψ, j ])

In general the inclusion (Subscript[V, r]) into (Subscript[V, m])⊗ (Subscript[V, n]) is given by
 (({{r}, {k}}))(Subscript[τ, k ])↦  (  (∑_( 0≤ i≤m-s
 0≤j≤n-s
 i+j=k))   (({{m-s}, {i}})) (Subscript[φ, i])  ⊗   (({{n-s}, {j}}))(Subscript[ψ, j ])) ((( (Subscript[φ, 0])⊗ (Subscript[ψ, 1 ])-  (Subscript[φ, 1])⊗ (Subscript[ψ, o ])))^s)
      where 0 ≤ k ≤ r,  s = (m + n - r) /2 and multiplication is internal mutiplication of tensors.

The following is a long list of functions which iteratively composes two of these injections to yield an injection into a triple tensor.

alpha[a_,b_,c_]:=alpha[a,b,c]=(1/2)*(b+c-a);
beta[a_,b_,c_]:=beta[a,b,c]=(1/2)*(a-b+c);
gamma[a_,b_,c_]:=gamma[a,b,c]=(1/2)*(a+b-c);
delta[a_,b_,c_]:=delta[a,b,c]=(1/2)*(a+b+c);

alpha, beta, gamma, and delta are the terms in which the coefficients in the Clebsch-Gordan injection given above.

IncludeSummandRank3[k_][{a_,b_,c_,d_},i_]:=IncludeSummandRank3[k][{a,b,c,d},i,{pp,qq,rr}]/;Rank3AdmQ[{a,b,c,d}]&& 0≤k&&k≤d&&1≤i&&i≤ Mult[a,b,c][d]
IncludeSummandRank3[{a_,b_,c_,d_},i_]:=IncludeSummandRank3[{a,b,c,d},i]=Map[IncludeSummandRank3[#][{a,b,c,d},i]&,Range[0,d]]

IncludeDualSummandRank3List[k_][{a_,b_,c_,d_},i_]:=IncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}]/;Rank3AdmQ[{a,b,c,d}]&& 0≤k&&k≤d&&1≤i&&i≤ Mult[a,b,c][d]
IncludeDualSummandRank3List[{a_,b_,c_,d_},i_]:=IncludeDualSummandRank3List[{a,b,c,d},i]=Map[IncludeDualSummandRank3List[#][{a,b,c,d},i]&,Range[0,d]]

PreIncludeDualSummandRank3List[k_][{a_,b_,c_,d_},i_]:=PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}]/;Rank3AdmQ[{a,b,c,d}]&& 0≤k&&k≤d&&1≤i&&i≤ Mult[a,b,c][d]
PreIncludeDualSummandRank3List[{a_,b_,c_,d_},i_]:=PreIncludeDualSummandRank3List[{a,b,c,d},i]=Map[PreIncludeDualSummandRank3List[#][{a,b,c,d},i]&,Range[0,d]]

PreIncludeSummandRank3List[k_][{a_,b_,c_,d_},i_]:=PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}]/;Rank3AdmQ[{a,b,c,d}]&& 0≤k&&k≤d&&1≤i&&i≤ Mult[a,b,c][d]
IncludeSummandRank3[{a_,b_,c_,d_},i_]:=IncludeSummandRank3[{a,b,c,d},i]=Map[IncludeSummandRank3[#][{a,b,c,d},i]&,Range[0,d]]

IncludeSummandRank3List[{a_,b_,c_,d_},i_]:=IncludeSummandRank3List[{a,b,c,d},i]=Map[IncludeSummandRank3List[#][{a,b,c,d},i]&,Range[0,d]]

PreIncludeSummandRank3List[{a_,b_,c_,d_},i_]:=PreIncludeSummandRank3List[{a,b,c,d},i]=Map[PreIncludeSummandRank3List[#][{a,b,c,d},i]&,Range[0,d]]

IncludeSummandRank3List[k_][{a_,b_,c_,d_},i_]:=IncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}]/;Rank3AdmQ[{a,b,c,d}]&& 0≤k&&k≤d&&1≤i&&i≤ Mult[a,b,c][d]

IncludeDualSummandRank3[k_][{a_,b_,c_,d_},i_]:=IncludeDualSummandRank3[k][{a,b,c,d},i,{PP,QQ,RR}]/;Rank3AdmQ[{a,b,c,d}]&& 0≤k&&k≤d&&1≤i&&i≤ Mult[a,b,c][d]
IncludeDualSummandRank3[{a_,b_,c_,d_},i_]:=IncludeDualSummandRank3[{a,b,c,d},i]=Map[IncludeDualSummandRank3[#][{a,b,c,d},i]&,Range[0,d]]

Listijm[{a_,b_,c_},k_]:=Module[{dummy={}},For[i=0,i<beta[a,b,c]+1,i++,For[j=0,j<alpha[a,b,c]+1,j++,For[m=0,m<gamma[a,b,c]+1,m++,If[i+j==k,dummy=Join[dummy,{{i,j,m}}]]]]];Listijm[{a,b,c},k]=dummy]/;ClebschGordanQ[{a,b,c}]&& 0≤k&&k≤c

IncludeSummandRank2List[k_][{a_,b_,c_},{x_,y_}]:=Map[(Power[(-1),Listijm[{a,b,c},k][[#]][[3]]]/Binomial[c,k])*Binomial[beta[a,b,c],Listijm[{a,b,c},k][[#]][[1]]]*Binomial[alpha[a,b,c],Listijm[{a,b,c},k][[#]][[2]]]*Binomial[gamma[a,b,c],Listijm[{a,b,c},k][[#]][[3]]]*TensorProduct[x[gamma[a,b,c]+Listijm[{a,b,c},k][[#]][[1]]-Listijm[{a,b,c},k][[#]][[3]],a],y[Listijm[{a,b,c},k][[#]][[2]]+Listijm[{a,b,c},k][[#]][[3]],b]]&,Range[1,Length[Listijm[{a,b,c},k]]]]/;ClebschGordanQ[{a,b,c}]&&HHeadQ[x]&&HHeadQ[y]&&0≤k&&k≤c

IncludeDualSummandRank2List[k_][{a_,b_,c_}]:=IncludeDualSummandRank2List[k][{a,b,c},{pp,qq}]/;ClebschGordanQ[{a,b,c}]&& 0≤k&&k≤c
IncludeDualSummandRank2List[{a_,b_,c_}]:=IncludeDualSummandRank2List[{a,b,c}]=Map[IncludeDualSummandRank2List[#][{a,b,c}]&,Range[0,c]]

IncludeSummandRank2List[k_][{a_,b_,c_}]:=IncludeSummandRank2List[k][{a,b,c},{pp,qq}]/;ClebschGordanQ[{a,b,c}]&& 0≤k&&k≤c
IncludeSummandRank2List[{a_,b_,c_}]:=IncludeSummandRank2List[{a,b,c}]=Map[IncludeSummandRank2List[#][{a,b,c}]&,Range[0,c]]

IncludeDualSummandRank2List[k_][{a_,b_,c_},{x_,y_}]:=Map[(Power[(-1),Listijm[{a,b,c},k][[#]][[3]]]/Binomial[c,k])*Binomial[beta[a,b,c],Listijm[{a,b,c},k][[#]][[1]]]*Binomial[alpha[a,b,c],Listijm[{a,b,c},k][[#]][[2]]]*Binomial[gamma[a,b,c],Listijm[{a,b,c},k][[#]][[3]]]*TensorProduct[x[gamma[a,b,c]+Listijm[{a,b,c},k][[#]][[1]]-Listijm[{a,b,c},k][[#]][[3]],a],y[Listijm[{a,b,c},k][[#]][[2]]+Listijm[{a,b,c},k][[#]][[3]],b]]&,Range[1,Length[Listijm[{a,b,c},k]]]]/;ClebschGordanQ[{a,b,c}]&&HHeadQ[x]&&HHeadQ[y]&&0≤k&&k≤c

IncludeSummandRank2[k_][{a_,b_,c_},{x_,y_}]:=Apply[Plus,IncludeSummandRank2List[k][{a,b,c},{x,y}]]/;ClebschGordanQ[{a,b,c}]&&HHeadQ[x]&&HHeadQ[y]&&0≤k&&k≤c

IncludeDualSummandRank2[k_][{a_,b_,c_},{x_,y_}]:=Apply[Plus,IncludeDualSummandRank2List[k][{a,b,c},{x,y}]]/;ClebschGordanQ[{a,b,c}]&&HHeadQ[x]&&HHeadQ[y]&&0≤k&&k≤c

PointedListRank3Adm[a_,b_,c_]:=PointedListRank3Adm[a,b,c]=Module[{dummy={}},For[ii=1,ii<Length[ListRank2Adm[a,b]]+1,ii++,dummy=Join[dummy,{{ListRank2Adm[ListRank2Adm[a,b][[ii]],c],ListRank2Adm[a,b][[ii]]}}]]; dummy];

The elements of PointedListRank3Adm are of the form {{Set}, d} associtated to a triple of numbers (a, b, c).  The elements of Set are
 the irreps (dim - 1) in the C - G decomp of  the triple tensor V_aV _bV _c coming from the imbedding of V_dV _c in (V_aV _b) V_c.

Mult[a_,b_,c_][d_]:=0/;Rank3AdmQ[{a,b,c,d}]==False
Mult[a_,b_,c_][d_]:=Mult[a,b,c][d]=Module[{counter=0},For[j=1,j<Length[PointedListRank3Adm[a,b,c]]+1,j++,If[MemberQ[PointedListRank3Adm[a,b,c][[j]][[1]],d],counter=counter+1]];counter]

depth[d_][a_,b_,c_]:=Module[{depthdummy={}},For[idepth=1,idepth<Length[PointedListRank3Adm[a,b,c]]+1,idepth++,If[MemberQ[ PointedListRank3Adm[a,b,c][[idepth]][[1]],d],depthdummy=Join[depthdummy,{PointedListRank3Adm[a,b,c][[idepth]][[2]]}]]] ;depthdummy]/;Rank3AdmQ[{a,b,c,d}]

The function depth[d][a, b, c] returns a list of the irreps of Va⊗Vb indexed by i so Vd is an irrep of V_depth[[i]]⊗Vc.  Because of multiplicity issues this will often be a set whose order is greater than one.

The includesummandrank3list will be return a list to be summed which will be the generalization of the rank 2 includesummand.  k is the basic element of index k in Vd inside VaVbVc and i is the ith occurance of the inclusion in VaVbVc since it can occur more than once.

PreIncludeSummandRank3List[k_][{a_,b_,c_,d_},i_,{x_,y_,z_}]:=  
IncludeSummandRank2List[k][{depth[d][a,b,c][[i]],c,d},{cc,z}]/;HHeadQ[x]&&HHeadQ[y]&&HHeadQ[z]&&Rank3AdmQ[{a,b,c,d}]&&0≤k&&k≤d

PreIncludeDualSummandRank3List[k_][{a_,b_,c_,d_},i_,{x_,y_,z_}]:=    
IncludeDualSummandRank2List[k][{depth[d][a,b,c][[i]],c,d},{CC,z}]/;HHeadQ[x]&&HHeadQ[y]&&HHeadQ[z]&&Rank3AdmQ[{a,b,c,d}]&&0≤k&&k≤d

ccterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,Apply[List,PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]][[2]]][[1]]]/;Head[PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]]===Times
ccterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]][[1]]]/;Head[PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]]===TensorProduct

CCterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,Apply[List,PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]][[2]]][[1]]]/;Head[PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]]===Times
CCterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]][[1]]]/;Head[PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]]===TensorProduct

rrterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,Apply[List,PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]][[2]]][[2]]]/;Head[PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]]===Times
rrterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]][[2]]]/;Head[PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]]===TensorProduct

RRterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,Apply[List,PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]][[2]]][[2]]]/;Head[PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]]===Times
RRterm[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,Apply[List,PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]][[2]]]/;Head[PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]]===TensorProduct

CoefFirstIncl[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]][[1]]/;Head[PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]]===Times
CoefFirstIncl[k_][{a_,b_,c_,d_},i_][j_]:=1/;Head[PreIncludeSummandRank3List[k][{a,b,c,d},i,{pp,qq,rr}][[j]]]===TensorProduct

DualCoefFirstIncl[k_][{a_,b_,c_,d_},i_][j_]:=Apply[List,PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]][[1]]/;Head[PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]]===Times
DualCoefFirstIncl[k_][{a_,b_,c_,d_},i_][j_]:=1/;Head[PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{PP,QQ,RR}][[j]]]===TensorProduct

IncludeSummandRank3List[k_][{a_,b_,c_,d_},i_,{x_,y_,z_}]:=IncludeSummandRank3List[k][{a,b,c,d},i,{x,y,z}]=Module[{includedummy={}},For[includei=1,includei<Length[PreIncludeSummandRank3List[k][{a,b,c,d},i,{x,y,z}]]+1,includei++,For[includej=1,includej<Length[IncludeSummandRank2List[ccterm[k][{a,b,c,d},i][includei][[1]]][{a,b,ccterm[k][{a,b,c,d},i][includei][[2]]},{pp,qq}]]+1,includej++,includedummy=Join[includedummy,{
TensorTimes[Binomial[d,k]*CoefFirstIncl[k][{a,b,c,d},i][includei]*IncludeSummandRank2List[ccterm[k][{a,b,c,d},i][includei][[1]]][{a,b,ccterm[k][{a,b,c,d},i][includei][[2]]},{pp,qq}][[includej]],z[rrterm[k][{a,b,c,d},i][includei][[1]],rrterm[k][{a,b,c,d},i][includei][[2]]]]}]]];
includedummy]

IncludeDualSummandRank3List[k_][{a_,b_,c_,d_},i_,{x_,y_,z_}]:=IncludeDualSummandRank3List[k][{a,b,c,d},i,{x,y,z}]=Module[{includedummy={}},For[includei=1,includei<Length[PreIncludeDualSummandRank3List[k][{a,b,c,d},i,{x,y,z}]]+1,includei++,For[includej=1,includej<Length[IncludeDualSummandRank2List[CCterm[k][{a,b,c,d},i][includei][[1]]][{a,b,CCterm[k][{a,b,c,d},i][includei][[2]]},{PP,QQ}]]+1,includej++,includedummy=Join[includedummy,{
TensorTimes[DualCoefFirstIncl[k][{a,b,c,d},i][includei]*IncludeDualSummandRank2List[CCterm[k][{a,b,c,d},i][includei][[1]]][{a,b,CCterm[k][{a,b,c,d},i][includei][[2]]},{PP,QQ}][[includej]],z[RRterm[k][{a,b,c,d},i][includei][[1]],RRterm[k][{a,b,c,d},i][includei][[2]]]]}]]];
includedummy]

IncludeSummandRank3[k_][{a_,b_,c_,d_},i_,{x_,y_,z_}]:=IncludeSummandRank3[k][{a,b,c,d},i,{x,y,z}]=Apply[Plus,IncludeSummandRank3List[k][{a,b,c,d},i,{x,y,z}]]

IncludeDualSummandRank3[k_][{a_,b_,c_,d_},i_,{x_,y_,z_}]:=IncludeDualSummandRank3[k][{a,b,c,d},i,{x,y,z}]=Apply[Plus,IncludeDualSummandRank3List[k][{a,b,c,d},i,{x,y,z}]]

Computing Central Functions

EvalMatCoeff[TensorProduct[TensorProduct[pp[i1_,m_],qq[j1_,n_],rr[k1_,l_]],TensorProduct[PP[i2_,m_],QQ[j2_,n_],RR[k2_,l_]] 
],{pMatrix_,qMatrix_,rMatrix_}]:=EvalMatCoeff[TensorProduct[TensorProduct[pp[i1,m],qq[j1,n],rr[k1,l]],TensorProduct[PP[i2,m],QQ[j2,n],RR[k2,l]] 
],{pMatrix,qMatrix,rMatrix}]=pMatrix[[i2+1,i1+1]] qMatrix[[j2+1,j1+1]]rMatrix[[k2+1,k1+1]]

pqrMatrix[{a_,b_,c_,d_}]:=pqrMatrix[{a,b,c,d}]={SymmetricPower[T1,a],SymmetricPower[T2,b],SymmetricPower[T3[s],c]}

The sample matrices X, Y are defined over the quadratic extension  of the polynomial ring
C[t1,t2]] generated by  w1, w2 such that   w1 w2 = 1 and w1 + w2 = t12. The routine wPolyNew expresses a polynomial in t1, t2, w1, w2 which is symmetric in w1 and w2 as a polynomial in t1, t2, t12.

The (n^th) Chebychev polynomial Cheb[n,w] expresses (((w1))^n) + (((w2))^n) as a polynomial in w = w1 + w2 where w1 w2 = 1.

Cheb[0,w_]=2;
Cheb[1,w_] = w;
Cheb[n_,w_] := Expand[w* Cheb[n-1,w] - Cheb[n-2,w]]

SymmetryTest[f_,n_]:=Coefficient[f,w,n]===Coefficient[f,w,-n];
SymmetryTest[f_]:=Apply[And,Map[SymmetryTest[f,#]&,Range[0,Exponent[f,w]]]]

wPolyNew[f_]:= Module[{f1=f},Expand[Coefficient[f1,w,0]+Apply[Plus,Map[Cheb[#,t12] Coefficient[f1,w,#]&,Range[Exponent[f,w]]]]]]/;SymmetryTest[f]

Put it all together:

Casimir[{a_,b_,c_,d_},i_,j_]:=Casimir[{a,b,c,d},i,j]=Inner[TensorProduct,IncludeSummandRank3[{a,b,c,d},i],IncludeDualSummandRank3[{a,b,c,d},j]]

TensorContraction[{a_,b_,c_,d_},i_,j_]:=TensorContraction[{a,b,c,d},i,j]=Expand[EvalMatCoeff[Casimir[{a,b,c,d},i,j],pqrMatrix[{a,b,c,d}]]]

ContractionReduce[E_] := Expand[Simplify[E -FullSimplify[(Coefficient[E, s, Exponent[E, s]])/(Coefficient[Expand[Power[Det[T3[s]] - 1, (Exponent[E, s]/2)]], s, Exponent[E, s]])]*Expand[Power[Det[T3[s]] - 1, (Exponent[E, s]/2)]]]] /; IntegerQ[Exponent[E, s]/2] && Exponent[E, s] ≥2
ContractionReduce[E_] :=E/;Exponent[E, s] ≤1
RecursiveContractionReduce[E_]:=Nest[ContractionReduce,E,(Exponent[E,s]/2)]/; IntegerQ[Exponent[E, s]/2] && Exponent[E, s] ≥2
RecursiveContractionReduce[E_]:=E/;Exponent[E, s] ≤1

ContractionReduce[E_] := Expand[Simplify[E -FullSimplify[(Coefficient[E, s, Exponent[E, s]])/(Coefficient[Expand[s*Power[Det[T3[s]] - 1, ((Exponent[E, s]-1)/2)]], s, Exponent[E, s]])]*Expand[s*Power[Det[T3[s]] - 1, ((Exponent[E, s]-1)/2)]]]] /; IntegerQ[(Exponent[E, s]-1)/2] && Exponent[E, s] ≥2
ContractionReduce[E_] :=E/;Exponent[E, s] ≤1
RecursiveContractionReduce[E_]:=Nest[ContractionReduce,E,(Exponent[E,s])]/; IntegerQ[(Exponent[E, s]-1)/2] && Exponent[E, s] ≥2
RecursiveContractionReduce[E_]:=E/;Exponent[E, s] ≤1

The function ContractionReduce uses that fact that Det[T3[s]]-1=0 to ensure that TensorContraction does not have powers of s larger than 4 (they occur as powers of 2).  Beyond degree 4 the FullSimplify command has difficulty.  This also ensures that the largest power of t123 in the final result will be one; possible since there is a monic quadratic relation in t123.  (We are assuming that the coefficients divide each other, but in the function field this is valid and since we are adding zero in the ring we expect this will always work.)  The exponent is always even since the determinant is quadratic and we are symmetrizing over this, so the condition checking that the exponent is even is unnecessary.

PreCentralRank3[a_,b_,c_,d_,i_,j_]:=PreCentralRank3[a,b,c,d,i,j]=Expand[FullSimplify[RecursiveContractionReduce[TensorContraction[{a,b,c,d},i,j]]/.ssub]/.quadsub]

CentralRank3[a_,b_,c_,d_,i_,j_]:=(Expand[wPolyNew[PreCentralRank3[a,b,c,d,i,j]]]/.{t1->(Subscript[t, 1]),t2->(Subscript[t, 2]),t3->(Subscript[t, 3]),t12->(Subscript[t, 12]),t13->(Subscript[t, 13]),t23->(Subscript[t, 23]),t123->(Subscript[t, 123])})/;Rank3AdmQ[{a,b,c,d}]&&i≥1&&i≤Length[depth[d][a,b,c]]&&j≥1&&j≤Length[depth[d][a,b,c]]

CentralRank2[a_,b_,c_]:=CentralRank3[a,b,0,c,1,1]
CentralRank1[a_]:=CentralRank3[a,0,0,a,1,1]

Central[a_,b_,c_,d_,e_,f_]:=Central[a,b,c,d,e,f]=CentralRank3[a,b,c,d,e,f]
Central[a_,b_,c_]:=Central[a,b,c]=CentralRank2[a,b,c]
Central[a_]:=Central[a]=CentralRank1[a]

Computations

The following routine calculates all admissible six-tuples so a+b+c=n.  We then implement this to calculate the a+b+c=n cases.

PreTriples[n_]:=PreTriples[n]=Module[{dummypretrip={}},For[ipretrip=1,ipretrip<Length[IntegerPartitions[n,3]]+1,ipretrip++,If[Length[IntegerPartitions[n,3][[ipretrip]]]==1,dummypretrip=Join[dummypretrip,{{IntegerPartitions[n,3][[ipretrip]][[1]],0,0}}],If[Length[IntegerPartitions[n,3][[ipretrip]]]==2,dummypretrip=Join[dummypretrip,{{IntegerPartitions[n,3][[ipretrip]][[1]],IntegerPartitions[n,3][[ipretrip]][[2]],0}}],dummypretrip=Join[dummypretrip,{IntegerPartitions[n,3][[ipretrip]]}]]]]; dummypretrip]
Triples[n_]:=Triples[n]=Module[{dummytrip={}},For[itrip=1,itrip<Length[PreTriples[n]]+1,itrip++,dummytrip=Join[dummytrip,Permutations[PreTriples[n][[itrip]]]]];dummytrip]
Quadruples[n_]:=Quadruples[n]=Module[{dummyquad={}},For[iquad=1,iquad<Length[Triples[n]]+1,iquad++,For[jquad=1,jquad<Length[ListRank3Adm[Triples[n][[iquad]][[1]],Triples[n][[iquad]][[2]],Triples[n][[iquad]][[3]]]]+1,jquad++,dummyquad=Union[dummyquad,{{Triples[n][[iquad]][[1]],Triples[n][[iquad]][[2]],Triples[n][[iquad]][[3]],ListRank3Adm[Triples[n][[iquad]][[1]],Triples[n][[iquad]][[2]],Triples[n][[iquad]][[3]]][[jquad]]}}]]];dummyquad]
SixTuple[n_]:=SixTuple[n]=Module[{dummysixtup={}},For[isixtup=1,isixtup<Length[Quadruples[n]]+1,isixtup++,For[jsixtup=1,jsixtup<Mult[Quadruples[n][[isixtup]][[1]],Quadruples[n][[isixtup]][[2]],Quadruples[n][[isixtup]][[3]]][Quadruples[n][[isixtup]][[4]]]+1,jsixtup++,For[ksixtup=1,ksixtup<Mult[Quadruples[n][[isixtup]][[1]],Quadruples[n][[isixtup]][[2]],Quadruples[n][[isixtup]][[3]]][Quadruples[n][[isixtup]][[4]]]+1,ksixtup++,dummysixtup=Join[dummysixtup,{{Quadruples[n][[isixtup]][[1]],Quadruples[n][[isixtup]][[2]],Quadruples[n][[isixtup]][[3]],Quadruples[n][[isixtup]][[4]],jsixtup,ksixtup}}]]]];dummysixtup]
SixTuple[0]:={{0,0,0,0,1,1}}
Central[{a_,b_,c_,d_,e_,f_}]:=Central[a,b,c,d,e,f]

CentralStratum[n_]:=For[ics = Length[SixTuple[n]], ics > 0, ics--, Print[χ[SixTuple[n][[ics]][[1]],SixTuple[n][[ics]][[2]],SixTuple[n][[ics]][[3]],SixTuple[n][[ics]][[4]],SixTuple[n][[ics]][[5]],SixTuple[n][[ics]][[6]]] == Central[SixTuple[n][[ics]]]]]

We list examples by a+b+c=n

a+b+c=0

CentralStratum[0]

a+b+c=1

CentralStratum[1]

a+b+c=2

CentralStratum[2]

a+b+c=3

CentralStratum[3]

a+b+c=4

CentralStratum[4]