In canonical quantisation, it is taught that the propagator for the Klein-Gordon field is defined as $$\Delta_F(\vec x - \vec y) \equiv \left < 0 \right | \overleftarrow{\mathcal T} \phi(\vec x) \phi(\vec y) \left | 0 \right > = \overleftarrow{\mathcal T} \int \frac{\text{d}^3\vec p}{(2\pi)^3 2E_{\vec p}} e^{i \vec p \cdot (\vec x - \vec y)},$$ where $\overleftarrow{\mathcal T}$ is the time ordering operator. It is shown that this is equivalent to the integral form, $$I_F(x-y) = \int \frac{\text{d}^4 p}{(2 \pi)^4} \frac{i}{p^2-m^2 + i\epsilon} e^{ip \cdot(x-y)},$$ which is often proved using the contour integral trick: $$I_F(x-y) = \int \frac{\text{d}^3\vec p}{(2 \pi)^4} i e^{i\vec p \cdot (\vec x - \vec y)} \int \text{d} p_0 \frac{e^{-ip_0(x-y)^0}}{(p_0)^2 - E_{\vec p}^2 + i \epsilon},$$ which is evaluated by residues to show that it is equivalent to the above form $\Delta_F(\vec x - \vec y)$. Textbooks will often state that $I_F(x-y)$ is a Green's function of the Klein-Gordon equation by showing that, $$\left [ \partial^2 - m^2 \right ] I_F(x) = \delta(x).$$
My question is: is there a deeper reason / more general proof that the Green's function for the free field's equation of motion gives the propagator? In this case with the Klein-Gordon field, to me it feels like a freak coincidence that the contour integration trick tells us that the Green's function is the propagator, and I don't see why this is true in the general case? Can this be explained using canonical quantisation and without path integrals?
For example: for the EM field with Lagrangian, $\mathcal L = -\frac{1}{4} F^{\mu\nu}F_{\mu\nu} - \frac{1}{2 \xi} \left ( \partial_\mu A^\mu \right )^2$, it is not hard to show that in Fourier space, the Green's function of the equations of motion is given by $$G^{\mu\nu}(p) = \frac{1}{p^2} \left ( g^{\mu\nu} - (1-\xi) \frac{p^\mu p^\nu}{p^2} \right ),$$ but I really struggle to see how this is related (or equal) to the propagator, $$D^{\mu\nu}(\vec x - \vec y) = \left < 0 \right | \overleftarrow{\mathcal T} A^{\mu}(\vec x) A^{\nu} (\vec y) \left | 0 \right >.$$ Proofs I have seen that this is the propagator for the EM field (without choosing a value of $\xi$) rely on the path integral approach.
