Consider a two-dimensional harmonic oscillator with the potential \[\begin{equation} V (x, y) = \frac{1}{2} m\omega^2(x^2 + y^2)\tag{1} \end{equation}\] We want to consider a degenerate enegry state with the eigenvalue \(3 \hbar \omega\), and compute the first-order corrections to the above energy level due to the perturbation \(H_1 = K x y\).
The level with energy \(3 \hbar w\) is triple degenerate with the states \(|20\rangle\), \(|11\rangle\) and \(|02\rangle\). We need to find the matrix elements of \(H_1=k x y\) in this degenerate subspace basis. The matrix we have will be \(3\times 3\) with \(9\) entries. But not all of them are nonzero (since \(x\) and \(y\) can change the state by \(1\) unit only), furthermore as \(xy\) is an Hermitian and real operator, the matrix will be symmetric (with respect to usual diagonal line). And finally the symmetry between \(x\) and \(y\) will force the matrix to be symmetric with respect to the other diagonal line. To sum up we need to calculate only \(1\) entry, others will be dictated by the symmetry, and will be the same. Let’s calculate \[\begin{equation} \langle 11|k x y|20\rangle=k\langle 1 |x|0\rangle\langle 1|y|2\rangle=\frac{k}{\sqrt 2 m w}. \tag{2} \end{equation}\] The matrix we have is, \[\begin{equation} H_1=\frac{k}{m w}\left( \begin{array}{ccc} 0 & 1/\sqrt 2 & 0 \\ 1/\sqrt 2 & 0 & 1/\sqrt 2 \\ 0 & 1/\sqrt 2 & 0 \\ \end{array} \right)=\frac{k}{m w} J_x. \tag{3} \end{equation}\] The eigenvalues are easy to get since \(J_x\) is in the spin-\(1\) space, so it will have eigenvalues \(0,1,-1\). For our matrix the eigenvalues are \(0,\pm \frac{k}{m w}\), which are the first order corrections to the energy, \(3w\). Now if we want to find the states corresponding to these energy levels we just find the eigenvectors of the above matrix, which are also the eigenvectors of \(J_x\). Using the eigenvectors, we have the following first order corrected eigenstates of \(H\), \[\begin{equation} |0\rangle=\frac{|20\rangle-|02\rangle}{\sqrt 2},|-\rangle=\frac{|20\rangle-\sqrt 2 |11\rangle+|02\rangle}{ 2}, |+\rangle=\frac{|20\rangle+\sqrt 2 |11\rangle+|02\rangle}{ 2} \tag{4} \end{equation}\] with the energies \(3w\), \(3w-\frac{k}{m w}\) and \(3w+\frac{k}{m w}\), respectively. It is important to note that this problem could have been solved exactly by diagonalizing the potential \(V=\frac{1}{2} m w^2(x^2+y^2)+k x y\) which can be done transforming the coordinates \(x,y\) to the coordinates \(\tilde{x}, \tilde{y}\). The transformation is found by diagonalizing the matrix,
\[\begin{equation} V_{i j}=\frac{\partial^2}{\partial x^i \partial x^j} V=\left(\begin{array}{cc} m w^2 & k \\ k & m w^2 \end{array}\right) \tag{5} \end{equation}\] The eigenvalues are \(m w^2+k\) and \(m w^2-k\) with the eigenvectors \(\frac{1}{\sqrt 2}\left(\begin{array}{cc} 1 \\ 1 \end{array}\right)\) and \(\frac{1}{\sqrt 2}\left(\begin{array}{cc} 1 \\ -1 \end{array}\right).\) What this algebra tells us is that if we define new coordinates \(\tilde{x_1}=\frac{x+y}{\sqrt 2}\) and \(\tilde{x_2}=\frac{x-y}{\sqrt 2}\), then \[\begin{eqnarray} % \nonumber to remove numbering (before each equation) V(x,y) &=& \frac{1}{2}V_{i j} x_i x_j= \frac{1}{2} m w^2(x^2+y^2)+k x y\nonumber\\ &=&\frac{1}{2}\tilde{V}_{i j} \tilde{x}^i \tilde{x}^j=\frac{1}{2}( m w^2+k) \tilde{x_1}^2+\frac{1}{2}( m w^2-k) \tilde{x_2}^2 \tag{6} \end{eqnarray}\] You can check that the momentum part is transformed simply by renaming \(p_i\) by \(\tilde{p_i}\). So the final Hamiltonian we have, is simply the Hamiltonian for two decoupled harmonic oscillators with different frequencies. The the total energy of the system is given by \(E=\sqrt{w^2+\frac{k}{m}}(\tilde{n_1}+\frac{1}{2})+\sqrt{w^2-\frac{k}{m}}(\tilde{n_2}+\frac{1}{2})\). Now as this is the exact answer, it must be equal to the approximate answer found before, at the first order in \(k\). To compare we need to set the values \(\tilde{n_1}\) and \(\tilde{n_2}\) to the values which produce \(E=3 w\) when \(k=0\). The results are \(|\tilde{2}\tilde{0}\rangle\), \(|\tilde{1}\tilde{1}\rangle\) and \(|\tilde{0}\tilde{2}\rangle\), but one must be careful that the states are not identical to the previous ones as the new pairs are in the new coordinates (that is why they have \(\;\;\tilde{}\;\;\)). The exact energy levels can be expanded at the first order as \[\begin{eqnarray} % \nonumber to remove numbering (before each equation) E_{\tilde{1}\tilde{1}} &=& \left[\sqrt{w^2+\frac{k}{m}}(1+\frac{1}{2})+\sqrt{w^2-\frac{k}{m}}(1+\frac{1}{2})\right]= 3 w+ O(k^2)\nonumber\\ E_{\tilde{2}\tilde{0}} &=& \left[\sqrt{w^2+\frac{k}{m}}(2+\frac{1}{2})+\sqrt{w^2-\frac{k}{m}}(\frac{1}{2})\right]= 3 w+ \frac{k}{mw}+ O(k^2)\nonumber\\ E_{\tilde{0}\tilde{2}} &=& \left[\sqrt{w^2+\frac{k}{m}}(\frac{1}{2})+\sqrt{w^2-\frac{k}{m}}(2+\frac{1}{2})\right]= 3 w-\frac{k}{mw}+ O(k^2))\tag{7} \end{eqnarray}\] We note that the agreement of the exact results for energy with the perturbation theory results at the first order.
PS: The image in the thumbnail is taken from [1].