Issue
The main issue here is you can't add symbolic tensors with explicit tensors. In your example, d is a symbolic tensor, but IdentityMatrix[3] is an explicit tensor. Note what happens when you add them:
$Assumptions = Element[d, Matrices[{3, 3}, Reals, Symmetric[{1,2}]]];
d + IdentityMatrix[3]
{{1 + d, d, d}, {d, 1 + d, d}, {d, d, 1 + d}}
The Listable attribute of Plus causes d to be added to all of the matrix elements of IdentityMatrix[3].
Symbolic identity matrix
Instead of using the explicit tensor IdentityMatrix[3], you can use either IdentityMatrix[n] or Inactive[IdentityMatrix][3]. For example:
d + IdentityMatrix[n]
d + Inactive[IdentityMatrix][3]
d + IdentityMatrix[n]
d + Inactive[IdentityMatrix][3]
In both cases no explicit matrix is generated. Now, for your example, it is better to use Inactive[IdentityMatrix][3] so that dimensions match up. If you had used Element[d, Matrices[{n, n}]] instead, you could use IdentityMatrix[n]. So:
dd = d - Tr[d] Inactive[IdentityMatrix][3]/3;
TensorReduce @ TensorContract[TensorProduct[dd, dd], {{1, 3}, {2, 4}}]
TensorContract[TensorProduct[d, d], {{1, 3}, {2, 4}}] -
2/3 TensorContract[
TensorProduct[d, Inactive[IdentityMatrix][3]], {{1, 3}, {2, 4}}] Tr[
d] + 1/9 TensorContract[
TensorProduct[Inactive[IdentityMatrix][
3], Inactive[IdentityMatrix][3]], {{1, 3}, {2,
4}}] Tr[d]^2
Unfortunately, at the present time, TensorReduce doesn't know what to do with these "symbolic" identity tensors. Therefore, we will need to help things along.
SymbolicTensors`IdentityTensor
One possibility is to make use of the internal symbol SymbolicTensors`IdentityTensor. When a single such symbol is involved in a tensor contraction, then TensorContract will properly simplify the tensor product. Here is an example:
tensor = TensorContract[TensorProduct[d, Inactive[IdentityMatrix][3]], {{1, 3}, {2, 4}}];
tensor /. Inactive[IdentityMatrix] -> SymbolicTensors`IdentityTensor
TensorContract[d, {{1, 2}}]
On the other hand, if the TensorProduct contains a "symbolic" identity tensor (i.e., Inactive[IdentityMatrix][3] or IdentityMatrix[n]) that is not contracted, then TensorContract produces an error:
tensor = TensorContract[TensorProduct[d, Inactive[IdentityMatrix][3], d], {{1, 5}}];
tensor /. Inactive[IdentityMatrix] -> SymbolicTensors`IdentityTensor
Part::pkspec1: The expression {SymbolicTensors`SymbolicTensorsDump`i} cannot be used as a part specification.
...
where I've suppressed most of the output. Similarly, if the TensorProduct contains multiple "symbolic" identity tensors, even if they are all contracted, TensorContract still produces an error:
tensor = TensorContract[TensorProduct[Inactive[IdentityMatrix][3], Inactive[IdentityMatrix][3]], {{1, 4}, {2, 3}}];
tensor /. Inactive[IdentityMatrix] -> SymbolicTensors`IdentityTensor
Part::partw: Part {2,3} of {3[1,1],2[1,2]} does not exist.
...
So, what is needed is a method to find "symbolic" identity tensors that are being contracted, and one at a time, replace those tensors with SymbolicTensors`IdentityTensor.
IdentityReduce
Update - The paclet now lives on GitHub, and can be installed using:
PacletInstall[
"TensorSimplify",
"Site" -> "http://raw.githubusercontent.com/carlwoll/TensorSimplify/master"
]
Here is a function that does this (soon to be a part of a TensorSimplify package I will be putting on GitHub):
IdentityReduce[expr_] := expr /. TensorContract -> ir
ir[t_TensorProduct, i_] := Module[{indices, imIndices, ids, dummy},
indices = tensorIndices[t];
If[indices === $Failed, Return[TensorContract[t, i]]];
imIndices = Position[
Replace[t, Inactive[IdentityMatrix][n_] -> dummy[n], {1}],
IdentityMatrix | dummy,
{2}
];
ids = Pick[
imIndices,
IntersectingQ[Flatten[i],indices[#]]& /@ imIndices[[All,1]]
];
If[ids==={},
TensorContract[t,i],
TensorContract[ReplacePart[t, ids[[1]]->SymbolicTensors`IdentityTensor],i] /. TensorContract->ir
]
]
ir[(Inactive[IdentityMatrix] | IdentityMatrix)[n_], {{1,2}}]=n;
ir[o_, i_]:=TensorContract[o,i]
tensorIndices[Verbatim[TensorProduct][t__]] := With[{r=Accumulate @* Map[TensorRank] @ {1,t}},
If[MatchQ[r, {__Integer}],
Association @ Thread @ Rule[
Range@Length[{t}],
Range[1+Most[r], Rest[r]]
],
$Failed
]
]
Let's try the OP example again:
e = TensorReduce @ TensorContract[TensorProduct[dd, dd], {{1, 3}, {2, 4}}];
IdentityReduce @ e
TensorContract[TensorProduct[d, d], {{1, 3}, {2, 4}}] -
2/3 TensorContract[d, {{1, 2}}] Tr[d] + Tr[d]^2/3
Notice how the symbolic identity tensors have been resolved. One can perform one final step, converting all TensorContract objects to Dot/Tr (if possible) using my FromTensor function:
FromTensor @ IdentityReduce @ e
-(1/3) Tr[d]^2 + Tr[Transpose[d].d]
It is possible to use TensorReduce on this expression, but it produces a result using Tr[MatrixPower[d, 2]] instead of Tr[d.d], so I won't do so.
