You can use the new in M10 Tensor* function to do this, although it is a bit long-winded. First, lets define your identity matrix:
id = IdentityMatrix[d];
Next, we can represent your individual terms using TensorProduct
and TensorTranspose
. For example,
$\delta _{a,b} \delta _{c,d} \delta _{e,f}$:
TensorProduct[id, id, id];
$\delta _{a,d} \delta _{e,b} \delta _{c,f}$:
TensorTranspose[TensorProduct[id, id, id], {1, 4, 3, 6, 5, 2}];
and similarly for the others. Finally, we want to contract the expression:
res = TensorContract[
TensorProduct[
(
TensorProduct[id, id, id] +
3 TensorTranspose[TensorProduct[id, id, id], {1, 4, 3, 6, 5, 2}]
),
(
TensorProduct[id, id, id] +
3 TensorTranspose[TensorProduct[id, id, id], {1, 4, 3, 6, 5, 2}]
)
],
{{3, 9}, {4, 10}, {5, 11}, {6, 12}}
];
Unfortunately, TensorReduce
is not able to simplify this further. However, we can make use of my TensorSimplify
package to help out. Install the paclet with::
PacletInstall[
"TensorSimplify",
"Site"->"http://raw.githubusercontent.com/carlwoll/TensorSimplify/master"
]
and then load the package with:
<<TensorSimplify`
Then:
TensorSimplify[res]
6 IdentityMatrix[d]\[TensorProduct]IdentityMatrix[d] +
d^2 IdentityMatrix[d]\[TensorProduct]IdentityMatrix[d] +
9 d TensorTranspose[
IdentityMatrix[d]\[TensorProduct]IdentityMatrix[d], {1, 3, 2, 4}]
In terms of kronecker deltas, the above is equivalent to:
$$d^2 \delta _{a,b} \delta _{i,j}+9 d \delta _{a,i} \delta _{b,j}+6 \delta_{a,b} \delta _{i,j}$$