Evaluating path integral I am having some trouble remembering how to evaluate path integrals involving multiple particles. Suppose that I have two interacting particles with Lagrangian
$$L= \frac{1}{2}m \dot y^2-\frac{1}{2}m \omega_0^2 (y-b)^2+\frac{1}{2}M \dot x^2-\frac{1}{2}M\omega_1^2 (x-a)^2 -\frac{1}{2} \mu \Omega^2 (x-y)^2,$$
where $a$ and $b$ are constants. Since the resulting path integral is Gaussian, I imagine it is exactly solvable, but I am nevertheless confused about how to proceed. Could someone please show me how to do this-- if you just show me how to do the path integral over the $y$ coordinates, I should be able to figure the rest out.
 A: There are, of course, many ways you could go about evaluating path integrals. The process described below is just an example of how I would proceed. To make it easier to read this answer quickly, here's an outline:


*

*In the introduction, I explain how the problem of "evaluating path integrals" for Gaussian Lagrangians essentially reduces to the problem of computing the time-ordered two-point correlation function.

*In the "solution" section, I explain how to compute the time-ordered correlation function for a general Lagrangian of the type that you describe.


Introduction
The first step in any physics problem is to decide what you want to calculate. Your question is very generic, but since the underlying distribution is Gaussian it is possible in this case to give a completely general answer. The naive goal is to find some 'universal' quantity that can at least formally be used to solve any other problem that we think of. 
A good candidate from the path-integral approach is the time-ordered two-point correlation function, $\langle 0 |T [\hat x(t_2) \hat x(t_1)]|0\rangle$, where $|0\rangle$ is the ground state of the system, and $T[O(t_2)O(t_1)]\equiv O(\max(t_2,t_1))O(\min(t_2,t_1))$. Note that the time-ordering operator $T$ is only defined in the Heisenberg picture, so the expression $T(O_1 O_2)$ doesn't make sense until $O_1$ and $O_2$ are associated with specific times $t$ and $t'$. The time-ordered correlation function is just one choice of elementary operator, that turns out to give full information about Gaussian systems. In your choice of elementary operator, however, you could have for example chosen to omit the time-ordering symbol. In a sense this would give you a larger class of elementary operators to work with, but there are subtleties that make it tricky to use the extra operators in the path integral formalism.
To see how the time ordered two-point correlation function (shortened to $\langle x_{t_2} x_{t_1}\rangle$) is useful, consider the basic problem of computing transition amplitudes of the form $A(t)=\langle\phi|\psi(t)\rangle=\langle\phi|U(t,0)|\psi(0)\rangle$, with $t>0$. Now, it is reasonable to assume that $|\psi(0)\rangle$ can be expanded as a $\hat x$-dependent power series $P(\hat x)$ applied to the ground state $|0\rangle$. For example, in the case of the one dimensional harmonic oscillator, we have $\hat x^n|0\rangle\propto a^{\dagger n}|0\rangle=|n\rangle$. The $\hat x$ operator is associated with a time in the Heisenberg picture, so we have $|\psi(0)\rangle = P(\hat x(0))|0\rangle$. Similarly, $|\phi\rangle = P_\phi(\hat x(t))|0\rangle$ for some other power series $P_\phi(\hat x)$. The full transition amplitude is then given by $\langle 0| (P_\phi(\hat x(t)))^\dagger P(\hat x(0))|0\rangle$, which can be decomposed into a sum of the form $\sum_{k_1,k_2} c_{k_1,k_2}\langle 0|\hat x(t)^{k_2}\hat x(0)^{k_1}|0\rangle$. Now each expectation value $\langle 0| x(t)^{k_2} x(0)^{k_1}|0\rangle$ can be computed readily in the path integral formalism, because they are already time-ordered. However, when the Lagrangian is Gaussian, these expectations can be decomposed even further using an extremely powerful algebraic property of Gaussian integrals, known as Wick's theorem: assuming $\langle \vec x\rangle_0=0$, we have
\begin{align*}
\langle x_{\alpha_1}x_{\alpha_2}\cdots x_{\alpha_n}\rangle_0 = \sum_{j>1} \langle x_{\alpha_1}x_{\alpha_j}\rangle_0\langle\prod_{k>1,k\neq j} x_{\alpha_k}\rangle_0
\end{align*}
(Here $\langle \cdot \rangle_0$ denotes an expectation value with respect to an arbitrary (possibly 'signed') Gaussian measure. This identity can be applied recursively to decompose any moment into a product of two-point correlation functions).
Hence, for Gaussian Lagrangians a good place to start is with the two-point correlation function.
Solution
As a first step, it helps to rewrite your Lagrangian in a more symmetric form. 
Setting $v=\left[\matrix{a \\ b}\right]$ and $q=\left[\matrix{x \\ y}\right]$, $\Lambda=\left[\matrix{m\omega_0^2 & 0 \\ 0 & M\omega_1^2}\right]$ and $u=\Omega\sqrt\mu\left[\matrix{1\\-1}\right]$ the potential energy term can be written as 
\begin{align*}
\frac{1}{2}((q-v)^T \Lambda (q-v)+q^T(uu^T)q)=&\frac{1}{2}(Q^T \Lambda Q+(Q+v)^T(uu^T)(Q+v)),\quad q=Q+v,\\
=&\frac{1}{2}(Q^T (\Lambda+uu^T) Q+2(u\cdot v)(v\cdot Q) + (v\cdot u)^2)
\end{align*}
To proceed, we neglect the last term $(v\cdot u)^2$ since both $v$ and $u$ are constant. 
By setting $\Lambda+uu^T\equiv \Gamma$, we can complete the square in the above expression to obtain 
\begin{align*}
\frac{1}{2}([Q+u\cdot v \Gamma^{-1}v]^T \Gamma [Q+u\cdot v \Gamma^{-1}v]),
\end{align*}
where we have neglected another constant term ($(u\cdot v)^2 v^T\Gamma^{-1} v$. Note that both $\Gamma$ and $\Gamma^{-1}$ are symmetric matrices, so that $v^T\Gamma^{-1}=(\Gamma^{-1}v)^T$). Finally, we set $\eta = Q+u\cdot v\Gamma^{-1}v$, to obtain $\frac{1}{2}(\eta^T\Gamma\eta)$. 
Now to express the Lagrangian in terms of $\eta$, we recall that $\eta = Q+u\cdot v\Gamma^{-1}v=q-v+u\cdot v\Gamma^{-1}v=q+\text{Const}$, so that $\dot{\eta}=\dot{q}$, and the kinetic term is of the form $T=\frac{1}{2}\dot{\eta}^T K \dot{\eta}$, with $K=\left[\matrix{m & 0 \\ 0 & M}\right]$. In physics, it's common to replace this with $-\eta^T K \ddot\eta$, assuming that $\eta$ 'vanishes' on the boundaries. In terms of $\eta$, your action is given by 
\begin{align*}
S&=\int dt\bigg[\frac{1}{2}\eta^T K\ddot\eta - \frac{1}{2}\eta^T\Gamma\eta\bigg]\\
&\equiv \langle \eta , \mathcal{S} \eta\rangle.
\end{align*}
Two-point correlation functions is computed by introducing a fictitious external driving force $h(t)$ that couples linearly to $x$, and 'completing the square':
\begin{align*}
\langle \eta,\mathcal S \eta\rangle+\langle\eta,h\rangle = \langle \eta + \frac{1}{2}\mathcal S^{-1}h, \mathcal S\eta+\frac{1}{2}h\rangle-\frac{1}{4}\langle h, \mathcal S^{-1} h\rangle.
\end{align*}
The two-point function is then typically computed using a generalization of the identity $\langle x_1 x_2 \rangle = \frac{1}{Z}\partial_{\alpha_1}\partial_{\alpha_2}\int d\mu(\vec x) e^{\alpha_1 x_1+\alpha_2x_2}\Big|_{\vec\alpha=0}$, which generalizes in the path integral setting to $\langle x(t_1) x(t_2)\rangle = \frac{1}{\mathcal Z}\frac{\delta}{\delta h(t_1)}\frac{\delta}{\delta h(t_2)}\int d\mu(x(t)) e^{\int dt h(t) x(t)}\Big|_{h(t)=0}$ (here $Z$ and $\mathcal Z$ are partition functions for the measures $\mu(x)$).
The result, you can check, is $\langle \vec\eta(t_1)\vec\eta(t_2)\rangle = \mathcal S^{-1}(t_1,t_2)$ (the integral kernel associated with the operator $\mathcal S^{-1}$). 
To invert $\mathcal S^{-1}$, it is easiest to exploit the time-translation invariance of $S$, and look at how $\mathcal S^{-1}$ acts on the Fourier transform of $\eta(t)$. In Fourier space, $\mathcal S = \frac{1}{2}\omega^2 K+ \frac{1}{2}\Gamma$. In your example, $K$ and $\Gamma$ are both symmetric $2\times 2$ matrices, so the inverse of $\omega^2 K+\Gamma$ can be computed analytically in $\omega$. 
Once you have inverted $\mathcal S$ in Fourier space, you need to know where the singularities (if any) of $\mathcal S^{-1}(\omega)$ occur, so that you can use these with contour integration to invert the Fourier transform and recover $\mathcal S^{-1}(t_2-t_1)$. Once you have obtained $S^{-1}(t)$, you can go to town using Wick's theorem to compute whatever (polynomial) expectation values you desire.
