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