Probably not the most efficient way to do so, but for descently sized `n`, this should work well: ArrayReshape[ T[ρ, #] & /@ Tuples[Range[3], n], ConstantArray[3, n] ]