Three spinless particles with two-body harmonic forces Problem
Consider a system of three spinless particles, each of mass $m$, moving on a straight line and bound to each other by two-body harmonic forces. The Hamiltonian of the system is
$$
H = − \frac{\hbar^2}{2m}
  \sum_{j = 1}^{3}
    \frac{\partial^2}{\partial x_j^2}
  +
  \frac{k}{2}
  \sum_{1 \leq j < k \leq 3} (x_j − x_k)^2
  \qquad
  (x_1, x_2, x_3 \in \mathbb{R},\;m, k > 0)
$$
acting on the Hilbert space $\mathcal{H} = L^2(\mathbb{R}^3, dx_1 dx_2 dx_3)$ and defined on its natural domain of self-adjointness.
Find a linear change of coordinates $R : (x_1, x_2, x_3) \mapsto (y_1, y_2, y_3)$ where, say, $y_3$ is the coordinate of the centre of mass, such that
– R induces a unitary transformation $U_R : L^2(\mathbb{R}^3, dx_1 dx_2 dx_3) \mapsto L^2(\mathbb{R}^3, dy_1 dy_2 dy_3)$,
– the transformed Hamiltonian $U_R H U_R^{−1}$ is the sum of three independent one-
body Hamiltonians.
Solution attempt
It should be noted that the term $\sum_{1 \leq j < k \leq 3} (x_j − x_k)^2$ is homogeneous, and can be written in the form $ \sum_{j, k=1}^3 Q_{jk} x_j x_k $, where we have the freedom to choose $Q$ to be a symmetric matrix.
The eigenvalues of Q are $3, 3, 0$, so $3$ has a multiplicity of 2, the system is underdetermined. Only one eigenvector is fixed, the other two could be arbitrarily chosen. Say we choose one eigenvector $\alpha := (\alpha_1, \alpha_2, \alpha_3)$ and select the remaining eigenvector with the Gram-Schmidt procedure.
So we have a parameterized family of orthogonal transformations $R_\alpha: (x_1, x_2, x_3) \mapsto (y_1, y_2, y_3)$ where now each $y_i$ depends on $\alpha$.
Now, we can acquire a family of Jacobi matrices $J_\alpha$ again depending on alpha, which appears in the Hamiltonian as $J_\alpha^T J_\alpha$ in the kinetic term, and I "hope", that the $\alpha$ could be chosen to yield $J_{\alpha_0}^T J_{\alpha_0} = \mathbb{1}$ for some value $\alpha_0$. Otherwise, other cross-terms will reappear in the kinetic term.
I didn't actually get to this point during my calculations yet, since the calculations got really-really messy (I miscalculated it several times), and I feel that there should be a much easier solution to this problem.
The problem says

Find a linear change of coordinates $R : (x_1, x_2, x_3) \mapsto (y_1, y_2, y_3)$ where, say, $y_3$ is the coordinate of the centre of mass...

This makes me think that there should be an easier way. In the case of the 2-particle problem, you could solve the problem easily via $(x_1, x_2) \mapsto (\frac{x_1 + x_2}{2}, x_1 - x_2)$.
I have read another possible solution here, but I think this solution is still too circumstantial.
Question: Is there a similar transformation for the 3-particle problem? Is there a set of rules or a rule of thumb that can be used to find transformations for problems like this?
Source:

*

*https://www.math.sissa.it/sites/default/files/Entrance_Examinations_pdf/Scritto2016.pdf
 A: I think I worked out the answer, and it turns out to be much simpler. I have not found this calculation anywhere explained in detail, so I hope it will help somebody.
TL;DR: An orthogonal transformation $R: (x_1, x_2, x_3) \mapsto (y_1, y_2, y_3)$ doesn't change a kinetic term with no crossterms and identical masses.
Explanation:
Let me rewrite the Hamiltonian in a more compact form:
$$
H = - \frac{\hbar^2}{2m}
\left ( \frac{\partial}{\partial \vec{x}} \right )^T I_{3 \times 3} \left ( \frac{\partial}{\partial \vec{x}} \right )
+ \frac{k}{2}
\vec{x}^T Q \vec{x},
$$
where $I_{3 \times 3}$ is just the identity matrix, and the potential term is
$$
\frac{k}{2} \sum_{1 \leq j < k \leq 3} (x_j - x_k)^2 =
\frac{k}{2} (x_1, x_2, x_3)^T
\begin{bmatrix} 
    2 & -1 & -1 \\
    -1 & 2 & -1 \\
    -1 & -1 & 2 \\
\end{bmatrix}
(x_1, x_2, x_3) =:  \frac{k}{2} \vec{x}^T Q \vec{x}.
$$
The matrix $Q$ can be diagonalized to give the eigenvalues $3, 0$ with $3$ with multiplicity 2. So there will be an orthogonal matrix $R$ which can be used to write
$$
Q = R D R^T,
$$
with $D$ being the diagonalized matrix, with the eigenvalues in the diagonal. We can always do this when $Q$ is a symmetric matrix.
Let us write $\vec{y} = R \vec{x}$. Then we know that $ \frac{\partial}{\partial \vec{y}} = R^{-1} \frac{\partial}{\partial \vec{x}}$ by using that (since it is an orthogonal matrix) and that the partial derivatives transform with the inverse of $R$, which is called the Jacobian matrix. Then we can easily substitute everything into the Hamiltonian as
$$
H = - \frac{\hbar^2}{2m}
\left (R^{-1} \frac{\partial}{\partial \vec{y}} \right )^T I_{3 \times 3} \left( R^{-1} \frac{\partial}{\partial \vec{y}} \right )
+ \frac{k}{2}
(R \vec{y})^T Q (R \vec{y})
\\=
- \frac{\hbar^2}{2m}
\left (\frac{\partial}{\partial \vec{y}} \right )^T R R^{-1} \left(  \frac{\partial}{\partial \vec{y}} \right )
+ \frac{k}{2}
\vec{y}^T R^T Q R \vec{y}
\\=
- \frac{\hbar^2}{2m}
\left (\frac{\partial}{\partial \vec{y}} \right )^T I_{3 \times 3} \left(  \frac{\partial}{\partial \vec{y}} \right )
+ \frac{k}{2}
\vec{y}^T D \vec{y},
$$
which is the desired, decoupled Hamiltonian.
An orthogonal transformation $R$ (there are many) can be found with the Gram-Schmidt procedure. One example is
$$
R = \frac{1}{\sqrt{6}}
\begin{bmatrix}
    \sqrt{3} & 1 & \sqrt{2} \\
    0 & -2 & \sqrt{2} \\
    \sqrt{3} & 1 & \sqrt{2} \\
\end{bmatrix},
$$
where the third column transforms to the center-of-mass $\frac{x_1 + x_2 + x_3}{3}$ with some scaling.
Most of the time you don't actually need to calculate the transformation itself, i.e. the eigenvalues of the Hamiltonian could be calculated only by knowing the eigenvalues of $Q$, and the kinetic term won't change due to the argument above.
