How to Prove $(A \times B) \otimes C + (B \times C) \otimes A + (C \times A) \otimes B = [(A \times B) \cdot C] \bf{I}$? Question
Assume that $A$, $B$, $C$, and $D$ are four vectors in $\mathbb{R}^3$. Then I want to show that

$$ {\bf{M}} \equiv (A \times B) \otimes C + (B \times C) \otimes A + (C \times A) \otimes B = [(A \times B) \cdot C] {\bf{I}} \tag{1}$$

where $\otimes$ is the tensor product, $\cdot$ is the scalar product and $\times$ is the cross product. Also, ${\bf{I}}$ is the second order identity tensor.

Motivation
The motivation behind this question is the vector identity
$$ [(A \times B) \cdot C] D = (A \times B) (C \cdot D) + (B \times C) (A \cdot D) + (C \times A) (B \cdot D) \tag{2}$$
which can be seen as
$$ [(A \times B) \cdot C] D = {\bf{M}}D  \tag{3}$$
So it seems that proving $(1)$ is equivalent to proving $(2)$. I can prove $(2)$ directly which is given in this post. But I cannot prove $(1)$ directly.

My work
Some computations using Summation Convention and the properties of Kronecker's Delta and Permutation Symbol lead to
$$\eqalign{
  & \,\,\,\,\,{\left[ {(A \times B) \otimes C + (B \times C) \otimes A + (C \times A) \otimes B} \right]_{km}}  \cr 
  &  = {A_i}{B_j}{C_m}{\varepsilon _{ijk}} + {B_i}{C_j}{A_m}{\varepsilon _{ijk}} + {C_i}{A_j}{B_m}{\varepsilon _{ijk}}  \cr 
  &  = {A_i}{B_j}{C_n}{\delta _{nm}}{\varepsilon _{ijk}} + {B_i}{C_j}{A_n}{\delta _{nm}}{\varepsilon _{ijk}} + {C_i}{A_j}{B_n}{\delta _{nm}}{\varepsilon _{ijk}}  \cr 
  &  = {A_i}{B_j}{C_n}{\delta _{nm}}{\varepsilon _{ijk}} + {B_j}{C_n}{A_i}{\delta _{im}}{\varepsilon _{jnk}} + {C_n}{A_i}{B_j}{\delta _{jm}}{\varepsilon _{nik}}  \cr 
  &  = \left( {{\delta _{nm}}{\varepsilon _{ijk}} + {\delta _{im}}{\varepsilon _{jnk}} + {\delta _{jm}}{\varepsilon _{nik}}} \right){A_i}{B_j}{C_n}  \cr 
  &  = \left( {{\delta _{nm}}{\varepsilon _{ijk}} - {\delta _{im}}{\varepsilon _{njk}} - {\delta _{jm}}{\varepsilon _{ink}}} \right){A_i}{B_j}{C_n}  \cr 
  & \mathop  = \limits^{???} \left( {{\varepsilon _{ijn}}{\delta _{km}}} \right){A_i}{B_j}{C_n}  \cr 
  &  = \left( {{A_i}{B_j}{C_n}{\varepsilon _{ijn}}} \right){\delta _{km}}  \cr 
  &  = \left[ {\left( {A \times B} \right) \cdot C} \right]{\delta _{km}}  \cr 
  &  = {\left\{ {\left[ {\left( {A \times B} \right) \cdot C} \right]{\bf{I}}} \right\}_{km}} \cr} $$
And the only thing remains to prove is
$${\delta _{nm}}{\varepsilon _{ijk}} - {\delta _{im}}{\varepsilon _{njk}} - {\delta _{jm}}{\varepsilon _{ink}} = {\varepsilon _{ijn}}{\delta _{km}} \tag{4}$$
Is there an elegant way to prove $(4)$ or $(1)$ directly?
 A: Linear Independent Case
A nice way of proving this is to use non-orthogonal basis for $\mathbb{R}^3$. Hence, consider the following definitions for the non-orthogonal basis
$$\matrix{
   {{{\bf{g}}_1} = A} & {{{\bf{g}}_2} = B} & {{{\bf{g}}_3} = C}  \cr 
 } \tag{1}$$
and then the dual basis will be
$$\matrix{
   {V = \left( {{{\bf{g}}_1} \times {{\bf{g}}_2}} \right) \cdot {{\bf{g}}_3}} \hfill & {} \hfill & {} \hfill  \cr 
   {{{\bf{g}}^1} = {{{{\bf{g}}_2} \times {{\bf{g}}_3}} \over V},} \hfill & {{{\bf{g}}^2} = {{{{\bf{g}}_3} \times {{\bf{g}}_1}} \over V},} \hfill & {{{\bf{g}}^3} = {{{{\bf{g}}_1} \times {{\bf{g}}_2}} \over V}} \hfill  \cr 
   {{{\bf{g}}^k} = {1 \over {2V}}{\varepsilon ^{ijk}}{{\bf{g}}_i} \times {{\bf{g}}_j}} \hfill & {} \hfill & {} \hfill  \cr 
   {{{\bf{g}}_i} \times {{\bf{g}}_j} = V{\varepsilon _{ijk}}{{\bf{g}}^k}} \hfill & {} \hfill & {} \hfill  \cr 
 } \tag{2}$$
next, we write the second order identity tensor in the basis ${{\bf{g}}^i} \otimes {{\bf{g}}_j}$ as follows
$${\bf{I}} = I_i^j{{\bf{g}}^i} \otimes {{\bf{g}}_j} \tag{3}$$
then we dot product by ${{\bf{g}}_m}$ from left and by ${{\bf{g}}^n}$ from right and use the property ${{\bf{g}}^i} \cdot {{\bf{g}}_j} = {{\bf{g}}_i} \cdot {{\bf{g}}^j} = \delta _j^i$  to obtain
$${{\bf{g}}_m}{\bf{I}}{{\bf{g}}^n} = I_m^n = \delta _m^n \tag{4}$$
Now, we combine $(2)$, $(3)$, and $(4)$ to get
$${\bf{I}} = I_i^j{{\bf{g}}^i} \otimes {{\bf{g}}_j} = \delta _i^j{{\bf{g}}^i} \otimes {{\bf{g}}_j} = {{\bf{g}}^j} \otimes {{\bf{g}}_j} = {{\bf{g}}^k} \otimes {{\bf{g}}_k}{\mkern 1mu} {\mkern 1mu} {\mkern 1mu}  = {1 \over V}\left( {{1 \over 2}{\varepsilon ^{ijk}}{{\bf{g}}_i} \times {{\bf{g}}_j}} \right) \otimes {{\bf{g}}_k}{\rm{ }} \tag{6}$$
or equivalently
$$\left[ {\left( {{{\bf{g}}_1} \times {{\bf{g}}_2}} \right) \cdot {{\bf{g}}_3}} \right]{\bf{I}} = \left( {{{\bf{g}}_1} \times {{\bf{g}}_2}} \right) \otimes {{\bf{g}}_3} + \left( {{{\bf{g}}_2} \times {{\bf{g}}_3}} \right) \otimes {{\bf{g}}_1} + \left( {{{\bf{g}}_3} \times {{\bf{g}}_1}} \right) \otimes {{\bf{g}}_2} \tag{7}$$
which is the relation $(1)$ in the question with a slightly different notation. 
Linear Dependent Case
Consider the case when  the vectors $\bf{g}_1$, $\bf{g}_2$, and $\bf{g}_3$ are linearly dependent and hence they cannot form a basis for $\mathbb{R}^3$. This means that there exists real numbers $a$ and $b$ not both zero such that
${{\bf{g}}_1} = a{{\bf{g}}_2} + b{{\bf{g}}_3}$
will hold. This will lead to $V=0$ and turns equation $(7)$ into a trivial identity
$$\eqalign{
  & V{\bf{I}} = \left( {{{\bf{g}}_1} \times {{\bf{g}}_2}} \right) \otimes {{\bf{g}}_3} + \left( {{{\bf{g}}_2} \times {{\bf{g}}_3}} \right) \otimes {{\bf{g}}_1} + \left( {{{\bf{g}}_3} \times {{\bf{g}}_1}} \right) \otimes {{\bf{g}}_2}  \cr 
  & {\bf{0}} = \left( {\left( {a{{\bf{g}}_2} + b{{\bf{g}}_3}} \right) \times {{\bf{g}}_2}} \right) \otimes {{\bf{g}}_3} + \left( {{{\bf{g}}_2} \times {{\bf{g}}_3}} \right) \otimes \left( {a{{\bf{g}}_2} + b{{\bf{g}}_3}} \right) + \left( {{{\bf{g}}_3} \times \left( {a{{\bf{g}}_2} + b{{\bf{g}}_3}} \right)} \right) \otimes {{\bf{g}}_2}  \cr 
  & {\bf{0}} = b\left( {{{\bf{g}}_3} \times {{\bf{g}}_2}} \right) \otimes {{\bf{g}}_3} + a\left( {{{\bf{g}}_2} \times {{\bf{g}}_3}} \right) \otimes {{\bf{g}}_2} + b\left( {{{\bf{g}}_2} \times {{\bf{g}}_3}} \right) \otimes {{\bf{g}}_3} + a\left( {{{\bf{g}}_3} \times {{\bf{g}}_2}} \right) \otimes {{\bf{g}}_2}  \cr 
  & {\bf{0}} =  - b\left( {{{\bf{g}}_2} \times {{\bf{g}}_3}} \right) \otimes {{\bf{g}}_3} - a\left( {{{\bf{g}}_3} \times {{\bf{g}}_2}} \right) \otimes {{\bf{g}}_2} + b\left( {{{\bf{g}}_2} \times {{\bf{g}}_3}} \right) \otimes {{\bf{g}}_3} + a\left( {{{\bf{g}}_3} \times {{\bf{g}}_2}} \right) \otimes {{\bf{g}}_2}  \cr 
  & {\bf{0}} = {\bf{0}} \cr} \tag{8}$$
Conclusion
Relation $(7)$ will hold for every three vectors $\bf{g}_1$, $\bf{g}_2$, and $\bf{g}_3$ in $\mathbb{R}^3$. Finally, replacing ${{\bf{g}}_1}$ by ${{\bf{g}}_i}$, ${{\bf{g}}_2}$ by ${{\bf{g}}_j}$, and ${{\bf{g}}_3}$ by ${{\bf{g}}_k}$ in $(7)$ and doing some further computations will prove identity $(4)$ in the question.
A: This answer uses representation theory. Let $V=\mathbb{R}^3$ be the standard representation of $\mathrm{SO}(3)$, and suppose that $T:V\to V$ is a $\mathrm{SO}(3)$-equivariant map. Schur's lemma would say that $T$ must be a scalar multiple of the identity because $V$ is irreducible, but this is a real representation rather than a complex one so we can't say that automatically.
Instead, observe $T$'s characteristic polynomial is cubic, so it has a real zero, so for some $\lambda$ the operator $T-\lambda\cdot\mathrm{Id}$ has nontrivial kernel containing some line $L$ on which $T$ acts by the scalar $\lambda$. Pick rotations $R,S\in\mathrm{SO}(3)$ so that $V=L\oplus RL\oplus SL$, then $T(R\ell)=RT\ell=\lambda(R\ell)$ for $\ell\in L$ tells us that $T$ acts by the scalar $\lambda$ on $RL$, and similarly for the last line $SL$, so we conclude $T=\lambda\cdot\mathrm{Id}$.
Using the dot product we can identify $V$ with its dual space $V^*$. If $v\in V$, the corresponding dual vector $v^*$ is defined by $v^*(w):=v\cdot w$. In general $V\otimes V^*$ may be identified with $\mathrm{End}(V)$, as representations of $\mathrm{SO}(3)$.
Denote $T_{abc}=(a\times b)\otimes c^*+(b\times c)\otimes a^*+(c\times a)\otimes b^*\in\mathrm{End}(V)$. This defines an antisymmetric trilinear map $V\times V\times V\to\mathrm{End}(V)$, which induces a linear map $\bigwedge^3V\to\mathrm{End}(V)$. In fact notice this map is $\mathrm{SO}(3)$-equivariant. Since $\dim\bigwedge^3V=1$, the map $\mathrm{SO}(3)\to\mathrm{GL}(\bigwedge^3 V)\cong\mathbb{R}^\times$ must be trivial because $\mathrm{SO}(3)$ is compact and connected, so $\bigwedge^3V$ is the trivial representation and its image in $\mathrm{End}(V)$ must be contained in the space of $\mathrm{SO}(3)$-equivariant endomorphisms, so $T=\lambda\cdot\mathrm{Id}$ for some $\lambda$.
Using the fact $\mathrm{tr}(v\otimes w^*)=v\cdot w$, we see 
$$3\lambda=\mathrm{tr}(T_{abc})=(a\times b)\cdot c+(b\times c)\cdot a+(c\times a)\cdot b=3\,(a\times b)\cdot c$$
and therefore $T_{abc}=((a\times b)\cdot c)\,\mathrm{Id}$.
A: Using clifford algebra, the identity (4) can be rewritten as (using vectors $m, n, i, j, k$)
$$(m \cdot i)(j \wedge k \wedge n) + (m \cdot j)(i \wedge k \wedge n) + (m \cdot k)(i \wedge j \wedge n) - (m \cdot n)(i \wedge j \wedge k) = 0$$
$m$ is arbitrary, and therefore it can be canceled to yield
$$i(j \wedge k \wedge n) + j (i \wedge n \wedge k) - k(j \wedge i \wedge n) - n(i \wedge j \wedge k) = 0$$
As a geometric identity, this reads, "Given a set of four vectors, form a volume from each combination of three. Then, find the plane orthogonal to the fourth and weight it such that that plane and that vector span the volume with the magnitude of that volume. The four such oriented, weighted planes formed from each combination add to zero."
We can solve this problem (prove the identity) by using clifford products and grade projection. Associativity of the product operation allows us to quickly derive true relations from grade projection, and these relations involve the terms that appear in the desired identity.
We'll do this:
$$\begin{align*}
\langle ijkn \rangle_2 &= \color{red}{i(j \wedge k \wedge n)} + i \wedge [j (k \cdot n) - k (j \cdot n) + n(j \cdot k)] \\
&= \color{red}{(i \wedge j \wedge k)n} + [i(j \cdot k) + k (i \cdot j) - j(i \cdot k)] \wedge n\\
\langle jink \rangle_2 &= \color{red}{j (i \wedge n \wedge k)} + j \wedge [i(n \cdot k) + k (i \cdot n) - n(i \cdot k)] \\
&= \color{red}{(j \wedge i \wedge n) k} + [j(i \cdot n) + n(j \cdot i) - i(j \cdot n)] \wedge k
\end{align*}$$
The requisite terms for identity (4) are in red (mind that with the wedge products, there are apparent, but not substantial, differences in minus signs).
Rearrange these identities to yield equations of the form
$$\begin{align*}
i(j \wedge k \wedge n) - (i \wedge j \wedge k)n + \ldots &= 0 \\
j(i \wedge n \wedge k) - (j \wedge i \wedge n)k + \ldots &= 0
\end{align*}$$
Then add those two equations to arrive at the desired identity.  All the other terms (...) will cancel: the $i \wedge n$ terms and $j \wedge k$ terms cancel during the subtraction; the other terms cancel during the addition.
A: Rewriting:
$${\delta _{nm}}{\varepsilon _{ijk}} - {\delta _{im}}{\varepsilon _{njk}} - {\delta _{jm}}{\varepsilon _{ink}} = {\varepsilon _{ijn}}{\delta _{km}}$$
as:
$${\delta _{nm}}{\varepsilon _{ijk}} = {\delta _{im}}{\varepsilon _{njk}} + {\delta _{jm}}{\varepsilon _{ink}} + {\delta _{km}}{\varepsilon _{ijn}}$$
makes it easier to see the symmetry. Assuming $n,m,i,j,k\in\{1,2,3\}$ we can easily prove the case for $i=j=k$, namely that every $\varepsilon=0$.
Assume two of $i,j,k$ are equal, say $i=j$. Then the left-hand side is $0$, and the $3^{rd}$ term on the right-hand side is zero. But we also have that:
$${\delta _{im}}{\varepsilon _{njk}} + {\delta _{jm}}{\varepsilon _{ink}}$$
$$= {\delta _{im}}{\varepsilon _{njk}} + {\delta _{im}}{\varepsilon _{jnk}}$$
$$= {\delta _{im}}{\varepsilon _{njk}} - {\delta _{im}}{\varepsilon _{njk}}$$
$$=0$$
And so the equations holds for this case as well.
So we need to consider $i,j,k$ distinct.
If $n\ne m$, we have the LHS disappearing, and only one term to consider on the RHS, for example if $i=m$ so that $\delta_{im}=1$. However as $n\ne m$ and the $i,j,k$ are distinct, we have that either $j=n$ or $k=n$ and so this term is $0$ too.
Finally, if $n=m$, the LHS is $\varepsilon_{ijk}$ and if again for example $i=m$ then we have:
$${\delta _{im}}{\varepsilon _{njk}}$$
$$=\varepsilon _{njk}$$
But $n=m=i$ so:
$$=\varepsilon _{ijk}$$
A: The following argument just occured to me. (I probably half-remembered it from Cvitanovic's Birdtrack book.)
We use the relation
$$\varepsilon_{iab}\varepsilon_{ipq}=(\delta_{ap}\delta_{bq}-\delta_{aq}\delta_{bp})$$
to expand
$$\varepsilon_{ani}\varepsilon_{abm}\varepsilon_{bjk}$$
in two different ways.
First
$$\begin{align*}
\varepsilon_{ani}\varepsilon_{abm}\varepsilon_{bjk}&=(\delta_{nb}\delta_{im}-\delta_{nm}\delta_{ib})\varepsilon_{bjk}\\
&=\delta_{im}\varepsilon_{njk}-\delta_{nm}\varepsilon_{ijk}.
\end{align*}$$
Second
$$\begin{align*}
\varepsilon_{ani}\varepsilon_{abm}\varepsilon_{bjk}&=\varepsilon_{ani}(\delta_{mj}\delta_{ak}-\delta_{mk}\delta_{aj})\\
&=\delta_{mj}\varepsilon_{kni}-\delta_{mk}\varepsilon_{jni}.
\end{align*}$$
Equating these gives
$$\delta_{im}\varepsilon_{njk}-\delta_{nm}\varepsilon_{ijk}=\delta_{mj}\varepsilon_{kni}-\delta_{mk}\varepsilon_{jni}$$
and rearranging gives $(4)$.
