Let $\vec{a}$ and $\vec{b}$ be two arbitrary vectors in 3-space $\Bbb R^3$, and $R$ be a real $3\times 3$matrix representing a proper rotation. That is, $R$ satisfies the following properties: $$ \sum_j R_{ij}R_{kj} = \sum_j R_{ji}R_{jk} = \delta_{ik} \\ \sum_{ijk} \epsilon_{ijk} R_{i1}R_{j2}R_{k3} = 1 $$ that is, $R$ is orthogonal and has determinant $+1$.
It is "obvious" that if you apply the rotation to each vector and then take the cross product, you get the same result as if you take the cross product first and then apply the rotation to that.
$$
(R \vec{a}) \times (R \vec{b}) = R (\vec{a} \times\vec{b} )
$$
there are several ways to prove this, for example, see
Rotations preserve cross products?
.
What I would like to see is a proof that only uses index manipulations. In terms of indices, what we are trying to prove is that assuming $R$ satisfies the two conditions given above, $$ \sum_p R_{ip} \sum_{jk} \epsilon_{pjk} a_j b_k = \sum_{jk} \epsilon_{ijk}\sum_n R_{jn}a_n \sum_m R_{km}b_m $$
One thing I tried is multiplying both sides by $R_{qi}$ and summing over $i$, because I can undo that by left-multiplying by $R$. This comes down to needing to prove that $$ \sum_{qjk}\epsilon_{ijk} R_{qi} R_{jn} R_{km} = \epsilon_{qnm} $$ which looks tantalizingly close to the determinant = 1 condition, but I'm not quite there.
I'll use lower and upper indices to keep track of which indices are (respectively) covariant and contravariant, to make notationally clearer some later steps where we use the standard inner product to convert indices between the two types. In particular, we write a rotation $R$ as $R^i{}_j$ and the tensor defining the cross product operation as $\epsilon_{ij}{}^k$. In this notation, the desired property (after some relabeling) is $$R^q{}_k \epsilon_{ij}{}^k = \epsilon_{ij}{}^p R^i{}_m R^j{}_n .$$ Here we use the Einstein summation convention, that is, we implicitly sum over any index that appears twice (once lower and once upper).
We write the inverse of a rotation as $(R^{-1})^i{}_j = (R^\top)^i{}_j = R_i{}^j$, in which notation the first property becomes $R_j{}^k R^i{}_k = \delta^i{}_k$. Then, as you've already done, we can contract both sides of the previous display equation with $R_q{}^r$ to show that the desired identity is equivalent to $$\epsilon_{ij}{}^r = \epsilon_{mn}{}^q R^m{}_i R^n{}_j R_q{}^r .$$
Now, raising and lowering indices with the standard inner product $g$ gives the more symmetric equivalent identity $$\phantom{(\ast)} \qquad \epsilon_{ijk} = \epsilon_{mnp} R^m{}_i R^n{}_j R^p{}_k , \qquad (\ast)$$ which is just the assertion that volume form $\epsilon$ is invariant under pullback by rotations.
If we use (the inverse of) the inner product to contract any two of the free indices $i, j, k$ in $\epsilon_{mnp} R^m{}_i R^n{}_j R^p{}_k$, say, $i, j$, we find $$g^{ij} \epsilon_{mnp} R^m{}_i R^n{}_j R^p{}_k = - g^{ij} \epsilon_{nmp} R^n{}_i R^m{}_j R^p{}_k = -g^{ji} \epsilon_{nmp} R^n{}_j R^m{}_i R^p{}_k ,$$ which together with symmetry of notation gives that $\epsilon_{mnp} R^m{}_i R^n{}_j R^p{}_k$ is totally skew (this particular statement holds for all linear transformations, not just rotations). Thus, to verify $(\ast)$ it's enough to check that it holds for $(i, j, k) = (1, 2, 3)$. But this is exactly the content of the second of the given identities for $R$, namely, $\epsilon_{mnp} R^m{}_1 R^n{}_2 R^p{}_3 = 1$.