First, the term "propagator" is usually defined as the Green's function of the first type, not the second type, i.e. as a solution to the diffential equation $\hat L G = \delta$.
At any rate, those definitions are ultimately equivalent – when the details are correctly written down – because the Green's function defined as the correlator in the second definition obeys the first differential equation.
The differential operator $\hat L$ is what appears in the linearized equations of motion for the field, in this case $\psi(x_1,t_1)$, and it only acts on $\psi(x_1,t_1)$, not $\psi^\dagger (x_2,t_2)$.
The time-ordering operator $T$ may be written in terms of the step function
$$ T(\psi \psi^\dagger) = \psi \psi^\dagger\cdot \theta(t_1-t_2) - \psi^\dagger \psi\cdot\theta(-t_1+t_2)$$
where $\psi$ is always at $x_1,t_1$ and $\psi^\dagger$ is at $x_2,t_2$. Now, ask what happens when you act with $\hat L$ on the right hand side of the displayed equation above.
By the Leibniz rule, there are the terms with $\hat L \psi = 0$. It vanishes by the equations of motion. But there are extra terms where $\hat L$ acts on the step functions.
The operator $\hat L$ contains the term that differentiates with respect to $t_1$ multiplied by a coefficient $C$. This turns $\theta(t_1-t_2)$ to $\delta(t_1-t_2)$. The same occurs in the next term, but with the opposite sign which cancels the sign that was already there. So the extra terms are
$$ \hat L T(\psi \psi^\dagger) = C\delta(t_1-t_2) (\psi \psi^\dagger + \psi^\dagger \psi)$$
I got two terms because there were two terms. However, these two terms exactly combine to the anticommutator of $\psi$ and $\psi^\dagger$ which only needs to be evaluated for $t_1=t_2$, the equal-time anticommutator, and the result is $D\cdot \delta(x_1-x_2)$.
That's why the action of $\hat L$ on the correlator ends up being $CD\delta(t_1-t_2)\delta(x_1-x_2)$ where the constants $C,D$ are mostly just factors of $i$ etc.
For bosonic fields, $\hat L$ has the second derivative with respect to time. One of the derivatives has the fate as above, the other one turns the other $\phi$, which plays the role of $\psi^\dagger$, into $\partial_t \phi$ which is the canonical momentum, and it's the right variable that has the $\delta$-function-like commutator. Also, the intermediate sign is the opposite one but the result is the same, some $CD\cdot\delta\cdot\delta$.
If ${\scr D}$ is some differential operator, a Green's function is a bidistribution $G(x,y)$ with the property that it solves the equation $${\scr D}_xG(x,y)=\delta(x,y)\tag{1}.$$
That is a definition. Now the thing is that to uniquely define $G(x,y)$ one must prescribe extra conditions. In particular, when ${\scr D}=\Box-m^2$ is the Klein-Gordon operator, the standard $i\varepsilon$ prescription is one such condition when solving (1) which has the property that the obtained solution is the time-ordered two-point function $\langle 0|T\{\phi(x)\phi(y)\}|0\rangle$ in free theory. That would be the Feynman propagator. The retarded/advanced propagators follow from other conditions added to (1).
Some authors may define $G(x,y)$ by including some numeric factor to the RHS of (1), but then one is changing the definition, and one must be consistent the whole way through.
So the short answer is that the retarded propagator is one possible Green's function, as is the Feynman propagator. They are all solutions to (1) which differ by some subsidiary condition which picks one out of the space of inverses of ${\mathscr{D}}$.
Best Answer
Let me change notation a little. I'll consider the general non-interacting one-particle Green's function--you can easily generalize to the retarded case if you wish, it's the same basic idea:
$$G_\sigma^0(r_1,\,t_1;\,r_2,\,t_2)=-i\langle \psi_0|\left\{ \Theta(t_1-t_2)c_\sigma(r_1,\,t_1)c_\sigma^\dagger(r_2,\,t_2)-\Theta(t_2-t_1)c_\sigma^\dagger(r_2,\,t_2)c_\sigma(r_1,\,t_1)\right\}|\psi_0\rangle$$
Now, we take a spatial Fourier transform:
\begin{align} G_\sigma^0(k,\,t)& =\int dt \int d^3 r e^{-ik\cdot r}G_\sigma(r,\,t)\notag\\ &=-i\langle \psi_0|\{\Theta(t)c_{k\sigma}(t)c_{k\sigma}^t)-\Theta(-t)c_{k\sigma}^\dagger(t)c_{k\sigma}(t)\}|\psi_0\rangle\end{align}
From the Heisenberg equation of motion, we can take the time dependence out of the operators:
$$ \frac{d}{dt} c_{k\sigma}(t)=\frac{i}{\hbar}[H,\,c_{k\sigma}(t)]$$
I won't give you the full solution, but I will give you the end result:
$$c_{k\sigma}(t)=e^{-\frac{i}{\hbar}\epsilon_k t}c_{k\sigma}$$
With this, the one-particle Green's function is given by
\begin{align} G_\sigma^0(k,\,t)&=-i\langle \psi_0|\left\{ \Theta(t)c_{k\sigma}c_{k\sigma}^\dagger -\Theta(t)c_{k\sigma}^\dagger c_{k\sigma}\right\}|\psi_0\rangle e^{-\frac{i}{\hbar}\epsilon_k t}\notag\\ &=-i\left\{ \Theta(t)(1-n_{k\sigma})-\Theta(-t)n_{k\sigma} \right\}e^{-\frac{i}{\hbar}\epsilon_k t} \end{align}
We now perform a temporal Fourier transform, all the while attaching an extra $e^{-\delta |t|}$ term to ensure good behavior at infinity:
\begin{align} G_\sigma^0(k,\,\sigma)&=\int dt e^{i\omega t-\delta t}G_\sigma(k,\,t) \end{align}
I think you can solve this integral with what I gave you--the solution is
$$G_\sigma^0(k,\,\omega)=\frac{1-n_{k\sigma}}{\omega-\epsilon_{k}+i\delta}+\frac{n_{k\sigma}}{\omega-\epsilon_{k}-i\delta}$$
Where in the above I have taken units where $\hbar$ is unity. Defining $\delta_k=\textrm{sgn}(\epsilon_{k}-\epsilon_F)\delta$ (where $\epsilon_F$ is the Fermi energy), we find a nice concise form of the above:
$$G^0(k,\,\omega)=\frac{1}{\omega-\epsilon_k+i\delta_k}$$
In the interacting system, we find something similar, except now with a term known as the self-energy $\Sigma(k,\,\omega)$:
$$G(k,\,\omega)=\frac{1}{\omega-\epsilon_k-\Sigma(k,\,\omega)+i\delta_{k}}$$
This self-energy can be thought of as a renormalization of the energy that comes about when we turn on interactions. In fact, we can rewrite the interacting Green's function in the form
$$G(k,\,\omega)=\frac{Z_k}{\omega-\widetilde{\epsilon_k}+i\delta_k}$$
where $Z_k$ is the quasiparticle weight
$$Z_k=\frac{1}{1-\frac{\partial \Sigma(k,\,\omega)}{\partial \omega}\bigg|_{\omega=\epsilon_k}}$$
and we have expanded around the pole at $\omega=\widetilde{\epsilon_k}$.
If you want to learn more about the details of the Green's function approach to many body physics, I'd suggest Zagoskin or Dickhoff and Van Neck's books. They are both excellent and up-to-date resources for the subject.
EDIT: If you wish to consider the Green function as an operator, look back to when I first took out the exponential term:
\begin{align} G_\sigma^0(k,\,t)&=-i\langle \psi_0|\left\{ \Theta(t)c_{k\sigma}c_{k\sigma}^\dagger -\Theta(t)c_{k\sigma}^\dagger c_{k\sigma}\right\}|\psi_0\rangle e^{-\frac{i}{\hbar}\epsilon_k t}\notag\\ \end{align}
We can ignore the sandwiched states and perform the temporal Fourier transform immediately to get the Green function in terms of the number operator, which is directly proportional to the non-interacting Hamiltonian.