Coupled quantum harmonic oscillators (exact $\neq$ perturbative)

Question

Suppose we are given the Hamiltonian $$\hat H = \hat H_0 + \hat H_p(\varepsilon) = \frac 1 {2m}(\hat p_1^2 + \hat p_2^2) +\frac 1 2 m \omega^2(\hat x_1^2 + \hat x_2^2) + \varepsilon m\omega^2\hat x_1\hat x_2 $$ and, after switching to CM coordinates, $$ \begin{split} \hat X= \frac {\hat x_1 + \hat x_2}{2}, &\qquad \hat \rho = \hat x_1 -\hat x_2, \\ \hat P = \hat p_1 + \hat p_2, &\qquad \hat \pi= \frac {\hat p_1-\hat p_2}{2}, \end{split} $$ we are asked to solve the corresponding eigenvalue problem $\hat H \psi = E \psi$ first exactly and then perturbatively. I managed to find the following exact (normalized) solution: $$ \begin{split} \psi_{00}(X,\rho) &= \Gamma_{00} \exp\left\{-\frac{m\omega}{2\hbar}\left(X^2\sqrt{1+\varepsilon} + \rho^2\sqrt{1-\varepsilon} \right) \right\}, \\ \Gamma_{00} &= \sqrt{\frac{m\omega}{\pi\hbar}\sqrt[4]{(1+\varepsilon)(1-\varepsilon)}}, \\ \psi_{n\ell}(X,\rho) &= \frac{1}{\sqrt{2^{n+\ell}n!\ell!}}\operatorname{He}_n\left(X\sqrt{\frac{2m\omega}{\hbar}\sqrt{1+\varepsilon}} \right) \operatorname{He}_\ell\left(\rho \sqrt{\frac{m\omega}{2\hbar}\sqrt{1-\varepsilon}} \right) \psi_{00}(X,\rho), \\ E_{n\ell} &= \hbar \omega \left\{\left(n+\frac 1 2 \right) \sqrt{1+\varepsilon} + \left(\ell + \frac 12 \right)\sqrt{1-\varepsilon}\right\}, \end{split} $$ where $\operatorname{He}_q(\xi)$ is the usual $q$-th Hermite polynomial. Treating the problem perturbatively, again in the new coordinates, the analysis leads to the first-order energy corrections $$\delta^{(1)}E_{00} = 0, \qquad \delta^{(1)}E_{01} = -\frac 1 2 \varepsilon\hbar \omega, \qquad \delta^{(1)}E_{10} = \frac12 \varepsilon \hbar \omega, $$ which are perfectly compatible with the first-order expansion of the exact energies. However, the first-order correction to the unperturbed eigenstate $\psi_{00}^{(0)}$ (I haven't tackled the others) does not match what I get when expanding $\psi_{00}$ above to first order in $\varepsilon$. Why is this?

---

Calculations. The theory of stationary perturbations gives $$\delta^{(1)}\psi_{00} = \sum_{n,\ell\neq 0} \frac{\langle \psi_{n\ell}^{(0)}| \hat H_p \psi_{00}^{(0)} \rangle}{E_{00}^{(0)}-E_{n\ell}^{(0)}} \psi_{n\ell}^{(0)}; $$ the matrix element in the numerator can be calculated through appropriate ladder operators $\hat A$ and $\hat \alpha$ (respectively for the unperturbed oscillator in $\hat X$ and $\hat \rho$) to be, knowing $\hat H_p = \varepsilon m \omega^2(\hat X^2 - \frac 1 4 \hat \rho^2)$ and remembering $(\hat a^\dagger + \hat a)^2 = (\hat a^\dagger)^2 + \hat a^2 + 2\hat a^\dagger \hat a + \hat 1$ when $[\hat a,\hat a^\dagger] = \hat 1$, $$\begin{split} \langle \psi_{n\ell}^{(0)}| \hat H_p \psi_{00}^{(0)} \rangle &= \varepsilon m \omega^2 \Bigg\{ \frac{\hbar}{4m\omega} \left( \langle \psi_{n\ell}^{(0)}| (\hat A^\dagger)^2 \psi_{00}^{(0)} \rangle + \underbrace{\langle \psi_{n\ell}^{(0)}| (\hat A)^2 \psi_{00}^{(0)} \rangle}_{=0} + 2\underbrace{\langle \psi_{n\ell}^{(0)}| \hat A^\dagger \hat A \psi_{00}^{(0)}\rangle}_{=0} + \underbrace{\langle \psi_{n\ell}^{(0)}| \psi_{00}^{(0)} \rangle}_{=0} \right) \\ &\qquad\qquad - \frac{\hbar}{4m\omega} \left( \langle \psi_{n\ell}^{(0)}| (\hat \alpha^\dagger)^2 \psi_{00}^{(0)} \rangle + \underbrace{\langle \psi_{n\ell}^{(0)}| (\hat \alpha)^2 \psi_{00}^{(0)} \rangle}_{=0} + 2\underbrace{\langle \psi_{n\ell}^{(0)}| \hat \alpha^\dagger \hat \alpha \psi_{00}^{(0)}\rangle}_{=0} + \underbrace{\langle \psi_{n\ell}^{(0)}| \psi_{00}^{(0)} \rangle}_{=0} \right) \Bigg\} \\ &= \frac{1}{4}\varepsilon \hbar \omega \left\{\sqrt 2 \langle \psi_{n\ell}^{(0)}| \psi_{20}^{(0)} \rangle - \sqrt 2 \langle \psi_{n\ell}^{(0)}| \psi_{02}^{(0)} \rangle \right\}, \end{split} $$ so that $$\delta^{(1)}\psi_{00} = \frac{\varepsilon \hbar \omega}{2\sqrt 2} \left[ \frac{1}{-\hbar\omega(2+0)}\psi_{20}^{(0)} - \frac{1}{-\hbar\omega(0+2)} \psi_{02}^{(0)} \right] = -\frac 1 8 \varepsilon \sqrt 2 (\psi_{20}^{(0)} - \psi_{02}^{(0)}).$$ However, $$ \begin{split} \psi_{00}(X,\rho) &= \left( \sqrt{\frac{m\omega}{\pi\hbar}} + \mathscr O(\varepsilon^2) \right) \exp \left\{ - \frac{m\omega}{2\hbar} \left[X^2\left(1+\frac 1 2 \varepsilon\right) + \rho^2\left(1-\frac 1 2 \varepsilon\right) + \mathscr O(\varepsilon^2) \right] \right\} \\ &= \sqrt{\frac{m\omega}{\pi\hbar}} \exp \left\{ - \frac{m\omega}{2\hbar} \left(X^2 + \rho^2 \right) \right\} \exp \left\{ - \frac{m\omega}{4\hbar} \varepsilon (X^2-\rho^2) \right\} + \mathscr O(\varepsilon^2) \\ &= \psi_{00}^{(0)}(X,\rho) \left\{1 - \frac{m\omega}{4\hbar} \varepsilon (X^2-\rho^2) + \mathscr O(\varepsilon^2) \right\} + \mathscr O(\varepsilon^2); \end{split}$$ hence the correction must be a linear combination of the unperturbed states $\alpha_{20}\psi_{20}^{(0)}+ \alpha_{02}\psi_{02}^{(0)}$ (which agrees with the perturbative result), in such a way that (substituting those states explicitly) $$ - \frac{m\omega}{4\hbar} \varepsilon (X^2-\rho^2) = \alpha_{20} X^2 \frac{m\omega}{\hbar} \sqrt{8} - \frac{1}{\sqrt 2}\alpha_{20} + \alpha_{02} \rho^2 \frac{m\omega}{\hbar} \frac{1}{\sqrt 2} - \frac{1}{\sqrt 2}\alpha_{02}. $$ But comparing coefficients leads to an (overdetermined) system of equations for the $\alpha$s, giving $$\alpha_{20}= - \frac{1}{8\sqrt 2}\varepsilon, \qquad \alpha_{02} = \frac{1}{2\sqrt 2} \varepsilon, $$ which disagrees with the perturbative prediction (apart from being inconsistent with the requirement $\alpha_{20} = -\alpha_{02}$ coming from the constant terms above). Where have I gone wrong?

Answer 1

The problem is that you haven't used a consistent normalization of $X$ and $\rho$. To see this explicitly, set $\epsilon = 0$. Then by symmetry, the ground state should only depend on the combination $x_1^2 + x_2^2$. But your ground state depends only on $X^2 + \rho^2$. This would have been correct, had you defined $$X = \frac{x_1 + x_2}{2}, \quad \rho = \frac{x_1 - x_2}{2}$$ but you instead defined one of them with a $1/2$ and one of them without it. This missing factor of $2$ squares into a factor of $4$, which is exactly the factor of $4$ you're missing at the end.

As a aside, it's generally a good idea to set $m = \omega = \hbar = 1$ at the beginning of calculations like these. The simpler form makes it much easier to spot issues, and you can always restore the factors with dimensional analysis.