For 3-particle systems the relevant permutation group is $S_3$, which is significantly different from $S_2$ in that not all representations of $S_3$ are 1-dimensional.
Whereas 2-particle states can always be written so they transform back to themselves under permutation of labels, i.e. always possible to write
$$
\psi_a(x_1)\psi_b(x_2)=\frac{1}{2}\psi^+_{ab}(x_1,x_2) + \frac{1}{2}\psi^-_{ab}(x_1,x_2)
$$
with
$$
P_{12}\psi^{\pm}_{ab}(x_1,x_2)=\pm \psi^{\pm}_{ab}(x_1,x_2)
$$
the symmetric and antisymmetric combinations, there are 3-particle states that cannot be made to transform back to a multiple of themselves. Already this is apparent for the triple coupling of spin-1/2 particles: there are two copies of total spin $1/2$ in the decomposition and you cannot arrange the states in $S=1/2$ to be symmetric or antisymmetric under all permutations of particles.
The situation does not improve for 3 spin-1 particles. The fully symmetric states will have total $L=3,1$, the fully antisymmetric states will have $L=0$, but there are states of mixed symmetry with $L=1,2$ (there are two copies of these mixed symmetry states).
Thus, you have no obvious way of knowing if a state with $L=1$ is symmetric or of mixed symmetry.
Of course the spin-statistics theorem does not forbid spin-states of mixed symmetry: the must be combined with spatial stages (also of mixed symmetry) to produce a properly symmetrized total state.
You have mixed indices in the end of your line. Correctly:
$$-\epsilon_{iab}\epsilon_{jcd}(x_ap_d\delta _{bc} - x_cp_b\delta _{ad}).$$
So further, $-\epsilon_{iab}\epsilon_{jcd}(x_ap_d\delta _{bc} - x_cp_b\delta _{ad})=\epsilon_{iab}\epsilon_{jcd}x_cp_b\delta _{ad}-\epsilon_{iab}\epsilon_{jcd}x_ap_d\delta _{bc}=\\=\epsilon_{idb}\epsilon_{jcd}x_cp_b - \epsilon_{iac}\epsilon_{jcd}x_ap_d=\epsilon_{iac}\epsilon_{jdc}x_ap_d-\epsilon_{ibd}\epsilon_{jcd}x_cp_b=\\=(\delta _{ij}\delta _{ad}-\delta _{id}\delta _{aj})x_ap_d-(\delta _{ij}\delta _{bc}-\delta _{ic}\delta _{bj})x_cp_b=\\=\delta _{ij}x_ap_a-x_jp_i-\delta _{ij}x_bp_b+x_ip_j=x_ip_j-x_jp_i=\epsilon_{ijk}L_k.$
Best Answer
There are two things going on in the last step. Let's take care of the simpler one first. The third term reduces thanks to a trigonometric identity: $$ -\hbar^2 \left ( 1 + \frac{\cos^2 \theta}{\sin^2 \theta} \right ) \frac{\partial^2}{\partial \varphi^2} = -\hbar^2 \frac{\sin^2 \theta + \cos^2 \theta}{\sin^2 \theta} \frac{\partial^2}{\partial \varphi^2} = -\hbar^2 \frac{1}{\sin^2 \theta} \frac{\partial^2}{\partial \varphi^2} $$
The first two terms are more interesting, and simplifying them involves recognizing something really cool: operator notation. $L^2$ is, as you know, an operator, not a number or a function, because it takes in functions and returns functions. If a mathematician were to classify it, they might say $L^2 \in \mathrm{Map}(\mathrm{Map}(\mathbb C,\mathbb C),\mathrm{Map}(\mathbb C,\mathbb C)) $. (Of course a mathematician could classifiy it more strictly, but that's not important now.) Now, in a lot of ways the rules for algebra with operators are the same as for numbers or functions, but there is one important caveat. You can't deal with other operators, particularly differential operators, in the same way. To get things right, it helps immensely to write the expression as if it were applied to a function. For example, the first two terms, $$ -\hbar^2 \frac{\partial^2}{\partial \theta^2} -\hbar^2 \frac{\cos\theta}{\sin\theta}\frac{\partial}{\partial \theta} $$ become $$ -\hbar^2 \frac{\partial^2 \Psi}{\partial \theta^2} -\hbar^2 \frac{\cos\theta}{\sin\theta}\frac{\partial \Psi}{\partial \theta} $$ Now this expression describes a function, not an operator, and for people used to normal calculus notation, this is much easier to deal with without making mistakes. Let's factor some things out. $$ -\hbar^2 \frac{\partial^2 \Psi}{\partial \theta^2} -\hbar^2 \frac{\cos\theta}{\sin\theta}\frac{\partial \Psi}{\partial \theta} = -\hbar^2 \frac 1{\sin\theta} \left( \sin \theta \frac{\partial^2 \Psi}{\partial \theta^2} + \cos\theta \frac{\partial \Psi}{\partial \theta} \right ) $$ Do you see it yet? The last step is just the product rule! $$ -\hbar^2 \frac 1{\sin\theta} \left( \sin \theta \frac{\partial^2 \Psi}{\partial \theta^2} + \cos\theta \frac{\partial \Psi}{\partial \theta} \right ) = -\hbar^2 \frac 1{\sin\theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial \Psi}{\partial \theta} \right ) $$ Then, to get back to operator notation, one simply doesn't write $\Psi$, because it was a generic function all along.
Indulge me for one more simple example, so that I may illustrate how subtle the difference between operators and functions can be. Take two operators, $A = \frac{\partial}{\partial \theta}$ and $B = \sin \theta$. What is their composite, $AB$? The naïve answer is "$\cos \theta$", which is incorrect. To arrive at the correct answer, we must recognize that $B$ is not a function, it is an operator. It is always implied that $B \Psi = \sin \theta \Psi$, even though no operator appears to be written. The correct answer to my question follows from the product rule: $$ AB = \cos \theta + \sin \theta \frac{\partial}{\partial \theta} $$ Go ahead and verify it for yourself by writing it with a test function instead of as an operator.