The brute-force method to find the symmetric polynomials of a point group is to start with the generic polynomial and perform a symmetrization over group elements. This can be done for small finite groups and does not require much knowledge about the representation theory.
Let $G$ be the point group (can be $O_h$, $T_d$, $D_{6h}$ etc) in consideration. For any element $g\in G$, first write down its three-dimensional representation $R(g)$, as a $3\times 3$ matrix acting on the vector $\Gamma=(\Gamma_x,\Gamma_y,\Gamma_z)$, s.t.
$$(g\Gamma)_i=\sum_{j}R(g)_{ij}\Gamma_j.$$
Then take a generic polynomial expansion of the free energy:
$$F = A^{(0)}_{ijk\cdots}\Gamma_{i}\Gamma_{j}\Gamma_{k}\cdots,$$
symmetrize the coefficient tensor $A^{(0)}_{ijk\cdots}$ by
$$A^{(0)}_{ijk\cdots}\to A^{(1)}_{ijk\cdots}=\frac{1}{|G|}\sum_{g\in G}A^{(0)}_{i'j'k'\cdots}R(g)_{ii'}R(g)_{jj'}R(g)_{kk'}\cdots ,$$
where $|G|$ is the order of the group $G$ (the number of elements in $G$). Then the symmetrized coefficients $A^{(1)}_{ijk\cdots}$ will produce all (and only) symmetric polynomials to the desired order.
Example: Let us consider a simple example of the tetrahedral group $T_d$, which has 24 group elements. They are generated by two generators:
- the two-fold reflection (denoted as $a$) about the mirror plane normal to the $(1,1,0)$ vector,
- the three-fold roation (denoted as $b$) around the $(1,1,1)$ axis by $2\pi/3$.
With these, the tetrahedral group has the presentation $T_d=\langle a, b|a^2=b^3=(ab)^4=1\rangle$. To obtain the matrix representations $R(g)$ for all group elements $g\in T_d$, we first find the representations for the generators $a$ and $b$, which can be done in Mathematica quite easily:
a = ReflectionMatrix[{1, 1, 0}];
b = RotationMatrix[2 Pi/3, {1, 1, 1}];
Then the whole $T_d$ group can be generated from $a$ and $b$ by consecutive multiplications until closure, which can be formulated as a problem of finding the fixed-point under right multiplication of the generators.
G = FixedPoint[DeleteDuplicates@Join[#, #.a, #.b] &, {IdentityMatrix[3]}];
$G$ now contains the 24 elements of $T_d$ group. Now we can use it to symmetrize the coefficient tensors $A^{(0)}_{ijk\cdots}$. Let us start from rank-2 tensors (corresponding to quadratic polynomials). How should we initialize the tensor $A^{(0)}_{ij}$? One nice trick is to initialize $A^{(0)}_{ij}$ by random numbers!
A0 = RandomReal[{}, {3, 3}];
d = ArrayDepth[A0];
The tensor rank (depth) is kept as $d$, which will be useful later. Given the group $G$ and the tensor rank $d$, we can construct the symmetrization operator (also known as the projection operator) $P_G$ for that rank:
PG = Sum[TensorProduct @@ ConstantArray[SparseArray@g, d], {g, G}];
The final step is to apply the symmetrization operator $P_G$ to the random tensor $A^{(0)}$ by tensor contractions:
A1 = TensorContract[TensorProduct[PG, A0], {2 # - 1, 2 d + #} & /@ Range@d];
To read out the symmetric polynomials from the symmetrized tensor $A^{(1)}$, we create an order parameter $r$ to contract with the coefficient tensor.
r = {x, y, z};
Dot @@ Prepend[A1]@ConstantArray[r, d] // Expand
The result will be something like
11.6381 x^2 + 11.6381 y^2 + 11.6381 z^2
After factoring out the common coefficient (which is random), we know that $x^2+y^2+z^2$ is the only $T_d$ symmetric polynomial on the quadratic level. The above procedure can be extended to higher order polynomials by increasing the rank of the initial random tensor $A^{(0)}$. For example, to obtain symmetric quartic polynomials, we need to set
A0 = RandomReal[{}, {3, 3, 3, 3}];
and rerun the program. The result will look like
18.2274 x^4 + 66.6814 x^2 y^2 + 18.2274 y^4 + 66.6814 x^2 z^2
+ 66.6814 y^2 z^2 + 18.2274 z^4
Different symmetric polynomials can be identified by different coefficients. So we know there are two $T_d$ symmetric polynomials on the quartic level: $x^4+y^4+z^4$ and $x^2 y^2+x^2 z^2+y^2 z^2$. In this way, one can find the $T_d$ symmetric polynomials for the lowest few orders:
$$\begin{array}{cc}\text{order} & \text{symmetric polynomials}\\ 2 & x^2+y^2+z^2\\ 3 & xyz\\ 4 & x^4+y^4+z^4\\ & x^2 y^2+x^2 z^2+y^2 z^2\\ 5 & xyz(x^2+y^2+z^2)\\6 & x^6+y^6+z^6\\ & x^4(y^2+z^2)+y^4(x^2+z^2)+z^4(x^2+y^2)\\ & x^2 y^2 z^2\\ \cdots & \cdots \end{array}$$
However, it worth mention that the computation complexity grows exponentially with the polynomial order. So very quickly we will get stuck and not able to proceed to higher orders.
A more clever and more systematic construction will require some knowledge of group representation theory. The idea is to start with the polynomial representation of each irreducible representation of the group and build up higher order polynomials by tensoring the irreducible representations. The symmetric polynomials will be found as the trivial representation. Useful tables for these can be found at this website: http://www.webqc.org/symmetrypointgroup-td.html.
