Tricky operator identity: $[L^2,[L^2,\vec{r}]]=2 \hbar ^2 \{ L^2, \vec{r}\}$? This operator identity showed up in a course I was taking, and it was given without proof.
$$[L^2,[L^2,\vec{r}]]=2 \hbar ^2 \{ L^2, \vec{r}\}$$
The curly brackets denote the anticommutator, $AB+BA$. The $\vec{r}$ operator is the position operator. The $L^2$ operator is given by: 
$$L^2 = -\hbar ^2 \left( \frac{1}{\sin \theta}  {\partial\over\partial\theta} (\sin \theta {\partial\over\partial\theta}) + \frac{1}{\sin^2  \theta} {\partial^2\over\partial\phi^2}\right)$$ 
Is there a way of proving this identity without tediously expanding all the commutators? I've been trying to find one but was unable to.
 A: The symbol $r$ in the identity represents (and will represent in the text below) the whole three-component vector of operators $\hat{\vec r} = (\hat x, \hat y, \hat z)$.
The simple way I found to prove the identity is to verify that all matrix elements of both sides match. Let's calculate the matrix elements of the operators $LHS,RHS$ between
$$\langle j,m,a| LHS| k,n,b\rangle$$
and similarly  for the right hand side. Here, $j,m$ and $k,n$ are the usual total angular momenta (which I will assume to be integers, just the orbital angular momentum case) and the $z$-component and $a,b$ represent the other quantum numbers that won't matter.
The advantage is that $\vec L$ combine to $L^2$ almost everywhere. The left hand side operator is
$$ L^2 L^2 r - 2 L^2 r L^2 + r L^2 L^2 $$
so the matrix element (because $L^2$ acts either on the bra or ket vector in a simple way) is the same as the matrix element of
$$ \hbar^4 r[ j(j+1)j(j+1)  - 2j(j+1)k(k+1) + k(k+1)k(k+1)] $$
The coefficient in the parenthesis is equal to a complete square,
$$ \hbar^4 r [j(j+1)-k(k+1)]^2 $$
Note that $\hbar^4 r$ is in all terms. The right hand side has the same matrix elements as the operator
$$ 2\hbar^4 r [j(j+1) + k(k+1)] $$
They don't look "obviously" equal: one is quartic, one is quadratic. But we must realize that the operators on both sides are $j=1$ vector operators, from the $\vec r$ factor, so they only change the angular momentum by zero or $\pm 1$.
So it is enough to compare the expressions for these three choices; for higher changes of $j$, the matrix elements on both sides clearly vanish (and are therefore equal). For $j=k$, the matrix element vanishes because of parity: $r$ carries the negative parity while the parities $(-1)^l$ are $(-1)^j$ or $(-1)^k$ for the bra/ket vectors.
For $j=k+1$, the LHS is
$$\hbar^4 r (k+1)^2 (k+2 - k)^2 = 4\hbar^2 r (k+1)^2 $$
while the RHS is
$$2\hbar^4 r[(k+1)(k+2)+k(k+1)]= 4\hbar^4 r(k+1)^2$$
so it works. The same verification applies to the case $k=j+1$, too, just $j,k$ are interchanged.
There are many other ways to calculate or verify the identity but I found this one easiest. Note that I am not assuming any coordinates; the abstract calculation above works in any coordinates.
A: I) For the record, here is the operator calculation that OP wants to avoid. The benefit of the calculation is that the operators are not sandwich with any bra/ket representation, and hence we do not have to worry about whether the bra/ket representation is faithful. Let us put $\hbar=1$ for simplicity. The starting point is the CCR
$$ [x^i, p_j]~=~\mathrm{i}\delta^i_j .\tag{1}$$
The CCR (1) ensures that the definition of the orbital angular momentum operator
$$ L_i~:=~\varepsilon_{ijk}~ x^jp_k~\stackrel{(1)}{=}~-\varepsilon_{ijk}~ p_jx^k \tag{2}$$
is Hermitian and that it does not suffer from operator ordering ambiguities. In particular, the position and the angular momentum are mutually perpendicular as operators
$$ x^i L_i~\stackrel{(2)}{=}~0~\stackrel{(2)}{=}~L_i x^i \tag{3}.$$
Einstein's summation convention is implicitly assumed everywhere in this answer.
II) Now let us calculate the LHS of OP's identity.
$$ [L_i,x^j]~=~\mathrm{i}\varepsilon_{ijk}~ x^k \tag{4}$$
$$ \Downarrow $$
$$[L^2,x^j]~=~\{L_i ,[L_i,x^j]\} ~\stackrel{(4)}{=}~\mathrm{i} \varepsilon_{ijk}~\{L_i ,x^k\} \tag{5}$$
$$ \Downarrow $$
$$[L^2,[L^2,x^j]]~\stackrel{(5)}{=}~
\mathrm{i} \varepsilon_{ijk}~\{L_i ,[L^2,x^k]\} 
~\stackrel{(5)}{=}~ 
\varepsilon_{ijk}\varepsilon_{k\ell n}\{L_i , \{L_{\ell} ,x^n\}  \}$$ 
$$~=~ \left(\delta_{i\ell}\delta_{jn} -\delta_{j\ell}\delta_{in}  \right)\{L_i , \{L_{\ell} ,x^n\}  \} ~=~ \{L_i , \{L_i ,x^j\}  \} -A_j, \tag{6}$$
where 
  $$A_j~:=~ \{L_i , \{L_j ,x^i\}  \} ~\stackrel{(3)}{=}~ [L_i , [L_j ,x^i] ] ~\stackrel{(4)}{=}~ \mathrm{i}\varepsilon_{ijk}[L_i ,x^k  ] ~\stackrel{(4)}{=}~  \varepsilon_{jik}\varepsilon_{ik\ell}~x^{\ell}~=~2x^j .\tag{7}$$
III) On the other hand, the RHS yields
$$ 2\{L^2, x^j\}~=~L_i \left( \{ L_i, x^j\} + [ L_i, x^j] \right)+\left( \{ x^j, L_i \} + [ x^j, L_i  ] \right) L_i ~=~ \{L_i , \{L_i ,x^j\}  \} +B_j,\tag{8} $$
where 
$$ B_j ~:=~ L_i[ L_i, x^j] +[ x^j, L_i  ]L_i
~\stackrel{(4)}{=}~\mathrm{i}\varepsilon_{ijk}[L_i,x^k] 
~\stackrel{(4)}{=}~-\varepsilon_{ijk} \varepsilon_{ik\ell}~x^{\ell}~=~-2x^j. \tag{9}$$
IV) Comparing the LHS and the RHS, we get OP's sought-for identity
$$ [L^2,[L^2,x^j]]~=~2\{L^2, x^j\} . \tag{10}$$
