I want to find the wave functions of $N$ coupled quantum harmonic oscillators having the following hamiltonian:
\begin{eqnarray}
H &=& \sum_{i=1}^N \left(\frac{p^2_i}{2m_i} + \frac{1}{2}m_i\omega^2 x^2_i + \frac{\kappa}{2} (x_i-x_{i+1})^2 \right)\,, \qquad x_{N+1}=0\,,\\
&=& \frac{1}{2}p^T Mp + \frac{1}{2}x^TKx\,,
\end{eqnarray}
where $M=\text{diag}(\frac{1}{m_1}, \cdots,\frac{1}{m_N})$ and $K$ is a real symmetric $N\times N$ matrix with positive eigenvalues,
\begin{equation}
K=
\begin{pmatrix}
k'_1& -\kappa & 0 & \cdots & 0 \\
-\kappa & k'_2& -\kappa & \ddots & \vdots \\
0 & -\kappa & \ddots& \ddots & 0\\
\vdots & \ddots &\ddots & k'_{N-1}&-\kappa \\
0&\cdots & 0 & -\kappa & k'_N
\end{pmatrix}
\end{equation}
with $k'_i = m_i\omega^2+2\kappa$ but $k'_{1,N} = m_{1,N}\omega^2+\kappa$. By choosing a basis which diagonalizes the matrix $K$, the hamiltonian can be express as the sum of uncoupled harmonic oscillators hamiltonian.
As an example, consider two coupled quantum harmonic oscillators with hamiltonian
\begin{equation}
H = \frac{p^2_1}{2m_1} + \frac{p^2_2}{2m_2} + \frac{1}{2}m_1\omega^2 x^2_1 + \frac{1}{2}m_2 \omega^2 x^2_2 + \frac{\kappa}{2} (x_1-x_2)^2 \,.
\end{equation}
We make the following changes of variables (normal coordinates)
\begin{eqnarray}
x &=& \frac{x_1 – x_2}{\sqrt{2}} \,, \\
X &=& \frac{m_1 x_1 + m_2 x_2}{M\sqrt{2}}\,,
\end{eqnarray}
or equivalently,
\begin{eqnarray}
x_1 &=& \frac{1}{\sqrt{2}}\left(X + \frac{m_2}{M}x\right) \,, \\
x_2 &=& \frac{1}{\sqrt{2}}\left(X – \frac{m_1}{M}x\right) \,,
\end{eqnarray}
where $M=(m_1+m_2)/2$. Then the hamiltonian becomes
\begin{equation}
H = \frac{p^2_x}{2\mu} + \frac{1}{2}\mu\omega_-^2 x^2 + \frac{p^2_X}{2M} + \frac{1}{2}M\omega_+^2 X^2 \,.
\end{equation}
where $\displaystyle\mu = \frac{m_1m_2}{M}$ and $\omega_+^2=\omega^2$ and $\omega_-^2 = \omega^2 + 2\kappa/\mu$.
The wave functions are
\begin{equation}
\Psi_{mn}(x_1,x_2) = \frac{1}{\sqrt{\pi x_0X_0}}\frac{e^{-x^2/2x_0^2}}{\sqrt{m!\,2^m}}\frac{e^{-X^2/2X_0^2}}{\sqrt{n!\,2^n}}H_m\left(\frac{x}{x_0} \right)H_n\left(\frac{X}{X_0} \right) \,,
\end{equation}
where $x=x(x_1,x_2)$ and $X=X(x_1,x_2)$ and $\displaystyle x_0=\sqrt{\frac{\hbar}{\mu\omega_-}}$ and $\displaystyle X_0=\sqrt{\frac{\hbar}{M\omega_+}}$.
How does all this work using matrices "formalism"? And how to extend it to $N$ CQHO?
Ultimately, I would like to redemonstrate (8) and (13) in http://arxiv.org/pdf/hep-th/9303048.pdf
Best Answer
I think it might be conceptually easier to start backwards. Let us assume we have $N$ decoupled harmonic oscillators: $$ H = \sum_{a=1}^N \left[ \frac{(\pi_a)^2}{2m} + \frac{1}{2} m \omega_i^2 \phi_a^2 \right] $$ Now imagine applying the unitary transformation $x_j = {U_j}^a \phi_a$ and $p_j = {U_j}^a \pi_a$ where $U^T U = 1$ in order that both $[x,p]$ and $[\phi,\pi]$ have the standard commutation relation. Such a transformation will act simply on the $\pi_i^2$ term, but the $\omega_i$'s in the second contribution to $H$ will lead to interesting quadratic couplings between the $x_j$. (If you further want unequal masses (I'm not sure you really do in the context of entanglement entropy of a free scalar field), you can perform the scaling $x_i \to \sqrt{m_i} x_i$ and $p_i \to p_i / \sqrt{m_i}$ without changing the canonical commutation relation between the $x$'s and $p$'s.)
Now the question is about a very particular nearest neighbor coupling between the $x_i$. The usual way to think about this transformation involves Fourier series. If I'm not mistaken, something like $$ x_j = \sum_{a=1}^N e^{2 \pi i j a /N} \phi_a $$
should act to diagonalize your Hamiltonian if I further make the assumption that $x_{N+1} = x_1$, i.e. I'm living on a circle.
In the context of the Srednicki paper and entanglement entropy, there are some more modern approaches involving two-point functions. You might want to look at section 2.2 of the review http://arxiv.org/abs/0905.2562 or section III of http://arxiv.org/abs/0906.1663.