Proof for the curl of a curl of a vector field
Yes, there's a more elegant way! It uses the language of differential forms, which has replaced the 19th-century language of gradients, divergences, and curls in modern geometry. You can appreciate the simplicity of this language even before learning how to read it:
For any 1-form $A$, $$\begin{align*} (\star d)(\star d)A & = (\star d \star)dA \\ \operatorname{curl} \operatorname{curl} A & = d^\dagger dA \end{align*}$$ Recalling that $\Delta = d d^\dagger + d^\dagger d$, we see that $$\begin{align*} \operatorname{curl} \operatorname{curl} A & = -d d^\dagger A + \Delta A \\ & = d(\star d \star)A + \Delta A \\ & = \operatorname{grad} \operatorname{div} A + \Delta A \end{align*}$$
This is the identity you wanted to prove, where $-\Delta$ is the vector Laplacian.
My favorite place to learn about differential forms is in Chapters 4 and 5 of Gauge Fields, Knots, and Gravity by John Baez and Javier Muniain.
Here's a rough glossary that should help you move between the language of differential forms and the old language of vector calculus. I'll start by telling you the various kinds of differential forms, and the basic operations on them.
- In $n$-dimensional space, there are $n+1$ kinds of differential forms, from 0-forms to $n$-forms. You can think of a $k$-form as a $k$-dimensional density. A 0-form is a function, and a 1-form is a row-vector field (in coordinate-free language, a dual-vector field).
- The exterior derivative is an operation $d$ that turns $k$-forms into $(k + 1)$-forms. As its name suggests, it generalizes the operation of differentiating a function.
- The Hodge star is an operation $\star$ that turns $k$-forms in to $(n - k)$-forms. It comes from the dot product between column vectors. In fact, the Hodge star encodes the same geometric information as the dot product: if you know one, you can reconstruct the other.
- The codifferential is an operation $d^\dagger$ that turns $k$-forms into $(k - 1)$-forms. In an odd-dimensional space, like ordinary 3-dimensional space, applying $d^\dagger$ to a $k$-form is the same as applying $(-1)^k \star d \star$. In an even-dimensional space, $d^\dagger$ always acts like $-\star d \star$.
If you keep in mind that a 0-form is a function and a 1-form is a row-vector field, all the familiar operations of vector calculus can be written in terms of the ones above.
- The gradient of a function $f$ is the 1-form $df$.
- The curl of a 1-form $A$ is the 1-form $\star dA$.
- The divergence of a 1-form $A$ is the function $\star d \star A$.
- The Laplacian of a function or 1-form $\omega$ is $-\Delta \omega$, where $\Delta = dd^\dagger + d^\dagger d$. The operator $\Delta$ is often called the Laplace-Beltrami operator.
With this glossary in hand, you should be able to follow the steps of the calculation above, which is mostly just translating back and forth between languages. The only tricky bit is getting the sign right when you rewrite $d^\dagger$ as $\pm \star d \star$: you have to figure out what kind of form $d^\dagger$ is being applied to.
I usually just grind through these types of things with the Einstein notation. The notational rule is that a repeated index is summed over the directions of the space. So,
$$ x_i x_i = x_1^2+x_2^2+x_3^2.$$
A product with different indices is a tensor and in the case below has 9 different components,
$$ x_i x_j = \left( \begin{array} \ x_1^2 & x_1x_2 & x_1 x_3 \\ x_2 x_1 & x_2^2 & x_2 x_3 \\ x_3 x_1 & x_3x_2 & x_3^2 \\ \end{array} \right).$$
Since we are dealing with the curle we also need the levi-cevita tensor $\epsilon_{ijk}$. This tensor has the property that it is equal to +1 when $ijk$ is an even permutation of 123, -1 when $ijk$ is an odd permutation of 123, and zero otherwise.
$$ \epsilon_{ijk} = \begin{cases} 1 \qquad ijk=123,312,231 \\ -1 \qquad ijk=213,132,321 \\ 0 \qquad \text{otherwise} \end{cases} $$
We can use the Levi-Cevita tensor to write the cross product of two vectors as,
$$ (\vec{u} \times\vec{v} )_k = \epsilon_{ijk} u_i v_j. $$ We will also need the Kronecker delta, $\delta_{ij}$, which is like an identity matrix; it is equal to 1 if the indices match and zero otherwise.
$$ \delta_{ij} = \begin{cases} 1 \qquad i=j \\ 0 \qquad i \neq j \end{cases} $$
Now that we have these tools defined we start by writing the curl of a vector field,
$$ (\nabla \times \vec{A})_k = \epsilon_{ijk} \partial_i A_j$$
Now to get the curl of the curl we write,
$$ (\nabla \times \nabla \times \vec{A} )_k = \epsilon_{ijk} \partial_i (\nabla \times \vec{A})_j = \epsilon_{ijk} \partial_i \epsilon_{abj} \partial_a A_b$$
$$ =\epsilon_{ijk} \epsilon_{abj} \partial_i \partial_a A_b $$
Now we need to consider this product of Levi-Cevita Symbols, $ \epsilon_{ijk} \epsilon_{abj} $. It is possible to express this product in terms of Kronecker delta's,
$$\epsilon_{ijk} \epsilon_{abj} = \delta_{ib} \delta_{ka} - \delta_{ia} \delta_{kb},$$
I won't prove this here because its outside the scope of the question and it isn't difficult to check.
With this identity in hand we get a simpler form for $\nabla \times \nabla \times \vec{A}$,
$$ (\nabla \times \nabla \times \vec{A} )_k = (\delta_{ib} \delta_{ka} - \delta_{ia} \delta_{kb}) \partial_i \partial_a A_b $$
$$ = \delta_{ib} \delta_{ka}\partial_i \partial_a A_b - \delta_{ia} \delta_{kb}\partial_i \partial_a A_b $$
$$ = \partial_b \partial_k A_b - \partial_i \partial_i A_k $$
$$ = ( \nabla( \nabla \cdot \vec{A} ) - \nabla^2 \vec{A} )_k $$
which proves the identity.
This approach probably doesn't seem to simple to you if you aren't familiar with the tools I am using, but honestly this approach has really grown on me as the easiest way to get the result quickly. This identity comes up all the time in the theory of Electricity and Magnetism and this is the formalism that we use in Physics.
Use geometric calculus, the calculus of clifford algebra. Let $A$ be some $k$-vector field. If $k=1$, then this is a conventional vector field. If $k=0$, then this is a scalar field. If $k=2$, then $A$ is a "bivector field"--with each point, we associate some oriented plane and magnitude. Compare with the vector field case, in which we associate with each point an oriented line (direction) and magnitude.
Now, define the action of $\nabla$:
$$\nabla A = \nabla \cdot A + \nabla \wedge A$$
where $\nabla \cdot A$ is a $k-1$ vector field, and $\nabla \wedge A$ is a $k+1$ vector field. If $k=0$, then $\nabla A = \nabla \wedge A$ by definition. Similarly, if $k=n$, then $\nabla A = \nabla \cdot A$.
Useful identity: $\nabla \cdot (\nabla \cdot A) = 0$ and $\nabla \wedge (\nabla \wedge A) = 0$. Both of these follow from the equality of mixed partial derivatives.
Thus, you can infer that
$$\nabla^2 A = \nabla (\nabla \cdot A + \nabla \wedge A) = \nabla \cdot (\nabla \cdot A) + \nabla \wedge (\nabla \cdot A) + \nabla \cdot (\nabla \wedge A) + \nabla \wedge (\nabla \wedge A)$$
The first and last terms are, as we just established, zero, and we're left with
$$\nabla^2 A = \nabla \cdot (\nabla \wedge A) + \nabla \wedge (\nabla \cdot A)$$
All well and good, but how do we connect this back to vector calculus? Relate these dot and wedge products to the cross product like so:
$$a \wedge b = i(a \times b)$$
where $i$ is the right-handed unit trivector. This allows us to relate $\nabla \times A$, for $A$ a vector field, to $\nabla \wedge A$ like so:
$$\nabla \wedge A = i (\nabla \times A)$$
Also note a similar identity involving a vector $a$ and a bivector $B$:
$$i (a \cdot B )= a \wedge (iB)$$
$iB$ is always a vector. This allows us to turn dots to crosses, and the identity becomes
$$\begin{align*}\nabla^2 A &= \nabla \cdot (\nabla \wedge A) + \nabla \wedge (\nabla \cdot A) \\ &= -i \nabla \wedge (i \nabla \wedge A) + \nabla (\nabla \cdot A) \\ &= -\nabla \times (\nabla \times A) + \nabla (\nabla \cdot A)\end{align*}$$
And you're done.