Multiplying propagators The amplitude for a particle going from $x$ to $y$ is $G(x,y)$.
So why isn't the amplitude for going $x$ to $y$ to $z$ 
$$ G(x,z) \neq \int G(x,y)G(y,z) dy^4 $$ 
but instead
$$ G(x,z) = \int G(x-y) (\Box_y +m^2) G(y-z) dy^4 = -(\Box_z +m^2) \int G(x,y)G(y,z) dy^4 $$
I know this comes from the action. But in terms of probabilities how does this make sense in simple terms? What is the physical meaning of this operator in the middle? i.e. can it be thought of as the amplitude for the particle turn turn round and head the other way?
 A: $\newcommand{\bphi}{\boldsymbol{\phi}}
\newcommand{\bD}{\mathbf{D}}
\newcommand{\phiu}{\bphi_\textrm{unforced}}
\newcommand{\force}{\mathbf{f}}
\newcommand{\phis}{\bphi_s}
\newcommand{\bG}{\mathbf{G}}$
To answer this question, I will first review some background information and introduce some compact notation. Then I will discuss why the correct expression is correct. Then I will discuss why the incorrect expression is incorrect.
Introduction/Review
You can figure this out just by looking at the situation classically. In classical field theory, the concepts of a wave operator and a green's function are important. I will review those concepts in turn. 
Wave operator
Classically, you have a field $\bphi$ which, when no external forcing is appplied, satisfies a wave equation $\bD \bphi=0$, where $\bD$ is the wave operator, which is some differential operator, in your case it is given by $\bD =\Box + m^2$. It is important to note that there are many solutions $\phiu$ to the wave equation $\bD \bphi =0$.
Green's function
Now suppose that the field $\bphi$ does feel a force $\force$. Then the field will satisfy the wave equation $\bD \bphi = \force$. Then you can always find a solution $\phis$ to the equation, but this solution is not unique, because if $\bD \phis = \force$, then you also have $$\bD \left(\phis + \phiu\right) = \bD \phis + \bD\phiu = \bD \phis + \mathbf{0}  =\force,$$ so $\phis + \phiu$ is also a solution whenever $\phiu$ is an unforced solution.
Now it is often the case that we are given an $\force$, and we want to find a $\phis$. One tool to do this is called the green's function. The green's function is one solution to the wave equation when the force $f(y)$ is given by a delta function $\delta(y)$. The green's function is denoted $G(y)$. Notice that there are many choices of green's function since solutions to the wave equation are not unique. But having picked a particular green's function (that is, a particular solution $\bD \bG = \boldsymbol{\delta}$, we can write the solution for the wave equation for arbitrary forcing ($\bD \bphi=\force$) as $\phi_s(x)= \int G(x-y)f(y) dy$, which I will write as $\phis = \bG * \force$. 
Why is the correct expression correct?
To explain why the expression is correct I will first prove it abstractly using symbols, then I will show how the proof works out with an example. 
Proof
Now that we understand the wave operator $\bD$ and green's functions $\bG$, we can ask why it is that $$ G(x,z) = \int G(x-y) (\Box_y +m^2) G(y-z) dy^4, $$ which in my notation is written $\bG = \bG * \bD \bG$. To think about this, it helps to consider a more general question, which is given a solution $\phis$, how can we obtain another solution $\phis'$, which experiences the same force $\force$ as $\phis$ (in symbols, $\bD \phis' = \force = \bD \phis$)? One way to do this is to apply $\bD$ to $\phis$ to measure the force $\force$ producing $\phis$, and then use our green's function procedure to produce another solution $\phis'$ corresponding to this force: $\phis' = \bG * \force = \bG*\bD \phis$. Notice that the right hand side depends on the choice of $\bG$, so unless we are very lucky, we will not have $\phis'$ be the same as $\phis$. 
But there is one case where $\phis'$ will be equal to $\phis$, and that is when $\phis$ itself came from our green's function $\bG$ (in symbols, $\phis = \bG * \force$ for some force $\force$). Then since we still have $\phis' = \bG * \force$, we find that $\phis' = \phis$. 
Now consider in particular the case where $\phis$ does come from a green's function, and in fact it comes from apply the greens function to a delta function force (in symbols, $\phis = \bG*\boldsymbol{\delta} = \bG$. Then, as just explained, we have $\phis' = \phis$. Now remembering how we originally found $\phis'$ by measuring the force and applying our green's function method to the force, the equation $\phis' = \phis$ becomes $\bG * \bD \phis = \phis$. Then plugging in our particular choice $\phis = \bG$, we obtain $\bG * \bD \bG = \bG$, which is what we wanted to show. 
Example
To make this more concrete, it would help to look at an example. Consider a forced simple harmonic oscillator with frequency $\omega$, where the (zero-dimensional) field $\bphi$ represents the oscillator position, and the wave operator is $\bD = \partial_t^2 + \omega^2$. What is the greens function? Well if we hit the oscillator with a hammer at $t=0$, we expect it to oscillate with frequency $\omega$, so $G(t) = 0$ for $t<0$ and $G(t) = \sin(\omega t)$ for $t>0$, up to an overall constant. (This constant is $\dfrac{1}{\omega^2}$, so our green's function is actually $\dfrac{1}{\omega^2} \sin(\omega t)$).
Now let's go through the argument with this example. We take our green's function $\bG = \dfrac{1}{\omega^2} \sin(\omega t)$ for $t>0$ and $\bG =0$ for $t<0$, and we find out what force created it. We do this by applying the wave operator given by $\bD = \partial_t^2 + \omega^2$. Upon applying the operator we find that the motion satisfies the wave equation for $t \ne 0$, but that at $t=0$, something funny happened, specficially a $\delta$-function force was applied. Now that we have this $\force$, we can ask what motion we get when we apply our green's function to this force. So we do the convolution $ \int \dfrac{1}{\omega^2} \sin(\omega(t-t')) \delta(t') dt'$ and we get the answer $\dfrac{1}{\omega^2} \sin(\omega(t))$, which is exactly the green's function we started with, so everything checks out. 
Why is the incorrect expression incorrect?
Now we can also ask the question what is wrong with $\bG = \bG*\bG$? Well notice first that the integral implied by the convolution '$*$' goes over all time. This is probably not what you meant. You probably meant $G(t_3,\vec{x}_3;t_1,\vec{x}_1)= \int G(t_3,\vec{x}_3;t_2,\vec{x}_2) G(t_2,\vec{x}_2;t_1,\vec{x}_1) d \vec{x}_2$, where only the intermediate spatial coordinate $\vec{x}_2$ is integrated over and not the time $t_2$. 
There is an obvious problem with this, which is that $\bG$ converts from force to field configuration, not field configuration to field configuration. So the equation  would fail on dimensional grounds. 
Notice, however, that even besides those two points, this has no chance of working because the future configuration $\phi(t_3)$ of the field depends not only on the past field configuration, but also on the past field momentum. Now you may counter that "ok, so there are many future field configurations corresponding to a given past field configuration, so my $\tilde{\bG}$ operator will not be unique, but your $\bG$ also wasn't unique, so what is the big deal?" Well you could define a non-unique $\tilde{G}$, but there would be no way for it to satisfy $\tilde{G}(t_3-t_1) = \tilde{G}(t_3-t_2)\tilde{G}(t_2-t_1)$. 
To see this, let's go back to the example of the simple harmonic oscillator. The natural choice of $\tilde{G}$ is $\cos(\omega t)$. As stated above, we needed to make a choice about the past ($t=0$) momentum of the oscillator, and here we chose the momentum to be zero. Now suppose we pick some past time $t_1$ and future time $t_3$. Then $\tilde{\bG}(t_3 -t_1) = \cos(\omega(t_3-t_1))$, Now suppose we pick some intermediate time $t_2$. Does $\tilde{G}(t_3-t_1) = \tilde{G}(t_3-t_2)\tilde{G}(t_2-t_1)$? Well no it doesn't. What goes wrong? We first we can compute $\tilde{G}(t_2-t_1)$ and we get $\cos(\omega (t_2-t_1))$, which is what we were supposed to get for the amplitude at the intermediate time $t_2$. But now comes the crucial problem. When we evolve this forward by multiplying by $\tilde{G}(t_3-t_2)$, we have to make an assumption about the momentum. Our choice was to assume the momentum is zero, but in this case it is not. This means we wind up getting $\cos(\omega(t_3 - t_2))\cos(\omega t_2)$ instead of $\cos(\omega(t_3-t_1))$ like we were supposed to. This is the basic problem with trying to get $\tilde{\bG} = \tilde{\bG}*\tilde{\bG}$.
A: The Feynman propagator is the green's function of the klein gordon equation therefore $G(x,z) = \int G(x-y) (\Box_y +m^2) G(y-z) dy^4$ is equivalent to $G(x,z) = \int G(x-y) \delta(y-z) dy^4$ if you use the delta function to perform the integral you get $G(x,z)=G(x-z)$. Now as far as intuition goes I would say that the reason that you cannot simply multiply the propagators is that you would then get the the probability that the the particle went from z to y and from y to x whereas the Feynman propagator only gives the probability that the particle went from z to x independent of whether or not it went through y. This is why you need the delta function to enforce that y and z are the same point.
