I am trying to derive the form of the squared angular momentum using the index notation but I am stuck at the very last step.
First, I defined the orbital angular momentum using the index notation: $\textbf{L} = \textbf{r} \times \textbf{p}$, so that $L_k = \epsilon_{kmn} r_m p_n$, where $\epsilon_{kmn}$ is the Levi-Civita $\epsilon$.
$$ L^2 = \textbf{L} \cdot \textbf{L} = L_kL_k = \epsilon_{kmn} \ r_m \ p_n \times \epsilon_{k \bar m \bar n} \ r_{\bar m} \ p_{\bar n} = \epsilon_{kmn} \ \epsilon_{k \bar m \bar n} \ r_m \ p_n \ r_{\bar m} \ p_{\bar n} \tag{1} $$
At this stage, I used the identity: $$ \epsilon_{kmn} \ \epsilon_{k \bar m \bar n} = \delta_{m \bar m} \ \delta_{n \bar n} - \delta_{m \bar n} \ \delta_{\bar m n} \tag{2} $$ where $\delta$ is the Kronecker $\delta$. Substituting Eq. (2) into (1): $$ L^2 = (\delta_{m \bar m} \ \delta_{n \bar n} - \delta_{m \bar n} \ \delta_{\bar m n}) \ r_m \ p_n \ r_{\bar m} \ p_{\bar n} \tag{3} $$ $$ L^2 = \delta_{m \bar m} \ \delta_{n \bar n} \ r_m \ p_n \ r_{\bar m} \ p_{\bar n} - \delta_{m \bar n} \ \delta_{\bar m n} r_m \ p_n \ r_{\bar m} \ p_{\bar n} \tag{4} $$
The expression $\delta_{m \bar m} \ \delta_{n \bar n}$ in the first term of Eq. (4) implies that $m = \bar m$, and $n=\bar n$. Otherwise, $\delta=0$ and the term vanishes. We don't want this to happen. Similarly, in the second term, $m = \bar n$ and $\bar m = n$. Eq. (4) becomes: $$ L^2 = r_m \ p_n \ r_{m} \ p_{n} - r_m \ p_n \ r_{n} \ p_{m} \tag{5} $$ Attacking Eq. (5) by repeated application of the relation $rp-pr=i\hbar$ with the appropriate indices:
$$ L^2 = r_m (r_m p_n - i\hbar) p_n - r_m (r_n p_n - i\hbar) p_m \tag{6} $$ $$ L^2 = r_m r_m p_n p_n - i\hbar r_m p_n - r_m r_n p_n p_m + i\hbar r_m p_m \tag{7} $$ $$ L^2 = r^2 p^2 - i\hbar r_m p_n - r_m r_n p_n p_m + i\hbar (\textbf{r} \cdot \textbf{p}) \tag{8} $$ $$ L^2 = r^2 p^2 - i\hbar r_m p_n - r_n r_m p_n p_m + i\hbar (\textbf{r} \cdot \textbf{p}) \tag{9} $$ I have made use of the fact that $r_m$ and $r_n$ commute.
$$ L^2 = r^2 p^2 - i\hbar r_m p_n - r_n (p_n r_m + i\hbar) p_m + i\hbar (\textbf{r} \cdot \textbf{p}) \tag{10} $$ $$ L^2 = r^2 p^2 - i\hbar r_m p_n - r_n p_n r_m p_m - i\hbar r_n p_m + i\hbar (\textbf{r} \cdot \textbf{p}) \tag{11} $$ $$ L^2 = r^2 p^2 - i\hbar r_m p_n - (\textbf{r} \cdot \textbf{p})^2 - i\hbar r_n p_m + i\hbar (\textbf{r} \cdot \textbf{p}) \tag{12} $$ $$ L^2 = r^2 p^2 - (\textbf{r} \cdot \textbf{p})^2 + i\hbar (\textbf{r} \cdot \textbf{p}) - i\hbar (r_m p_n + r_n p_m) \tag{13} $$
The correct answer, however, is: $$ L^2 = r^2 p^2 - (\textbf{r} \cdot \textbf{p})^2 + i\hbar (\textbf{r} \cdot \textbf{p}) \tag{14} $$ where the last term, $- i\hbar (r_m p_n + r_n p_m)$, should vanish. I cannot see why it should vanish! Or maybe I did some mistake in the above steps. Any help?