How to find the QM propagator for more general quadratic theories? (and why didn't my attempts work?) The Problem
I'm working on a problem using the Path Integral approach to quantum mechanics. Actually, more accurately I'm interested in the Wick rotated version of the problem, so I suppose technically it's a 1D statistical field theory, but that's not the point. The point is that I have integrated out one field in the problem, so that the effective Lagrangian for the remaining field, $x$, despite being quadratic, has non trivial energy dependance. It is presented simplest in Fourier space, where the action is of the form:
$$S[x]=\int \frac{d\omega}{2\pi} \frac{a(\omega)}{2b(\omega)} |x(\omega)|^2$$
where $a$ and $b$ are both polynomials in $\omega$. What I'm trying to compute is (equivelant to) the QM propagator (in the sense of the kernel) for this kind of theory.
My first instinct, since it's a quadratic Lagrangian, was to find the classical action given some boundary conditions. I think I'm correct in saying that this is
$$x(\omega)=\sum_ka_k\delta (\omega - \omega_k)$$
where $\omega_k$ are the roots of $a(\omega)$, and $a_k$ are chosen such that $x(t)$ goes through some set of points. Since this is basically a higher derivative theory, I was of course expecting to need to specify more points than just the start and end, but I'm struggling with the step of substituting this back in. It looks like I'd have to integrate the square of some delta functions multiplied by a function that is zero at all the points that the deltas are non-zero?
What Else I've Tried
I tried leaning into the SFT analogy. In that picture, the probability of a path is
$$P(path)\propto exp(-S[path]).$$
Because I have a quadratic Lagrangian, this looks suspiciously like a multivariate normal distribution over a continuum of variables corresponding to the different frequencies. The precision matrix would be $\frac{a(\omega)}{b(\omega)}2\pi \delta(\omega-\omega')$, and this is diagonal so the the covariance matrix would be $\frac{b(\omega)}{a(\omega)}2 \pi\delta(\omega-\omega')$ (=the propagator in the other sense!).
My understanding was that the probability that the path taken coincides with a few observations at a few points, i.e. $x(t_k)=x_k$, could be found up to a normalisation by marginalising out all other points in the multivariate normal. This is well known to just involve dropping them from the covariance matrix (and mean vector) so all I would have to do is Fourier transform back:
$$\Sigma(t)=FT^{-1}[\frac{b(\omega)}{a(\omega)}](t).$$
And I would find that $x_k$ would be distributed in a multivariate normal distribution with vanishing mean, and covariance matrix $\Sigma_{jk} = \Sigma(t_j-t_k)$, but I've tried this out, and I don't think it works. For instance, in the free particle case, I find that $\Sigma(t)\propto|t|$, so given only two points, I would expect
$$\Sigma_{jk} \propto \left(\begin{array}{cc}
0 & \Delta t\\
\Delta t & 0
\end{array}\right)$$
and $$\Sigma^{-1}_{jk} \propto \left(\begin{array}{cc}
0 & -1/\Delta t\\
-1/\Delta t & 0
\end{array}\right).$$
The well known heat kernel instead implies $P(xy)\propto exp(-\frac{m(x-y)^2}{2\Delta t})$. To it's credit, this means my method got the cross terms right, but fails for the variances. This seems to be true for the other examples I tried too, like the harmonic oscillator.
Any idea how to approach this properly, and why my attempt failed?
 A: You can write your action as
\begin{equation}
S = \frac{1}{2} \int  \frac{d\omega}{2\pi} \int \frac{d\omega'}{2\pi} x(\omega') K(\omega, \omega') x(\omega)
\end{equation}
where
\begin{equation}
K(\omega,\omega') = \frac{a(\omega)}{b(\omega)} \delta(\omega-\omega')
\end{equation}
This can be formally rewritten as a matrix equation, where $x$ is a vector and $K$ is a matrix
\begin{equation}
S = \frac{1}{2} x^T K x
\end{equation}
The propagator is then the inverse of $K$. One way to see this is by computing the path integral with a
\begin{equation}
Z[J] = \int Dx e^{i \left(S[x] + \int d\omega x(\omega) J(\omega)\right)} = e^{-\frac{i}{2} J^T K^{-1} J}
\end{equation}
where $K^{-1}$ appears in the place that the propagator appears in a quadratic theory.
So the propagator is
\begin{equation}
G(\omega) = \frac{b(\omega)}{a(\omega)}
\end{equation}
As an example, a standard theory would have $b=1$ and $a=\omega^2+m^2$. Using the formula above would yield the standard propagator, $G(\omega)=\left(\omega^2+m^2\right)^{-1}$.
However, you should be aware that for general $a(\omega)$ and $b(\omega)$, you can easily run into issues such as having a non-local theory, Ostragradski ghosts, or violating unitarity bounds on amplitudes such as the Froissart bound.
