Quantum harmonic oscillators with momentum-position coupling I have two coupled quantum harmonic oscillators given by the following Hamiltonian:
$$H=\frac{p_{x}^{2}}{2}+\frac{\omega^{2} x^{2}}{2}+\frac{p_{y}^{2}}{2}+\frac{\Omega^{2} y^{2}}{2}+\frac{C p_{x} y}{2}.$$
As you can see, the coupling is done over the position variable of one oscillator and the momentum of the other. I need to find the wavefunction of states where either of the two oscillators (or both) is excited.
What I tried to do:
In general when I have couple harmonic oscillators where the coupling term is of the form $C (x_1-x_2)^2$ I start by diagonalising the Hamiltonian, then define normal coordinates and write down the Hamiltonian using those. I.e., I write the Hamiltonian under the form
$$H=\frac{1}{2} \sum_{i=1}^{2} p_{i}^{2}+\frac{1}{2} \sum_{i, j=1}^{2} x_{i} K_{i j} x_{j}.$$
Then I diagonalise $K$,
$$K_D = UKU^T.$$
And define normal coordinates
$$\left(\begin{array}{l}
x_{+} \\
x_{-}
\end{array}\right) \equiv U\left(\begin{array}{l}
x_{1} \\
x_{2}
\end{array}\right),$$
leading to
$$H=\frac{1}{2}\left[p_{+}^{2}+p_{-}^{2}+\omega_{+}^{2} x_{+}^{2}+\omega_{-}^{2} x_{-}^{2}\right].$$
I am unable to use this here since I do not know how to diagonalise this Hamiltonian with the term $C p_x y$ present.
 A: Your problem is basically  that of a charged particle in a magnetic field in the Landau gauge, but
I'll explain the general theory for such systems as it's both interesting and suprisingly complicated --- and seldom found in textbooks. It's essentially  the theory of Bose Bogoluibov transformations.
Suppose we are given classical (or quantum) quadratic  Hamiltonian
$$
{\mathcal H}[{\bf p},{\bf x}]= \frac 12 M_{ij} p_ip_j+ \frac 12 V_{ij} x_ix_j +  K_{ij}\,p_i x_j \nonumber\\
=\frac 12  \left[\matrix{{\bf p}^T& {\bf x}^T}\right] H \left[\matrix{ {\bf p}\cr {\bf x}}\right], \nonumber
$$
where  and  $H$ is the real, positive definite, symmetric $2N$-by-$2N$ matrix
$$
H= \left[\matrix{M& K\cr K^T &V}\right].
$$
(If $H$ is not positive definite the diagonalization is not possible   as the energy of the system is not bounded below and the eigenfrequencies are pure imaginary.)
We seek a  transformation
$$
 \left[\matrix{ {\bf p}\cr {\bf x}}\right]= S \left[\matrix{ {\bf P}\cr {\bf X}}\right]
$$
that diagonalizes $H$ via  the congruence transformation ${H}\to S^T HS $ and preserves the classical Poisson brackets and/or the quantum commutation relations. This preservation   requires
$$
S^T JS=J, \quad J= \left[\matrix{0&- {\mathbb I}_n \cr {\mathbb I}_n&\phantom{-} 0}\right],
$$
so $S\in {\rm Sp}[2N, {\mathbb R}]$, the non-compact symplectic group of linear canonical transformations.   That positive definiteness of $H$ is sufficient to ensure that we can find such an $S$ is the statement  of Williamson's theorem. The resulting frequencies are called the symplectic spectrum of $H$. They are unrelated to the eigenvalues of $H$. Further  matrix $S$ is neither orthogonal nor unitary, and so does not obey $S^{-1}=S^\dagger$
We observe  that $H$ being positive definite ensures that the matrices $H^{\pm 1/2}$ are well defined. We then see that
the matrix
$$
\tilde J= H^{-1/2}JH^{-1/2}
$$
is skew symmetric and hence an element of $\mathfrak{so}[2N]$. Therefore  there exists an     $R\in {\rm O}(2N)$  that conjugates $\tilde J$ into the Cartan algebra  of $\mathfrak{so}(2n)$ --- i.e.
$$
R^TH^{-1/2} J H^{-1/2} R= \left[\matrix{0&- \Omega^{-1} \cr \Omega^{-1} &0}\right],
$$
with  $\Omega^{-1}$ diagonal. As we are allowing  $R\in {\rm O}[2N]$ rather than demanding  $  R\in {\rm SO}[2N]$,  we can ensure that   $\Omega^{-1}$   possesses  strictly positive entries. This being so, we define
$$
 D= \left[\matrix{\Omega^{1/2}&0 \cr 0& \Omega^{1/2}}\right],
 $$
so that
$$
 DR^T H^{-1/2} J H^{-1/2} RD=J,\nonumber\\
 DR^T H^{-1/2}H H^{-1/2} R D= D^2.\nonumber
$$
Thus
$$
 S= H^{-1/2} RD \in {\rm Sp}[2N, {\mathbb R}],
 $$
and
$$
 S^T HS= D^2 = \left[\matrix{\Omega&0 \cr 0& \Omega}\right]
 $$
where  $\Omega= {\rm diag}(\omega_1,\ldots \omega_N)$.  We have
$$
{\mathcal H}[{\bf P}, {\bf Q}]=\frac 12  \left[\matrix{ {\bf P}^T& {\bf X}^T}\right]\left[\matrix{ \Omega &0\cr 0&\Omega}\right]\left[\matrix{ {\bf P}\cr {\bf X}}\right]=
\sum_i \frac {\omega_i}{2} (P_i^2+X_i^2),
$$
which is a set of
decoupled harmonic oscillators with frequency $\omega_i$.
We can of course scale $X_i\to \sqrt{\omega} X_i$ and $P_i\to P_i/\sqrt{\omega}$ without changing the commutation relations, and so
$$
\sum_i \frac {\omega_i}{2} (P_i^2+X_i^2)\to \sum_i \frac {1}{2} (P_i^2+\omega_i^2X_i^2),
$$
which probably looks more familiar.
The eigenvalues $\omega_i$ are most easily found from those of   $JH$ which are  $\pm i\omega_i$.
Note added: I see that this same strategy is discussed in a math stack exchange answer: https://math.stackexchange.com/questions/1171842/finding-the-symplectic-matrix-in-williamsons-theorem
A: I am liking the answers which have been posted, but figure to offer an answer which uses more primitive means.
This system has a classical analog with Lagrangian of the form:
$$
L=\frac{1}{2}\cdot \dot{x}\cdot M\cdot \dot{x}-\frac{1}{2}\cdot x\cdot K\cdot x\mp \frac{1}{2}\cdot \dot{x}\cdot C\cdot x
$$
I say $\mp \frac{1}{2}\cdot \dot{x}\cdot C\cdot x$ since I am not 100 % certain that the $p_x y$ term is to be considered momentum - potential energy as is the case in electrodynamics.
Lagrange' s equation $\frac{d}{dt}\cdot \frac{\partial L}{\partial \dot{x}}=\frac{\partial L}{\partial x}$ produces:
$$
\left( M^T+M \right) \cdot \ddot{x}  \mp C\cdot \dot{x}=\mp C^T\cdot \dot{x}-\left(K^T+K\right)\cdot x
$$
Which produces an eigenvalue equation and I assume that all who are reading this know how to find the normal modes of this system such that the Lagrangian becomes:
$$
L=\frac{1}{2} \underset{\alpha }{\Sigma } \left(\dot{Q}_{\alpha }^2-\omega _{\alpha }^2 Q_{\alpha }^2\right)
$$
where the $\left\{\overset{\to }{Q}_{\alpha }\right\}$ and $\left\{\omega _{\alpha }\right\}$ are the eigenmodes and eigenfrequencies.
And so the classical Hamiltonian is:
$$
H=\frac{1}{2} \underset{\alpha }{\Sigma } \left(p_{\alpha }^2+\omega _{\alpha }^2 Q_{\alpha }^2\right)
$$
And the only step left is to use canonical quantization :
$$
Q_{\alpha }\to \hat{Q}_{\alpha },p_{\alpha }\to \hat{p}_{\alpha }=-i\hbar \cdot \frac{\partial }{\partial Q_{\alpha }}
$$
where I use the hats to indicate that $\hat{Q}$ and $\hat{p}$ are quantum operators.
A: I started off by assuming a general transformation of the form:
$$\left\{\begin{array}{l}
p_{x}=v p_{-}+c p_{+}+d x_{+}, \\
p_{y}=k p_{-}+l p_{+}, \\
x=n x_{+}+q x_{-}, \\
y=u x_{+}+ax_{-}+b p_{+} ,
\end{array}\right.$$
with ${p_+, p_-,x_+, x_-}$ the new variables. To avoid confusion, I am also renaming the coupling coefficient $C$ to $\gamma$. Next I plugged the above transformations in the Hamiltonian and demanded that no off-diagonal terms survive. This constraints the transformation coefficients. I also demand that the coefficients next to $p_+^2 $ and $p_-^2$ in the Hamiltonian be $1$. This leads to the set of transformations
$$\left\{\begin{array}{l}
p_{x}=-\sqrt{1-l^{2}} p_{+}+l p_{m}-\frac{a \gamma}{2} x_{-}, \\
p_{y}=k p_{+}+\sqrt{1-l^{2}} p_{-}, \\
x=n x_{+}, \\
y=a x_{-} .
\end{array}\right.$$
We are now left with the set of coefficients ${a, l, n}$ to be determined. This can be done by enforcing the canonical commutation relations on the new and old sets of coordinates. We then find
$$[x, p_y]=0 \Rightarrow l=0, \\
[x, p_x]=i \hbar \Rightarrow n=-1, \\
[y, p_y] = i \hbar \Rightarrow a=1.$$
The resulting set of transformations is hence
$$\left\{\begin{array}{l}
p_{x}=-p_{+}- \frac{\gamma}{ 2} x . \\
p_{y}=p_{-}, \\
x=-x_{+}, \\
y=x .
\end{array}\right.$$
leading to the Hamiltonian
$$H = p_{+}^{2}+p_{-}^{2}+\omega^{2} x_{+}^{2}+\left(\Omega^{2}-\gamma^{2} / 4\right) x_{-}^{2}.$$
The normal modes are obviously
$$\omega_+ = \pm \omega, \,\,\,\,\,\,\,\, \omega_- = \sqrt{\left(\Omega^{2}-\gamma^{2} / 4\right)}.$$
