I'm working with two neutrons bound in a 3D harmonic trap (so, living in a harmonic oscillator potential) which also interact with each other with a given potential. The states of the individual neutrons can be written as $|n_1l_1 s_1 \rangle$, $|n_2 l_2 s_2 \rangle$ (I have ommited the $m$ quantum number), and the product state as $|n_1 l_1 s_1, n_2 l_2 s_2\rangle$, with the convenient antisymmetrization.
One of the terms of the interaction potential between the neutrons is of the form $V(r)\frac{1}{2}(1 + P_r)$, where $P_r$ is an operator that exchanges the spatial degrees of freedom of the two particles ($P_r|n_1 l_1 s_1, n_2 l_2 s_2\rangle =|n_2 l_2 s_1, n_1 l_1 s_2\rangle$). So, $\frac{1}{2}(1 + P_r)$ projects onto the "even" part of the nuclear two body wavefunction (part that does not change sign on the exchange of the spatial degrees of freedom).
My question is: How do I apply this projector in the calculation of a matrix element of the form $\langle n_1l_1s_1, n_2l_2s_2 |V(r)P_r| n_3l_3s_3, n_4l_4s_4\rangle$. In other words, how is the state $|n_2 l_2 s_1, n_1 l_1 s_2\rangle$ related to $|n_1 l_1 s_1, n_2 l_2 s_2\rangle$ in general?
