I would like to calculate the following transformation of a fourth rank tensor,
$$ C_{ijkl}=\Sigma_{m=1}^{3}\Sigma_{n=1}^{3}\Sigma_{p=1}^{3}\Sigma_{q=1}^{3}a_{im}a_{jn}a_{kp}a_{lq}C_{mnpq} $$
where $a_{xy}$ are direction cosines, or as some other literature calls them, components of the rotation matrix (Are those same?)
But this fourth rank tensor can be represented as a 6 by 6 matrix, by using
$$ 11\rightarrow 1\qquad 22\rightarrow 2\qquad 33\rightarrow 3\qquad 23,32\rightarrow 4\qquad 31,13\rightarrow 5\qquad 12,21\rightarrow 6 $$
$$\left[ \begin{array}{cccccc} c_{11}&c_{12}&c_{13}&c_{14}&c_{15}&c_{16}\\ c_{21}&c_{22}&c_{23}&c_{24}&c_{25}&c_{26}\\ c_{31}&c_{32}&c_{33}&c_{34}&c_{35}&c_{36}\\ c_{41}&c_{42}&c_{43}&c_{44}&c_{45}&c_{46}\\ c_{51}&c_{52}&c_{53}&c_{54}&c_{55}&c_{56}\\ c_{61}&c_{62}&c_{63}&c_{64}&c_{65}&c_{66} \end{array}\right] $$
So my question is, is there any simpler form, maybe a matrix form of the equation $$ C_{ijkl}=\Sigma_{m=1}^{3}\Sigma_{n=1}^{3}\Sigma_{p=1}^{3}\Sigma_{q=1}^{3}a_{im}a_{jn}a_{kp}a_{lq}C_{mnpq} $$
Otherwise, can the above equation be represented in the form of a matrix?
This equation is too hard to conceive just by looking at it.
Or, maybe, is there any expansion of this equation?
I have no idea how you are getting a rank-4 tensor from a rotation problem--perhaps you can explain more about that?
Generally the term direction cosine refers to the fact that the rotation group on $\mathbb{R}^{n>1}$, e.g. $SO_n$, can always be written as an $n\times n$ matrix representation where $$A_{ij} = \mathbf{b}_i\cdot\mathbf{r}_j,\ \ \ 1\leq i,j \leq n ,\ \ \ A\in SO_n,$$ for two orthonomal bases $\mathscr{B} = \{\mathbf{b}_i\}_i^n$ and $\mathscr{R} = \{\mathbf{r}_i\}_i^n$. Since the bases are orthonormal, for the $n = 2,3$ cases we can meaningfully visualize the fact that $\mathbf{b}_i\cdot\mathbf{r}_j = \cos\theta_{ij},$ where $\theta_{ij}$ is the angle in between basis vectors $\mathbf{b}_i$ and $\mathbf{r}_j$. Thus the generic matrix elements of these rotation groups have a geometric interpretation as the cosine of the angle in between the basis vectors--the so-called "direction cosine".
While I suppose you can represent the 4-index tensor with a matrix--e.g. tabulate the various entries of the tensor this way--this "matrix" has no real advantages from a linear algebra perspective. Unlike for the 2-index case, summing along indices will not be equivalent to matrix/vector multiplication, and you certainly shouldn't expect anything about the spectrum or invariant spaces of the tensor as an operator to hold for this matrix representaion!
Your equation looks like it's trying to rotate each index of your tensor individually. For instance, let's fix all but the first index to the first entry, say 1. Then your map looks like $$C'_{111i} = \sum_{j=1}^{3}a_{ij}C_{111j}$$ which will effectively be like changing the coordinates of an $\mathbb{R^3}$ vector to new coordinates in a rotated basis. By doing this to each index you have the classic way of rotating the entire tensor--e.g. see here or here. Depending on specifically what you want to do you might be able to find a more economical rotation operator, say $\Phi_{ijkl}^{\ \ \ \ \ \ \ \ mpqr}$, for which $$C'_{ijkl} = \sum_{m,p,q,r}\Phi_{ijkl}^{\ \ \ \ \ \ \ \ mpqr}C_{mpqr},$$ but there is no general form for such $\Phi$ I know, save of course the trivial $$\Phi_{ijkl}^{\ \ \ \ \ \ \ \ mpqr} = a_i^ma_j^pa_k^qa_l^r,$$ but even supposing you were, I see no way this hypothetical 8-index operator would somehow possess a matrix representation with which meaningful calculations could be made.