How to generate isotropic tensors?

Consider a real tensor $bf{T}$ upon $mathopotimeslimits_{n}(mathbb{R}^3)$, with a definite order $n$, it could be isotropic, i.e., be invariant under the action of all elements in $SO(3)$(not $O(3)$, so pseudo-tensor is allowed), only if which matches a certain pattern.

e.g., with $n=1$, thus $bf{T}$ is a vector in $mathbb{R}^3$. The only possible isotropic vector must be zero vector, $bf0$.

$n=2$, is $lambda bf{I}$, where $bf{I}$ is identity tensor, or $lambdadelta_{ij}$, in component. $delta_{ij}$ is Kronecker delta.

$n=3$, is $lambda epsilon_{ijk}$.

$n=4$, is $alpha delta_{ij}delta_{kl}+betadelta_{ik}delta_{jl}+gammadelta_{il}delta_{jk}$.

I use this piece of code to get the pattern:

tensor[n_] := tensor[n] = Array[Symbol[
    StringJoin["t", Sequence @@ ToString /@ {##}]
]&, ConstantArray[3, n]];

this makes an abstract tensor in order $n$.

For order $n$, the result $T’_{i_1′ i_2’dots i_n’}$ after above transformation is
$$T’_{i_1′ i_2’dots i_n’}=R_{i_1’i_1}R_{i_2’i_2}dots R_{i_n’i_n}T_{i_1i_2dots i_n}$$
where $R$ is the rotation matrix in $SO(3)$.

Then use:

matrixZ[q_] := RotationMatrix[q, {0, 0, 1}]
matrixY[q_] := RotationMatrix[q, {0, 1, 0}]

to generate rotation matrices in two orthogonal axes. Tensor is isotropic iff its invariant under those matrices(for arbitrary q)(cause any rotation can be
decompose into three rotations in these two axes, first around $z$, then $y$, and finally $z$ again).

res[n_, mat_, ts_] := TensorContract[
    Outer[Times, ts, Sequence@@(ConstantArray[mat, n])],
    {#, 2#+n}& /@ Range[n]
]

give the result $mathbf{T}’$, corresponding to rotation matrix mat, tensor ts and order n.

Now I can solve it(e.g. order $3$):

tensor3 = tensor[3];
matZ = matrixZ[Pi/2];
matY = matrixY[Pi/2];
resZ = res[3, matZ, tensor3];
resY = res[3, matY, tensor3]
tensor3 /. First @ Solve[
    tensor3 == resZ == resY,
    Flatten @ tensor3
]
(*
{{{0, 0, 0}, {0, 0, t123}, {0, -t123, 0}}, 
 {{0, 0, -t123}, {0, 0, 0}, {t123, 0, 0}}, 
 {{0, t123, 0}, {-t123, 0, 0}, {0, 0, 0}}}
*)

Questions
This method is sometimes very slow, seems the performance depends on the value of q of the matrix.

For example, in order $4$:

tensor4=tensor[4];
matZ1=matrixZ[Pi/4];
matY1=matrixY[Pi/4];
matZ2=matrixZ[Pi/2];
matY2=matrixY[Pi/2];
resZ1=res[4,matZ1,tensor4];
resY1=res[4,matY1,tensor4];
resZ2=res[4,matZ2,tensor4];
resY2=res[4,matY2,tensor4];

If I chose matZ1, matY1, the symbolic way will take too long time to get solution(over five minute), but for matZ1, matY2 only 0.0711628s:

AbsoluteTiming@Solve[{
    tensor4 == resY2 == resZ1}, Flatten@tensor4, Method->Reduce]

{0.0779194, {{t1111 -> t1122 + t1212 + t1221, t1112 -> 0, 
    t1113 -> 0, t1121 -> 0, t1123 -> 0, t1131 -> 0, t1132 -> 0, 
    t1133 -> t1122, t1211 -> 0, t1213 -> 0, t1222 -> 0, 
    t1223 -> 0, t1231 -> 0, t1232 -> 0, t1233 -> 0, t1311 -> 0, 
    t1312 -> 0, t1313 -> t1212, t1321 -> 0, t1322 -> 0, 
    t1323 -> 0, t1331 -> t1221, t1332 -> 0, t1333 -> 0, 
    t2111 -> 0, t2112 -> t1221, t2113 -> 0, t2121 -> t1212, 
    t2122 -> 0, t2123 -> 0, t2131 -> 0, t2132 -> 0, t2133 -> 0, 
    t2211 -> t1122, t2212 -> 0, t2213 -> 0, t2221 -> 0, 
    t2222 -> t1122 + t1212 + t1221, t2223 -> 0, t2231 -> 0, 
    t2232 -> 0, t2233 -> t1122, t2311 -> 0, t2312 -> 0, 
    t2313 -> 0, t2321 -> 0, t2322 -> 0, t2323 -> t1212, 
    t2331 -> 0, t2332 -> t1221, t2333 -> 0, t3111 -> 0, 
    t3112 -> 0, t3113 -> t1221, t3121 -> 0, t3122 -> 0, 
    t3123 -> 0, t3131 -> t1212, t3132 -> 0, t3133 -> 0, 
    t3211 -> 0, t3212 -> 0, t3213 -> 0, t3221 -> 0, t3222 -> 0, 
    t3223 -> t1221, t3231 -> 0, t3232 -> t1212, t3233 -> 0, 
    t3311 -> t1122, t3312 -> 0, t3313 -> 0, t3321 -> 0, 
    t3322 -> t1122, t3323 -> 0, t3331 -> 0, t3332 -> 0, 
    t3333 -> t1122 + t1212 + t1221}}}

Can we find isotropic tensor patterns faster(like kill duplicate equations)? Is there any way to avoid choosing parameters(e.g., matX, matY etc.) manually?

What’s the universal formula of isotropic tensors?

generalIsotropicTensor[rank_, syms_ : {}, symbol_ : A] := 
 Block[
{a = symmetrizedTensor[rank, syms, symbol], 
[Epsilon] = LeviCivitaTensor[3],
isotropyRequirements},
  isotropyRequirements = 
   Normal[Plus @@ (MapIndexed[
       TensorTranspose[
         TensorContract[TensorProduct[[Epsilon], a], #1], 
         Permute[Range[rank + 1], 
          Cycles[{Range[2, First[#2] + 1]}]]] &, 
       Table[{3, 3 + i}, {i, rank}]])];

  a /. Quiet[
    First[Solve[Thread[Flatten[isotropyRequirements] == 0], 
      Variables[a]]]]
  ]

symmetrizedTensor[rank_, syms_ : {}, symbol_ : A] := 
 Normal[SymmetrizedArray[
   pos_ :> Subscript[symbol, StringJoin[ToString /@ pos]], 
   ConstantArray[3, rank], syms]]

MatrixForm@*generalIsotropicTensor /@ Range[4]

generalIsotropicTensor[4, {{Cycles[{{1, 2}}], 1}}]