I have what I think is a relatively simple problem, at least for the case at hand it should be (my group is only has 12 elements). I provide some specific details below for my problem, but skip to below if you don't care about geometric point groups, the question does not rely explicitly on this particular example.
Background
The group $D_{3d}$ is the symmetry group of a triangular anti-prism, which looks like this:
I have highlighted the $C_3$ axis (120 degree rotations) in black and the three $C_3$ axes (120 degree rotation) in blue. The red axes aren't important for this discussion. The set of all of the symmetries are the identity $\mathfrak{e}$, inversion through the midpoint, $\mathfrak{i}$, two $C_3$ rotations about the black axis, three $C_2$ rotations (one for each blue axis), two roto-inversions $S_6$ ($C_3$ rotation followed by inversion), and three mirror planes $\sigma_d$ ($C_2$ rotation followed by inversion). There are twelve total elements of the group.
I can form a 6-dimensional representation of this group which is the set of permutation matrices for the six outer corners, which looks like
I can then write them all as permutations, $$ \begin{array}{|c|c|c|} \hline \text{name} & \text{order} & \text{permutation} \\\hline \mathfrak{e} & 1 & \text{no permutation} \\\hline \mathfrak{i} & 1 & (14)(25)(36) \\\hline C_3 & 2 & { (135)(246)\\ (531)(642) } \\\hline C_2 & 3 & { (14)(23)(56)\\ (25)(16)(34)\\ (36)(12)(45) } \\\hline \sigma_d & 3 & { (26)(35)\\ (13)(46)\\ (15)(24) } \\\hline S_6 & 2 & { (123456)\\ (654321) } \\\hline \end{array} $$
This representation is reducible, and we can calculate how it decomposes using the character table of $D_{3d}$
$$ \begin{array}{|c|c|c|c|c|c|c|} \hline & \mathfrak{e} & 2C_3 & 3C_2 & \mathfrak{i} & 2S_6 & 3\sigma_d \\ \hline A_{1g} & 1 & 1 & 1 & 1 & 1 & 1 \\ \hline A_{1u} & 1 & 1 & 1 & -1 & -1 & -1 \\ \hline A_{2g} & 1 & 1 & -1 & 1 & 1 & -1 \\ \hline A_{2u} & 1 & 1 & -1 & -1 & -1 & 1 \\ \hline E_g & 2 & -1 & 0 & 2 & -1 & 0 \\ \hline E_u & 2 & -1 & 0 & -2 & 1 & 0 \\ \hline \end{array} $$ so that this representation $(R)$ decomposes as a sum of irreps $r$ (whose labels appear in the left column), $$R = \oplus_{r}\, n_r \,r$$ The numbers of time each irrep appears are found from the formula $$n_r = \frac{1}{12} \sum_{g\in D_{3d}} \chi_r^*(g) \chi_R(g)$$ where $g$ are the group elements, $\chi_r(g)$ is its character from the character table, and $\chi_R(g)$ is the trace (character) of the 6-dimensional representation matrices. I find from this that the 6-dimensional permutation representation decomposes as $$R = A_{1g}\oplus A_{2u} \oplus E_g \oplus E_u$$
Question
Now my question is: how do I obtain the linear combinations of corners which transform in the irreps of dimension greater than 1?
My attempt
More specifically, I start with the free vector space whose basis are the corners labelled $i$, so a vector takes the form $$\vert \psi \rangle = \sum_{i=1}^6 c_i \vert i \rangle \quad \, c_i \in \mathbb{R}$$ I want to find a basis for each of the invariant subspaces. This is easy to do for the 1d irreps, because I know that they act by simply multiplying by the respective characters: construct a generic linear combination of the basis vectors as above and ensure that under each permutation it changes by the appropriate sign.
But I cannot figure out a general way to do this for the 2d irreps, since I don't explicitly have the $2\times 2$ matrices for those irreps. My attempt to do this was by explicitly constructing the transformation: for a given permutation $p\in R$ which send $i\to p(i)$, I need to enforce that $$ \begin{pmatrix} \sum_i a_{p(i)} \vert i \rangle \\[10pt] \sum_i b_{p(i)} \vert i \rangle \end{pmatrix} = \begin{pmatrix} A & B \\[10pt] C & \chi_{r}(p)-A \end{pmatrix} \begin{pmatrix} \sum_i a_{i} \vert i \rangle \\[10pt] \sum_i b_{i} \vert i \rangle \end{pmatrix} \quad \forall \, p \in R $$ where $r$ is either $E_g$ or $E_u$, and the vectors $\sum_i a_i \vert i \rangle$ and $\sum_i b_i \vert i \rangle$ ideally form an orthonormal basis for the invariant subspace (of course the $a_i$ and $b_i$ are only determined up to a rotation within the subspace). But I don't know what the 2d representation matrix, i.e. what the constants $A,B,C$ are for the different symmetry elements. It seems this is not a good way to go about finding the basis vectors, is there a straightforward way to construct basis vectors for different subspaces?