Green functions, Propagators and probability amplitude that a particle propagates (The original post has been copiously edited to make it more clear but it still precisely corresponds to what I had intended to ask)
My QFT knowledge has very much rusted and i got confused by these few lines from Peskin and Schroeder:

p.27: " [..] the amplitude for a particle to propagate from $y$ to $x$ is $\langle 0| \phi(x) \phi(y) |0\rangle $. We will call this quantity $D(x — y)$."


(The relation with the commutator is explicitly calculated at (2.53) p.28, + bottom of p.29:
$$ i\hbar\, \Delta(x-y) := [\phi(x) , \phi(y)] = \cdots = D(x-y)- D(y-x) = \langle 0|[\phi(x) , \phi(y)]|0\rangle\tag{2.53}$$
the r.h.s. are implicitly understood as being proportional to the $\mathbb{1}$ operator)

Now the expressions of the retarded and Feynman propagators are given (2.55) p.30

$$D_R := \theta (x^0 -y^0) \langle 0|[\phi(x) , \phi(y)]|0\rangle \tag{2.55}$$

and (2.60) p.31  (without commutators)

$$D_F := \theta (x^0 -y^0) \langle 0|\phi(x) \cdot \phi(y)|0\rangle + \theta (y^0 -x^0) \langle 0|\phi(y) \cdot \phi(x)|0\rangle \tag{2.60}$$

which by definition of "propagator" or "Green's function" satisfy $$(\square +m^2) G(x,y)= -i\delta^4(x-y).$$

My confusion: propagators also seem to have the interpretation of  amplitude of propagation, cf. e.g. wikipedia or Peskin last 2 lines p.82, but the three different functions $D(x-y),\ D_R(x-y),\ D_F(x-y)$ obviously cannot have the exact same interpretation! so what is the interpretation of each?
If the first question is too easy, here is a second: Are propagators always combinations of the $D(x-y)$?

*

*for interacting fields?

*for more general PDEs?

Old Edit: I'm not trying to relate interacting theories to the free one, so $D(x-y)$ stands for the amplitude of propagation in each theory not only in the free one. If propagators are distinct from $D(x-y)$ in general, can they still be expressed in terms of it?

New comments:

*

*the interpretation of $\langle 0| \phi(x) \phi(y) |0\rangle $ as the amplitude for a particle to propagate from $y$ to $x$ assumes quite a lot of structure (as opposed to something that is defined just from the equation): the $\phi$ have to be operators acting in some space with a decomposition in terms of $a_{\mathbf{k}}, a^{\dagger}_{\mathbf{k}}$ so that one can have a "particle" interpretation.

*So in mathematics, when considering a linear PDE, one can just distinguish between "bi-solutions" (solution in each variable) of the homogeneous equation or "bi-solutions" of the equation with a non-zero r.h.s./a source. Now (as opposed to when I asked the question), I have decided to call a propagator any bisolution of the homogeneous equation and Green's function any bisolution of the eq. with $\delta$ on the r.h.s.

*There is a link between propagators and Green functions given by the Duhamel formula (when one interprets PDEs as an infinite system of ODE i.e. ODE with value in a space of functions of spatial variables), cf. e.g. this and this; cf. also this and related posts.

*For the Klein-Gordon equation, there is a 2 dimensional vector space of propagators with basis:
$$ D(x-y)=i\hbar\, \Delta_+ (x-y) =\int_{\mathbb{R}^3} \frac{e^{-i k\cdot (x-y)} \, d^3 \mathbf{k}}{(2\pi)^3\, 2\hspace{.9pt}\omega_{\mathbf{k}}} \quad\text{where}\quad \omega_{\mathbf{k}}=\sqrt{\mathbf{k}^2 + m^2}$$
and
$$ D(y-x)=i\hbar\, \Delta_- (x-y) =\int_{\mathbb{R}^3} \frac{e^{+i k\cdot (x-y)} \, d^3 \mathbf{k}}{(2\pi)^3\, 2\hspace{.9pt}\omega_{\mathbf{k}}} $$
The space of Green functions is an affine space (cf. rule: particular solution + solutions of homogeneous) with the same associated vector space. cf. e.g. Field Quantization (1996), Walter Greiner, Joachim Reinhardt, right after eq. (4.164) p.114 (he calls "commutator fct" what I call "propagator" while the latter is used by him for "Green's fct"; moreover he takes another "basis"); cf. also Introduction to the Theory of Quantized Fields (1976), N.N. Bogoliubov, D.V. Shirkov, § right after eq. (15.11) p.142.

*It is a completely different story to study non-linear PDEs because one cannot use propagators or Green functions to express solutions. What happens in practice is the following: let us consider for example the so-called $\lambda\, \phi^4$ theory (power 4 in the Lagrangian but 3 in the equation...)
$$ \big(\square +\, m^2\big) \phi = - \frac{\lambda}{3!}\, \phi^3 $$
The idea of perturbation is, if $\lambda <<1$, to consider r.h.s. as an exterious "source" (although it does contain the unknown function $\phi$) and treat this equation as a linear one (l.h.s.) with a source term...

 A: Propagators are independent of interactions. They only depend on the free part of the Lagrangian. For example, the KG Lagrangian reads
$$
\mathcal L_{KG}=\frac{1}{2}(\partial \phi)^2-\frac{1}{2}m^2\phi^2-\mathcal H_I(\phi)
$$
where $\mathcal H_I$ could be, eg, $\frac{g}{4!}\phi^4$. There are different definitions of the propagator (all equivalent, of course). For our purposes, the easiest definition is that it is the Green function of the KG PDE. The fastest way to find $D_F$ is to go to momentum space, and "solve" for the propagator:
$$
\mathcal L_{KG}(k)=\phi(k)(k^2+m^2)\phi(k)-\mathcal H_I(k)\equiv \phi(k) D_F^{-1}\phi(k)-\mathcal H_I(k)
$$
where we find
$$
D_F(k)=\frac{1}{k^2+m^2}
$$
(there are several conventions on the normalisation of $D_F$, so you can sometimes find some numerical factors in front of the expression above)
On the other hand, its not difficult to prove that $D_F$ agrees with your expression:
$$
D_F=\langle 0|\mathcal T\ \phi_0(x)\phi_0(y)|0\rangle
$$
where $\phi_0(x)$ is the solution of the free KG equation. If you add interactions, then the amplitude to go from $x$ to $y$ is no longer the propagator. You can call the amplitude $C(x,y)$, and we usually call this new quantity a correlator. In perturbation theory, it is not difficult to prove that $C$ can always be written in terms of the propagator (as an asymptotic expansion over some parameter that weights the interactions). The correlator depends on $\mathcal H_I$, but the propagator does not. You can prove that if $\mathcal H_I=0$, then $C(x,y)=D_F(x,y)$.
With this, we can answer your first two questions: 1) $D_F$ is not the propagation amplitude, but $C(x,y)$ is. They agree when there are no interactions. 2) The propagator is independent of interactions, but the correlator is not. The latter can be written as an asymptotic expansion where the different terms include powers of $D_F$, so that correlators can always we written in terms of propagators.
As for your third question, the propagator does depend on the PDE, as it depends on the free part of the Lagrangian. For example, if you take a different Lagrangian, as Dirac's
$$
\mathcal L_D=\bar\psi(i\not\partial-m)\psi-\mathcal H_I(\psi)
$$
then using the same line of reasoning as before we find that the propagator is
$$
S(k)=\frac{1}{\not k-m}
$$
I hope this answers your questions.
Addendum
Some people make no distinction between the concept of propagator and correlator. They call the former free propagator and the latter exact propagator (or even bare vs dressed propagators). Note that this is not very usual, for most people just call $D_F$ the propagator, and call $C$ correlator (or two-point function). IMO it is very important to understand that these are different objects, so I don't like to use the same word for them.
That being said, it seems that you just call both objects propagator. If this is the case, you should convince yourself that the free propagator is independent of interactions, and the exact one is not. So if you are going to use the same word for both, then at least you should not use the same notation. For example, in Srednicki's book he writes $\Delta$ for the free propagator, and $\boldsymbol \Delta$ for the exact one. By the way, Peskin & Schroeder do call $C$ correlator, and $D_F$ propagator (see page 82).
Anyway, if you see the links I posted, you can see that $D_F(k)=1/k^2+m^2$ (in KG's case), or $D_F(k)=1/\not k-m$ (in Dirac's). Then it is easy to see that there are no coupling constants in these expressions, so they are indeed independent of interactions. By the Dyson series you can calculate the correlators (which do depend on the coupling parameters), which are given, at fist order, by the propagators.
A: Let $P_1(x,y), P_2(x,y)$ and $G(x,y)$ be "bisolutions" of the Klein-Gordon equations (or more general linear PDEs)
$$ \left\lbrace \begin{aligned}
    \big(\square + m^2\big)\, P_1(x,y) &= 0\\
    P_1\left( x^0, \mathbf{x}, x^0, \mathbf{y} \right) & = \delta^{(3)}(\mathbf{x}-\mathbf{y})\\
    \lim_{t\to x_0} \partial_t\, P_1\left( t, \mathbf{x}, x^0,\mathbf{y}\right) & = 0
 \end{aligned}\right. \label{1}\tag{1}
$$
$$ \left\lbrace \begin{aligned}
    \big(\square + m^2\big)\, P_2(x,y) &= 0\\
    P_2\left( x^0, \mathbf{x}, x^0, \mathbf{y} \right) & = 0\\
    \lim_{t\to x_0} \partial_t\, P_2\left( t, \mathbf{x}, x^0,\mathbf{y}\right) & = \delta^{(3)}(\mathbf{x}-\mathbf{y})
 \end{aligned}\right. \label{2}\tag{2}
$$
$$ \left\lbrace \begin{aligned}
    \big(\square + m^2\big)\, G(x,y) &= \delta^{(4)}(x-y)\\
    G\left( x^0, \mathbf{x}, x^0, \mathbf{y} \right) & = 0\\
    \lim_{t\to x_0} \partial_t\, G\left( t, \mathbf{x}, x^0,\mathbf{y}\right) & = 0
 \end{aligned}\right. \label{3}\tag{3}
$$
(As mentionned in here, initial/boundary conditions are required to ensure unicity of a solution). They all have a "propagation" interpretation in the sense that they are the kernel of the linear maps $I_1, I_2$ (i.e. $I_i(f)(x)=\int_{\mathbb{R}^3} P_i(t,\mathbf{x};y^0,\mathbf{y})\, f(\mathbf{y})\, d\mathbf{y}$) and $J$:
$$ I_i : \left\lbrace \begin{aligned}
    E_0\, & \longrightarrow E\\
    f(\mathbf{y}) & \longmapsto I_i(f)(x)
\end{aligned}\right.\enspace i=1\ \text{or}\ 2\qquad J : \left\lbrace \begin{aligned}
    F\; & \longrightarrow E\\
    h(y) & \longmapsto G(h)(x)
\end{aligned}\right.$$
from some chosen space $E_0$ of initial conditions (at time $y^0$) which are functions of the spatial variables $\mathbf{y}$ only and $F$ of source which are functions of space-time $y$ to the space $E$ of functions/distributions of space-time within which one looks for solutions. To $f$ one associates $I_i(f)(x)$, the unique solution of
$$ \left\lbrace \begin{aligned}
    \big(\square + m^2\big)\, I_1(f)(x) &= 0\\
    I_1(f)(y^0,\mathbf{x}) & = f(\mathbf{x})\\
    \lim_{t\to y^0} \partial_t\, I_1(f)\left( t, \mathbf{x}\right) & = 0
 \end{aligned}\right.\quad\text{and}\quad \left\lbrace \begin{aligned}
    \big(\square + m^2\big)\, I_2(f)(x) &= 0\\
    I_2(f)(y^0,\mathbf{x}) & = 0\\
    \lim_{t\to y^0} \partial_t\, I_2(f)\left( t, \mathbf{x}\right) & = f(\mathbf{x})
 \end{aligned}\right. \label{4}\tag{4}
$$
while $G(h)(x)$ is the unique solution of
$$\left\lbrace \begin{aligned}
    \big(\square + m^2\big)\, G(f)(x) &= h(x)\\
    G(h)(y^0,\mathbf{x}) & = 0\\
    \lim_{t\to y^0} \partial_t\, G(h)\left( t, \mathbf{x}\right) & = 0
 \end{aligned}\right. \label{5}\tag{5}
$$
i.e. $P_i$ propagate initial conditions to solutions of the free equations while Green functions propagate initial conditions to solutions of the equation with source. These are linear maps in case the PDE is linear. This fails for non-linear equations, except in the approximations of perturbation theory...
