# 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),$$

$$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.

• Since you have posted a bounty asking for a more detailed answer, could you clarify what about Mike Stone's answer is insufficient? It looks rather comprehensive to me. Jun 27, 2021 at 4:43
• @RichardMyers MS’s answer is great but there are a few details in it that I do not understand how to implement exactly. For instance how to I get $R$? How exactly do I compute $H^{-1/2}$?
– Y2H
Jun 27, 2021 at 7:22
• @RichardMyers another issue is that, when I define $D$ as Mike Stones notes and assume no constraints on $S$ at all, Mathematica simply shows that the equation $S^T H S$ cannot be satisfied.
– Y2H
Jun 27, 2021 at 12:10

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

• Never seen such computations. Very cool! +1 Jun 13, 2021 at 13:43
• I am unsure If one could find a symmetric K matrix such that the proposed form of the Hamiltonian matches the question's Hamiltonian since there is only a $p_x y$ term. Also, is it 100% sure that the $p_k$ in the first equation should not be a $p_j$ ? Aside from that, thank you very much for posting this since it looks like it works beautifully when there are xipj and xjpi terms. Jun 13, 2021 at 17:31
• Thanks for catching typo (now fixed)! I think that you have $H=(1/2)(p_x+\lambda y)^2- (1/2)\lambda^2 y^2 +\ldots$ which is why I was suggesting that you regard $\lambda y=A_x$ as a gauge field. YOu can also add a $xp_y-yp_x$ angular momentum term to make symmetric but may nt help much. Jun 13, 2021 at 18:19
• I don't know of one. I chased this stuff down in papers. In particular Simon, Chaturvedi, Srinivasan "Congruences and Canonical Forms for a Positive Matrix: Application to the Schweinler-Wigner Extremum Principle" arXiv:math-ph/9811003 Jun 14, 2021 at 11:49
• Juts look at my notes oin this: people.physics.illinois.edu/stone/bose_bogoliubov.pdf Jun 26, 2021 at 13:10

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.

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.$$

$$H = p_{+}^{2}+p_{-}^{2}+\omega^{2} x_{+}^{2}+\left(\Omega^{2}-\gamma^{2} / 4\right) x_{-}^{2}.$$
$$\omega_+ = \pm \omega, \,\,\,\,\,\,\,\, \omega_- = \sqrt{\left(\Omega^{2}-\gamma^{2} / 4\right)}.$$