How to use Mathematica to prove that isotropic materials have only two independent parameters
Summary
There are 3 independent degrees of freedom for the tensor $A_{abcd}$ if it is invariant under Rayleigh transformation, and it also satisfied $A_{abcd}=A_{cdab}$. The independent components can be taken as $A_{2332},A_{3223},$ and $A_{3322}$. If we further require that $A_{abcd}=A_{bacd}$ (and/or $A_{abcd}=A_{abdc}$), the number of independent components reduce to 2, which can be taken as $A_{3223}$ and $A_{3322}$. By direct analysis, we can actually see that
$$A_{abcd}=A_{3322}\delta_{ab}\delta_{cd}+A_{3232}\left(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}\right)$$
Note: Below, I'll detail a Mathematica implementation to get this result via using invariance under $SO(3)$ through explicit transformation. However, in practice, we do not need Mathematica to get this result. As we are trying to construct a rank-4 invariant object, $\delta_{ab}$ is the only $SO(n)$ invariant tensor that we can use ($\epsilon_{abc}$ is no use here), hence one can immediately make the ansatz $$A_{abcd}=a\delta_{ab}\delta_{cd}+b\delta_{ac}\delta_{bd}+c\delta_{ad}\delta_{bc}$$ which is what we got through Mathematica anyway. We can get the result above from this by imposing $A_{abcd}=A_{bacd}$, hence the result that there are 2 independent components does not really need Mathematica computation.
Details
One straightforward, albeit not elegant, method is to simply utilize the invariance under $SO(3)$ explicitly.
Let us illustrate this with a simpler example by trying to find a $3\times3$ matrix invariant under rotations: $A_{ab}=M_{ac}A_{cd}(M^T)_{db}$, or analogous to the way you wrote down the Rayleigh product, $A_{ab}=M_{ac}M_{bd}A_{cd}$. To find this matrix, we first define RayleighProd2[i_, input_]
which does this transformation:
RayleighProd2[i_, input_] :=
TensorContract[
TensorProduct[RotationMatrix[Sequence @@ i],
TensorContract[
TensorProduct[RotationMatrix[Sequence @@ i], input],
{2, 4}]],
{2, 4}];
This basically takes an input matrix $A_{ab}$ and takes it into $M_{ac}M_{bd}A_{cd}$ where the rotation matrix $M$ is specified by the first argument of RayleighProd2
, which will be a list of the form {a,v}
where v
is the the vector pointing rotation axis and a
is the angle of rotation.
If we require that this transformation is an identity transformation, that means that the matrix $M_{ac}M_{bd}A_{cd}$ should be independent of a
and v
. For practical purposes, we consider this transformation for all three orthogonal axes for arbitrary angle a
, and then take the derivative wrt a
: the resultant tensor should be zero for any a
, which we will fix to 0 for convenience. We define this operation as
ClearAll[casesToBeChecked];
casesToBeChecked[input_] := Block[{\[Theta]},
Flatten[D[
Flatten[Table[
RayleighProd2[{\[Theta], v}, input],
{v, {{0, 0, 1}, {0, 1, 0}, {1, 0, 0}}}]],
\[Theta]]] /. \[Theta] -> 0];
We can now find the form of the matrix $A_{ab}$ invariant under these rotations:
Array[a, {3, 3}] /. Solve[casesToBeChecked[Array[a, {3, 3}]] == 0]
(* {{{a[2, 2], 0, 0}, {0, a[2, 2], 0}, {0, 0, a[2, 2]}}} *)
We immediately see that $A_{ab}=c\delta_{ab}$. QED.
We can extend this approach to OP's case as follows:
RayleighProd4[i_, input_] :=
TensorContract[
TensorProduct[RotationMatrix[Sequence @@ i],
TensorContract[
TensorProduct[RotationMatrix[Sequence @@ i],
TensorContract[
TensorProduct[RotationMatrix[Sequence @@ i],
TensorContract[
TensorProduct[RotationMatrix[Sequence @@ i], input],
{2, 6}]],
{2, 6}]],
{2, 6}]],
{2, 6}];
and
casesToBeChecked[input_] := Block[{\[Theta]},
Flatten[D[
Flatten[Table[
RayleighProd4[{\[Theta], v}, input],
{v, {{0, 0, 1}, {0, 1, 0}, {1, 0, 0}}}]],
\[Theta]]] /. \[Theta] -> 0];
For which
Array[a, {3, 3, 3, 3}] /. Solve[casesToBeChecked[Array[a, {3, 3, 3, 3}]] == 0]
gives the form of the tensor:
$$\left( \begin{array}{ccc} \left( \begin{array}{ccc} \left( \begin{array}{c} a(2,3,3,2)+a(3,2,3,2)+a(3,3,2,2) \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ a(3,3,2,2) \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ a(3,3,2,2) \\ \end{array} \right) \\ \left( \begin{array}{c} 0 \\ a(3,2,3,2) \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} a(2,3,3,2) \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) \\ \left( \begin{array}{c} 0 \\ 0 \\ a(3,2,3,2) \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} a(2,3,3,2) \\ 0 \\ 0 \\ \end{array} \right) \\ \end{array} \right) & \left( \begin{array}{ccc} \left( \begin{array}{c} 0 \\ a(2,3,3,2) \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} a(3,2,3,2) \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) \\ \left( \begin{array}{c} a(3,3,2,2) \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ a(2,3,3,2)+a(3,2,3,2)+a(3,3,2,2) \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ a(3,3,2,2) \\ \end{array} \right) \\ \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ a(3,2,3,2) \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ a(2,3,3,2) \\ 0 \\ \end{array} \right) \\ \end{array} \right) & \left( \begin{array}{ccc} \left( \begin{array}{c} 0 \\ 0 \\ a(2,3,3,2) \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} a(3,2,3,2) \\ 0 \\ 0 \\ \end{array} \right) \\ \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ a(2,3,3,2) \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ a(3,2,3,2) \\ 0 \\ \end{array} \right) \\ \left( \begin{array}{c} a(3,3,2,2) \\ 0 \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ a(3,3,2,2) \\ 0 \\ \end{array} \right) & \left( \begin{array}{c} 0 \\ 0 \\ a(2,3,3,2)+a(3,2,3,2)+a(3,3,2,2) \\ \end{array} \right) \\ \end{array} \right) \\ \end{array} \right)$$
We can see that it depends on only 3 independent variables:
DeleteCases[
Flatten[Array[a, {3, 3, 3, 3}] /. Solve[casesToBeCheckedNew[Array[a, {3, 3, 3, 3}]] == 0]],
0] // Union
(* {a[2, 3, 3, 2], a[3, 2, 3, 2], a[3, 3, 2, 2], a[2, 3, 3, 2] + a[3, 2, 3, 2] + a[3, 3, 2, 2]} *)
Now, we can impose further constraints regarding the symmetries of $A_{abcd}$. For example, we see that
A = Array[a, {3, 3, 3, 3}] /. Solve[casesToBeChecked[Array[a, {3, 3, 3, 3}]] == 0][[1]];
Union[DeleteCases[Flatten[Transpose[A, {2, 1, 3, 4}] - A], 0]]
(* {a[2, 3, 3, 2] - a[3, 2, 3, 2], -a[2, 3, 3, 2] + a[3, 2, 3, 2]} *)
which means the independent number of variables reduces to 2 if we require $A_{abcd}=A_{bacd}$. The same constraint appears for $A_{abcd}=A_{abdc}$, so the independent parameter number is still 2:
Union[DeleteCases[Flatten[Transpose[A, {1, 2, 4, 3}] - A], 0]]
{a[2, 3, 3, 2] - a[3, 2, 3, 2], -a[2, 3, 3, 2] + a[3, 2, 3, 2]}
Finally, the symmetry $A_{abcd}=A_{cdab}$ seems to be already satisfied:
Union[DeleteCases[Flatten[Transpose[A, {3, 4, 1, 2}] - A], 0]]
(* {} *)