Proof
Courtesy of this thread from PhysicsForums.
\begin{align}
\nabla \times (\vec{A} \times \vec{B})
&=\partial_l \hat{e}_l \times (a_i b_j \hat{e}_k \epsilon_{ijk}) \\
&=\partial_l a_i b_j \epsilon_{ijk} \underbrace{ (\hat{e}_l \times \hat{e}_k)}_{(\hat{e}_l \times \hat{e}_k) = \hat{e}_m \epsilon_{lkm} } \\
&=\partial_l a_i b_j \hat{e}_m \underbrace{\epsilon_{ijk} \epsilon_{mlk}}_{\text{contracted epsilon identity}} \\
&=\partial_l a_i b_j \hat{e}_m \underbrace{(\delta_{im} \delta_{jl} - \delta_{il} \delta_{jm})}_{\text{They sift other subscripts}} \\
&=\partial_j (a_i b_j \hat{e}_i)- \partial_i (a_i b_j \hat{e}_j) \\
&=\color{blue}{a_i \partial_j b_j \hat{e}_i + b_j \partial_j a_i \hat{e}_i} - (\color{green}{a_i \partial_i b_j \hat{e}_j + b_j \partial_i a_i \hat{e}_j}) \\
&= \vec{A}(\nabla \cdot \vec{B}) + (\vec{B} \cdot \nabla)\vec{A} - (\vec{A} \cdot \nabla)\vec{B} - \vec{B}(\nabla \cdot \vec{A}) \\
\end{align}
Edit: As is pointed out by @enzotib, the blue and green sums are derivatives of products.
Why did the deltas vanish?
Due to Kronecker $\delta$'s sifting property. Recall the definition of Kronecker delta:
$$\delta_{ij}=\begin{cases}
0,\quad \text{if } i\ne j, \\
1,\quad \text{if } i=j.
\end{cases}$$
Thus, for $j\in\mathrm Z$:
$$\sum\limits_{-\infty}^{\infty}a_i\delta_{ij} = a_j$$
This is just like filtering (or sifting), because only when $i=j$ does $\delta_{ij} = 1$. Other terms are zeroes. This also works for partial derivatives.
For example,
$$\partial_l a_i b_j \hat{e}_m \delta_{im} \delta_{jl} = \partial_l a_i b_j \hat{e}_i \delta_{jl}$$
- If $l\ne j$ then $\partial_l a_i b_j \hat{e}_i \delta_{jl} = \partial_l a_i b_j \hat{e}_i(0) = 0$;
- If $l=j$ then $\partial_l a_i b_j \hat{e}_i \delta_{jl} = \partial_j a_i b_j \hat{e}_i (1) = \partial_j a_i b_j \hat{e}_i$.
Thus,
$$\partial_l a_i b_j \hat{e}_m \delta_{im} \delta_{jl} =\partial_j (a_i b_j \hat{e}_i).$$
Note that $\hat{e}_i$ is a const std::vector<int>.
Some thoughts
- There are two cross products (one of them is Curl) and we use different subscripts (of partials and Levi-Civita symbol to distinguish them, e.g., $l$ for the curl and $k$ for $\vec{A} \times \vec{B}$.
- We move the variables around quite often.
- The cross product of two basis is explained in the underbrace.
- The contracted epsilon identity is very useful. We replaced them by Kronecker $\delta$-s.
- Kronecker $\delta$-s are known to select things efficiently.
- In many proofs of vector calculus identities (this one included), we add and substract extra terms.
- Do I love Levi-Civita symbols and Einstein Notation? I'm ambivalent.
So, your first step is indeed correct, but the one you are not sure about is wrong. You are probably confused by your own notation, I prefer to write $\partial_i=\frac{\partial}{\partial x_i}$ instead of $\partial x_i$. And remember that repeated indices are being summed over, so e.g. $u_ju_j=\mathbf{u}\cdot\mathbf{u}$. Given that, you have
\begin{equation}
u_j\partial _i u_j-u_j\partial_j u_i
=(1/2)\partial_i(u_j u_j)-(\mathbf{u}\cdot \nabla)u_i
=(1/2) \partial_i (\mathbf{u}\cdot\mathbf{u})-(\mathbf{u}\cdot \nabla)u_i,
\end{equation}
which (taking into account your earlier work) is the $i$-th component of the equation
\begin{equation}
\mathbf{u}\times(\nabla\times \mathbf{u})=\frac{1}{2}\nabla(\mathbf{u}\cdot\mathbf{u})-(\mathbf{u}\cdot \nabla)\mathbf{u}
\end{equation}
Best Answer
I am not sure if this helps, but I just found that you can actually use the contracted epsilon identity to get a pretty good vector/dyadic representation of the chain rule. The nice thing about this particular identity (at least to me) is that it does what any good chain rule should do, and applies the operator in question to the argument.
Please allow me to preempt this by saying that I am not that familiar with the conventions of dyadic notation, but I will present what I have figured out in index notation form, so that if anyone wants to go in, and fix my notation, they will know how to. Anyway, here is what I found:
\begin{equation} \nabla \times \left(\mathbf{A} \circ \mathbf{B}\right) = -\left(\nabla_\mathbf{B}\cdot \mathbf{A}\right)\left(\nabla\times \mathbf{B}\right) + \left(\frac{\partial \mathbf{A}}{\partial \mathbf{B}}\right)^T\left(\nabla \times \mathbf{B}\right) + \left(\frac{\partial \mathbf{A}}{\partial \mathbf{B}}\right)^T \!\!\!\begin{array}{c} _\cdot \\ ^\times\end{array}\!\!\!\left(\nabla \mathbf{B}^T\right) \end{equation}
The vertical operator notation means that I am crossing the second index of the left tensor with the 1st index of the right tensor, and then contracting the 1st index of the left tensor with the 2nd index of the right tensor.
Working it out in index notation, we start with \begin{equation} \epsilon_{ijk}\frac{\partial A_j}{\partial B_l}\frac{\partial B_l}{\partial X_k} \end{equation}
What we want is to have the levi-civita symbol apply to B. The first step for this is to "free" the l and k indices. We can do this, using the kronecker delta
\begin{equation} \epsilon_{ijk}\frac{\partial A_j}{\partial B_l}\frac{\partial B_l}{\partial X_k} = \delta_{k_0k_1}\delta_{l_0l_1}\epsilon_{ijk_0}\frac{\partial A_j}{\partial B_{l_0}}\frac{\partial B_{l_1}}{\partial X_{k_1}} \end{equation}
Now, we can get to the contracted epsilon identity by adding appropriate cancelling terms to the RHS:
\begin{equation} \left(\delta_{k_0k_1}\delta_{l_0l_1}\epsilon_{ijk_0}\frac{\partial A_j}{\partial B_{l_0}}\frac{\partial B_{l_1}}{\partial X_{k_1}} - \delta_{k_0l_1}\delta_{l_0k_1}\epsilon_{ijk_0}\frac{\partial A_j}{\partial B_{l_0}}\frac{\partial B_{l_1}}{\partial X_{k_1}}\right) + \delta_{k_0l_1}\delta_{l_0k_1}\epsilon_{ijk_0}\frac{\partial A_j}{\partial B_{l_0}}\frac{\partial B_{l_1}}{\partial X_{k_1}} \end{equation}
We can then use the contracted epsilon identity to rewrite this as
\begin{equation} \epsilon_{rk_0l_0}\epsilon_{rk_1l_1}\epsilon_{ijk_0}\frac{\partial A_j}{\partial B_{l_0}}\frac{\partial B_{l_1}}{\partial X_{k_1}} + \delta_{k_0l_1}\delta_{l_0k_1}\epsilon_{ijk_0}\frac{\partial A_j}{\partial B_{l_0}}\frac{\partial B_{l_1}}{\partial X_{k_1}} \end{equation}
Now, we are halfway there, for $\epsilon_{rk_1l_1}$ applies to $\frac{\partial B_{l_1}}{\partial X_{k_1}}$
Notice also, that when we juxtapose $\epsilon_{rk_0l_0}$ and $\epsilon_{ijk_0}$ we can apply the contracted epsilon identity once again. Combining these two facts gives us \begin{equation} \left(-\delta_{ri}\delta_{l_0j} + \delta_{rj}\delta_{l_0i}\right)\frac{\partial A_j}{\partial B_{l_0}}\epsilon_{rk_1l_1}\frac{\partial B_{l_1}}{\partial X_{k_1}} + \delta_{k_0l_1}\delta_{l_0k_1}\epsilon_{ijk_0}\frac{\partial A_j}{\partial B_{l_0}}\frac{\partial B_{l_1}}{\partial X_{k_1}} \end{equation}
The expression in vector/dyadic notation above follows directly from this.