Coordinate free proof for $a\times (b\times c) = b(a\cdot c) - c(a\cdot b)$

Since $b\times c$ is normal to the plane $b,\,c$ span, $a\times (b\times c)$, which is orthogonal to this vector, is in said plane. The coefficients $B,\,C$ for which the result is $Bb+Cc$ are invariant under rotations, and clearly $B$ must be linear in $a,\,c$ while $C$ is linear in $a,\,b$, so constants $B',\,C'$ exist with $a\times (b\times c) =B' (a\cdot c) b + C' (a\cdot b) c$. Since the left-hand side is antisymmetric, $C'=-B'$. Since $B'$ must be a constant (since both sides are linear in each vector), we can use any vectors we like for which both sides are non-zero to compute $B'$. Example: $a=b=i,\,c=j$ so $a\times (b\times c) = i\times k = -j$ and $(a\cdot c) b - (a\cdot b) c = -j$ as required.

Okay, I really hate the sign issues with the Hodge star, so I am gonna assume that $\star\star=1$, the end result will be good.

The fundamental relationship is that for $k$-vectors/forms, we have $\langle\omega,\eta\rangle\mu=\omega\wedge\star\eta$, where the angle brackets are the inner product on the exterior algebra and $\mu$ is the volume form/multivector.

Let's start with $x$ being an arbitrary vector and taking a look at $$ \langle x,a\times(b\times c)\rangle\mu=\langle x,\star(a\wedge\star(b\wedge c))\rangle\mu= \\=x\wedge(a\wedge\star(b\wedge c))=(x\wedge a)\wedge\star(b\wedge c)= \\=\langle x\wedge a,b\wedge c\rangle\mu=\det\left(\begin{matrix}\langle x,b\rangle & \langle x,c\rangle \\ \langle a,b\rangle & \langle a,c\rangle\end{matrix}\right)\mu= \\=(\langle x,b\rangle \langle a,c\rangle-\langle x,c\rangle \langle a,b\rangle)\mu, $$ comparing the LHS with the RHS we can "divide" (ofc not really divide) by $\mu$ since the coefficients need to agree, and because $x$ was arbitrary, and the inner product is nondegenerate, we can "" divide "" by $x$ and we have $$ a\times(b\times c)=b\langle a,c\rangle-c\langle a,b\rangle. $$

Let $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$ be vector fields on $\mathbb{R}^{3}$ (we could extend to $\mathbb{R}^{n}$ if we wish!), considered as a Riemannian manifold equipped with metric $g$ and induced Hodge map $\star$. Let $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$ have corresponding vector field representations $U,V,W$ on $\mathbb{R}^{n}$ respectively. Then

$$\begin{align} \mathbf{a}\times(\mathbf{b}\times \mathbf{c}) \,\equiv\, \star(\widetilde{U} \wedge \star(\widetilde{V} \wedge \widetilde{W})) \end{align}$$

where $\widetilde{X}$ denotes the metric dual of $X$ (i.e. $\widetilde{X}=g(X,-)$) and the equivalence is up to metric dual. Then

$$\begin{align} \star(\widetilde{U} \wedge \star(\widetilde{V} \wedge \widetilde{W})) &\,=\, \star(\widetilde{U} \wedge i_{W}\star \widetilde{V}) \,=\, \star ( i_{W}\widetilde{U} \wedge \star \widetilde{V} - i_{W} (\widetilde{U}\wedge \star \widetilde{V})) \\ &\,=\, (i_{W}\widetilde{U})\star\star \widetilde{V} - \star i_{W}(\widetilde{U}\wedge\star \widetilde{V}) \\ &\,=\, g(U,W)\widetilde{V} - \star(\widetilde{U}\wedge\star \widetilde{V}) \wedge \widetilde{W} \\ &\,=\, g(U,W)\widetilde{V} - g(U,V)\widetilde{W} \end{align}$$

where $i_{X}$ denotes the interior derivative with respect to $X$ and we have used the identities:

$$\begin{align} \star\star \alpha &\,=\, \alpha \\[0.2cm] \star(\widetilde{X}\wedge\star\widetilde{Y}) &\,=\, g(X,Y) \\[0.2cm] \star(\alpha \wedge \widetilde{X}) &\,=\, i_{X}\star \alpha \\[0.2cm] \widetilde{X} \wedge \star\alpha &\,=\, (-1)^{p+1}\star i_{W}\alpha \end{align}$$

for any $p$-form $\alpha$ and vector fields $X,Y$ (note that the first and second identities are specific to $\mathbb{R}^{3}$). Then note that

$$g(U,W)\equiv \mathbf{a}\cdot\mathbf{c}\quad\text{and}\quad g(U,V)\equiv\mathbf{a}\cdot\mathbf{b}$$

and so then the result follows after taking the metric dual of our expression.