Proof for the curl of a curl of a vector field For a vector field $\textbf{A}$, the curl of the curl is defined by
$$\nabla\times\left(\nabla\times\textbf{A}\right)=\nabla\left(\nabla\cdot\textbf{A}\right)-\nabla^2\textbf{A}$$
where $\nabla$ is the usual del operator and $\nabla^2$ is the vector Laplacian.
How can I prove this relation? I tried the ugly/unefficient/brute-force method, by getting an expression for the LHS and the RHS for an arbitrary vector field $$\textbf{A}=\left(a(x,y,z),b(x,y,z),c(x,y,z)\right)$$
It does work (duh), but is there a more elegant way of doing this? Using matrix notation maybe?

EDIT: I got very good answers, from various perspectives. I would say @Spencer's derivation is the one I was looking for, using Einstein notation - and as a physics student, this was very helpful. However, @Vectornaut's solution not only is short and elegant, but it also introduced me to a whole new range of mathematics - and as a theoretical physics student, I appreciate learning new mathematical theories and trying to see how we can use them in physics.
 A: 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.
A: 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.
A: Let's actually follow suggestion by @idm and use the Grassmann identity to make a proof.
$$a\times (b\times c) = b (a \cdot c) - (a\cdot b) c.\qquad (1)$$
Let us assume, that we know, that it is true for $a$, $b$, $c$ being vectors with real numbers as their elements. The rigorous proof of the identity
$$\nabla \times (\nabla \times A)=\nabla (\nabla \cdot A)-\nabla^2 A\qquad (2)$$
can be constructed from (1) by looking the basic facts about how polynomials work.
Polynomial $p$ in variables $x_1,\dots,x_n$ is a formal expression of the form
$p=\sum_{k_1,\dots,k_n} a_{k_1,\dots,k_n} x_1^{k_1}\dots x_n^{k_n}$, where there are only finitely many terms with nonzero $a_{k_1,\dots,k_n}$. In each sum $k_i$ goes over all nonnegative integers.
The terms in this sum are called its monomials, and real numbers $a_{k_1,\dots,k_n}$ are called coefficients. Denote with $p|_{x_n=b}$ the polynomial, obtained from $p$ by replacing $x_n$ with a real number $b$ (so it becomes a polynomial in $n-1$ variables).
Claim 1. Let $p$ be a polynomial in $x_1,\dots,x_n$, and let $p|_{x_n=b}=0$. Then $p=(x_n-b)q$ for some polynomial $q$.
Proof. Replace $x_n$ with $y+b$ in $p$. You will get a polynomial $\widetilde p$ in $x_1,\dots,x_{n-1},y$. By assumption $\widetilde p|_{y=0}=0$, so all monomials of $\widetilde p$ have now a nonzero power of $y$. Thus we can write $\widetilde p = y \widetilde q$ and replace back $y$ with $x_n-b$. $\square$
Claim 2. Let $p$ be a polynomial in $x_1,\dots,x_n$, which is not 0 (i.e. it has at least 1 nonzero coefficient). Then it is not zero at some point $x_1=a_1,x_2=a_2,\dots,x_n=a_n$.
In other words, nonzero polynomial is nonzero as a function.
Proof.
You prove this by induction in $n$. For $n=0$ the claim is trivial.
Induction step. Let $d$ be the degree of $p$ with respect to variable $x_n$ (i.e. the highest k_n in the term $ax_1^{k_1}\dots x_n^{k_n}$ with nonzero $a$, appearing in $p$). We want to find $a_n$, such that after substituting it for $x_n$ we get a nonzero polynomial of $x_1,\dots,x_{n-1}$. Then we can apply the induction hypothesis to find $a_1,\dots,a_{n-1}$. Assume to the contrary, that we can't find such $a_n$. Then for every real number $b$ we get $p|_{x_n=b}=0$. Find any $d+1$ different real numbers $b_0,\dots,b_d$. By applying claim 1 $d+1$ times to $p$, write $p=(x_n-b_0)(x_n-b_1)\dots(x_n-b_d)q$. By assumption $q$ is nonzero (otherwise $p$ will be zero). So degree of $p$ in $x_n$ is at least $d+1$. Contradiction. $\square$
Claim 3. If we open the brackets in each of the 3 components of (1) (by only using definitions of $\times$ and $\cdot$, distributive law, commutativity and associativity of multiplication and addition), then we will get the same expression on the left hand side, as on the right.
Proof. Suppose for one of the 3 components we get expression $L$ on the left, and $R$ on the right, and these are not the same. Then $L-R$ is a nontrivial polynomial in 9 variables. But according to Grassmann identity, it is zero for all the values of these 9 variables. So by claim 2 it should be identically zero. $\square$
Derivatives $\partial/\partial x$, $\partial/\partial y$, $\partial/\partial z$ and components of $A$ satisfy all the properties we needed to show (1), i.e. all the properties, listed in claim 3, except we can't permute components of $A$ with derivatives.
But note, that while expanding $(1)$ we can always keep $c$ on the far right in each product. So we don't need to use associativity of multiplication to permute components of $c$ with components of $a$ and $b$.
So we can apply (1) to show (2) rigorously.
A: 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. 
A: This is an intermediate answer, using symmetry and a little of brute force... the advantage is that this method is -I think- reachable for a (smart?) student at a vector calculus course.

*

*(the hardest part?) convince yourself that every differential operator in the formula is rotation-invariant. Actually, one can use another simpler version: the differential operators in the formula are invariant under simultaneous chance of variables and coordinates: you can permute permute variables and compute gradient, and permute back variables/coordinates and you will get the same. Similar for divergence (it is actually a dual computation). For curl, you get a sign depending on the sign of the permutation, but you need to compute the curl twice, so you are done. Notice that it is enough to show the cases $1\leftrightarrow 2$ and $1\leftrightarrow 3$.


*The formula is $\mathbb R$-linear on $A$, so, you need to show it for
$A=(a,0,0)$, $A=(0,b,0)$ and $A=(0,0,c)$. But from 1), it is enough to prove only one of this possibilities.


*use brute force to check the formula for $A=(a(x,y,z),0,0)$. It is notably  simpler than brute force for the general case.
A: You can apply any vector identity to $\nabla$ so long as you keep track of what you are differentiating. To this end, we write
$$
  \nabla\times(\nabla\cdot A) = \check\nabla\times(\dot\nabla\times\check{\dot A}).
$$
Whether we write $\check{\dot A}$ or $\dot{\check A}$ is immaterial because mixed partials commute. Now we know for any vectors $a, b, c$ that
$$
  a\times(b\times c) = (a\cdot c)b - (a\cdot b)c
$$
thus
$$
  \check\nabla\times(\dot\nabla\times\check{\dot A})
    = (\check\nabla\cdot\check{\dot A})\dot\nabla - (\check\nabla\cdot\dot\nabla)\check{\dot A}.
$$
In the more traditional notation where $\nabla$ only differentiates directly to its right:
$$
  \nabla\times(\nabla\cdot A) = \nabla(\nabla\cdot A) - (\nabla\cdot\nabla)A = \nabla(\nabla\cdot A) - \nabla^2A.
$$
This "only differentiate to the right" convention is exactly why people run into trouble manipulating $\nabla$ expressions. In single variable calculus, everything is commutative so you can always make this happen; but in vector calculus, you cannot arbitrarily commute everything so that $\nabla$ is directly to the left of what it differentiates because we have non-commutative operations like $\times$.

This idea isn't really well demonstrated by the above identity, so I would like to give a better example. Consider
$$
  \nabla\times(A\times B).
$$
What happens is that people try to naively apply the identity
$$
  a\times(b\times c) = (a\cdot c)b - (a\cdot b)c
$$
and then arrive at
$$
  \nabla\times(A\times B) \not= (\nabla\cdot A)B - (\nabla\cdot B)A.
$$
But of course this is wrong! In the LHS, $\nabla$ is differentiating both $A$ and $B$; how did it magically only differentiate one or the other in the RHS? The correct statement is
$$
  \dot\nabla\times(\dot A\times\dot B) = (\dot\nabla\cdot\dot A)\dot B - (\dot\nabla\cdot\dot B)\dot A.
\tag{$*$}
$$
$\nabla$ is differentiating both $A$ and $B$ in both expressions on the RHS. Now, we need one more general fact: the derivative of an expression is the sum of the derivatives of its subexpressions. To illustrate:
$$
  \dot\nabla\times(\dot A\times\dot B)
    = \dot\nabla\times(\dot A\times B) + \dot\nabla\times(A\times\dot B).
$$
We could go from here and derive the identity we want, or we could apply the subexpression rule to the RHS of ($*$):
$$
  (\dot\nabla\cdot\dot A)\dot B = (\dot\nabla\cdot\dot A)B + (\dot\nabla\cdot A)\dot B,
$$$$
  (\dot\nabla\cdot\dot B)\dot A = (\dot\nabla\cdot\dot B)A + (\dot\nabla\cdot B)\dot A.
$$
All together,
$$
  \dot\nabla\times(\dot A\times\dot B)
    = (\dot\nabla\cdot\dot A)B + (\dot\nabla\cdot A)\dot B
      - (\dot\nabla\cdot\dot B)A - (\dot\nabla\cdot B)\dot A.
$$
Returning to more traditional notation where $\nabla$ only differentiates directly to its right,
$$
  \nabla\times(A\times B)
    = (\nabla\cdot A)B + (A\cdot\nabla)B
      - (\nabla\cdot B)A - (B\cdot\nabla)A,
$$
exactly as it should be.
