Expanding in terms of energy eigenfunctions is not necessary.
The defining property of G (the function, not the operator) is, as I said, that if \psi(x,t) is a solution to Schrodinger's equation, then \psi(x,t) = \int dx_0 G(x,t, x_0, t_0) \psi(x_0, t_0). G(x,t,x_0, t_0) has a kind of self-referential property:
G(x,t,x_0, t_0) = \int dx_1 G(x,t, x_1, t_1) G(x_1, t_1, x_0, t_0)
(where t > t_1 > t_0)
So you can keep expanding, so that for any N > 0,
G(x,t,x_0,t_0) = \int dx_1 dx_2 ... dx_N \Pi_{j=0}^{N} G(x_{j+1}, t_0 + (j+1) \delta t, x_j, t_0 + j \delta t)
where \delta t = \frac{t-t_0}{N+1} and x_{N+1} = x.
This is an exact equality. But to go further toward the path integral, we need an approximation for G(x, t+\delta t, x', t) for when \delta t is small. I'm going to number the steps from here:
1. Express G(x, t+\delta t, x', t) in terms of the delta function.
So here is an expression for G for small \delta t:
G(x, t+\delta t, x', t) \approx e^{-i H(x, p, t) \delta t} \delta(x-x')
If H were time-independent, then this would be valid for all \delta t, not just when \delta t is small.
2. Use the delta function to get rid of the x-dependence of the Hamiltonian
Now, the next approximation is that since \delta(x-x') is zero away from x=x', we can replace x by x' in the exponential:
e^{-i H(x, p, t) \delta t} \delta(x-x') \Rightarrow e^{-i H(x', p, t) \delta t} \delta(x-x')
3. Use the integral representation of the delta function.
The next trick is to use the integral representation of the delta-function:
\delta(x-x') = \frac{1}{2\pi} \int dk e^{i k (x-x')}
So we need to evaluate:
\frac{1}{2\pi} e^{-i H(x',p, t) \delta t} \int dk e^{i k (x-x')}
We can bring the operator e^{-i H(x',p,t) \delta t} inside the integral:
\frac{1}{2\pi}\int dk e^{-i H(x',p, t) \delta t} e^{i k (x-x')}
4. Replace the operator p by its eigenvalue.
Since e^{i k (x-x')} is an eigenstate of momentum, we can replace p by its eigenvalue, k. This gives us:
\frac{1}{2\pi}\int dk e^{-i H(x',k,t) \delta t} e^{i k (x-x')} = \frac{1}{2\pi}\int dk e^{i (k \frac{x-x'}{\delta t} - H(x',k,t)) \delta t}
This is now a pure function expression, H(x',k,t) is only a function, not an operator.
5. Use method of stationary phase to evaluate the integral.
Let's expand the expression in the exponential in a power series, centered on k=k_0 (we'll decide k_0 a little later):
k \frac{x-x'}{\delta t} - H(x',k,t) \approx (k_0 \frac{x-x'}{\delta t} - H(x', k_0, t)) + (k-k_0) (\frac{x-x'}{\delta t} - \frac{\partial H(x,p,t)}{\partial p}|_{x=x', p=k_0}) + \frac{(k-k_0)^2}{2} \frac{\partial^2 H(x,p,t)}{\partial p^2}|_{x=x', p=k_0} + ...
The method of stationary phase chooses k_0 so that the linear term, proportional to k-k_0, vanishes. This means choosing k_0 so that
\frac{x-x'}{\delta t} - \frac{\partial H(x,p,t)}{\partial p}|_{x=x', p=k_0} = 0
6. Use classical Hamiltonian equations to rewrite the expression in terms of velocity and Lagrangian
Here's where the first miracle happens: The classical (non-quantum) Hamilton equations of motion are:
\frac{dx}{dt} = \frac{\partial H}{\partial p}
\frac{dp}{dt} = - \frac{\partial H}{\partial x}
So the stationary phase condition is just that \frac{x-x'}{\delta t} = v, where v is the classical velocity associated with the Hamiltonian H(x,p,t) whenx=x', p=k_0.
The quadratic term is proportional to \frac{\partial^2 H(x,p,t)}{\partial p^2}|_{x=x', p=k_0}, which can be written as \frac{\partial v}{\partial p}|_{p=k_0, v = \frac{x-x'}{\delta t}}. For a simple nonrelativistic Hamiltonian of the form \frac{p^2}{2m} + V(x,t), v = \frac{p}{m}, so \frac{\partial v}{\partial p} = \frac{1}{m}. I'm going to assume that that's the case.
Here's where the second miracle happens: The constant term in the Taylor expansion is:
k_0 \frac{x-x'}{\delta t} - H(x', k_0, t)
For our choice of k_0, this is the same as
(p v - H(x,p,t))|_{p = k_0, v = \frac{x-x'}{\delta t}, x = x'}
That's the classical expression for the Lagrangian L(x,v,t) in terms of the Hamiltonian H(x,p,t). So putting everything together, our exponent can be written as
[ L(x,v,t)|_{x=x', v= \frac{x-x'}{\delta t}} - \frac{1}{2m} (k-k_0)^2 + ...] i \delta t
So we can write (keeping only up to the quadratic term):
G(x,t+\delta t, x', t) \approx \frac{1}{2m} e^{i L(x,v,t) \delta t} \int dk e^{-i \frac{1}{2m} (k-k_0)^2 \delta t}
We still have an integral, and it's actually divergent, but if we pretend that \delta t has a small negative imaginary component, then we can evaluate it as a Gaussian: \int dk e^{-i \frac{1}{2m} (k-k_0)^2 \delta t} = \sqrt{\frac{2\pi m}{i \delta t}}. So our final expression for G:
G(x, t+\delta t, x', t) \approx \sqrt{\frac{m}{2 \pi i \delta t}} e^{i L(x',v,t) \delta t}
Putting this back into our expression for G(x, t, x_0, t_0) gives:
G(x,t,x_0, t_0) \approx (\sqrt{\frac{m}{2 \pi i \delta t}})^N \int dx_1 dx_2 ... dx_N \Pi_{j=0}^N e^{i L(x_j, v_j, t_j) \delta t} (where v_j = \frac{x_{j+1} - x_j}{\delta t} and where t_j = t + j \delta t)
= (\sqrt{\frac{m}{2 \pi i \delta t}})^N \int dx_1 dx_2 ... dx_N e^{i \sum_{j=0}^N L(x_j, v_j, t_j) \delta t}
Since \psi(x,t) = \int dx_0 G(x,t, x_0, t_0) \psi(x_0, t_0), we can immediately write an expression for \psi:
\psi(x,t) = (\sqrt{\frac{m}{2 \pi i \delta t}})^N \int dx_0 dx_1 dx_2 ... dx_N e^{i \sum_{j=0}^N L(x_j, v_j, t_j) \delta t} \psi(x_0, t_0)
To the extent that the limit N \Rightarrow \infty makes sense, we can write this as:
\psi(x,t) = \int \mathcal{D}[z(t)] e^{i \int_{x_0}^x L(z, \frac{dz}{dt}, t) dt} \psi(x_0, t_0)