Relation between the propagator and probability for the infinite well This may be an easy question, but I am really confused about it. 
For the infinite square well, the (time-dependent) energy eigenfunctions are (inside the well):$$\psi_n(x,t) = \sqrt{2/L}\:e^{-iE_nt/\hbar}\:sin\left(\frac{n\pi}{L}x\right)$$
with $E_n = \frac{n^2\pi^2}{2mL^2}\hbar^2$ the eigenvalues of energy, and $L$, the width of the well. So, the probability to find a particle of energy $E_m$ between $x = a$ and $x = b$ at time $t$ is given by Born's rule:$$P(a,b;t) = \int_{a}^{b}\psi_m(x)^*\psi_m(x)dx$$
This probability could be understood as the probability for a particle to be found between $a$ and $b$ at some time, and every $|\psi_m(x)^*\psi_m(x)|^2$ is a squared probability.
On the other hand, the propagator is the amplitude for the particle to travel $a\rightarrow b$ in time $t_a\rightarrow t_b$: $$Propagator = \langle x_b,t_b\rvert x_a,t_a\rangle
= \langle x_b,t_b\rvert(\sum_m\rvert m\rangle\langle m\rvert)\rvert x_a,t_a\rangle = 
\sum_m\psi_m(x_b,t_b)^*\psi_m(x_a,t_a) = \sum_m e^{iE_m(t_b-t_a)/\hbar}\psi_m(x_b)^*\psi_m(x_a)$$
So, my question is: if the propagator is an amplitude, then squaring it should give a probability. However, squaring equation (3): $$P_{a\rightarrow b} = |\langle x_b,t_b\rvert x_a,t_a\rangle|^2 = \left|\sum_m e^{iE_m(t_b-t_a)/\hbar}\psi_m(x_b)^*\psi_m(x_a)\right|^2$$
which obviously is not a probability, since a member like $|\psi_m(x)^*\psi_m(x)|^2$ is actually a squared probability. Then, how can I get a probability from the propagator?
 A: 
which obviously is not a probability, since a member like $|\psi_m(x)^*\psi_m(x)|^2$ is actually a squared probability. Then, how can I get a probability from the propagator?

Right, but only one of your x's is a variable and contributes to the probability amplitude. Recall that your initial state is $|x_0,t_0\rangle$, so if we wish to expand in eigenstates we write (allow me to take $t_0 = 0$)
$$
|\Psi(x_0, 0)\rangle = |x_0,t_0\rangle= \sum_m |m\rangle\langle m|x_0,t_0\rangle 
$$
We can similarly expand the final state in eigenstates
$$
|\Psi(x, t)\rangle = |x,t\rangle= \sum_n e^{iE_nt}|n\rangle\langle n|x,t\rangle 
$$
so that if we take the overlap the sums will collapse to one
$$
\langle x,t\rvert x_0,t_0\rangle = \sum_m e^{iE_mt/\hbar}  \langle x ,t|m\rangle\langle m|x_0,t_0\rangle=  \langle x ,t| \sum_m e^{iE_mt/\hbar} |m\rangle\langle m|x_0,t_0\rangle
$$
The point is that the "extra" terms aren't part of your density, they're constant with respect to the final position. 

Also, be careful. This won't give the probability. What will give the probability is 
$$P(a\leq x \leq b)= \int_{a}^b |\langle x,t\rvert x_0,t_0\rangle |^2 dx $$
A: Your calculations are perfectly correct--so I will directly address the basic issue you are raising. In particular, why is the dimensionality of $\langle x_b,t_b|x_a,t_a\rangle$ not so as for its modulus square to have the dimensionality of a linear probability density, i.e., $L^{-1}$. The issue is not sensitive to the picture that we use, so I will use the Schrödinger picture. 
The reason is that in a continuous space, position eigenstates are not normalizable and are Dirac-normalized. The position eigenstates are normalized as $\langle x|y\rangle=\delta(x-y)$. More explicitly, if a state is localized at a position $x_0$, its wavefunction, by definition, is $\psi(x)=\langle x|\psi\rangle=\langle x|x_0\rangle=\delta(x-x_0)$. Thus, the dimensionality of the wavefunction itself is that of an ordinary linear probability density. This is not an inconsistency precisely because such a wavefunction has been recognized to be non-normalizable and is defined as a Dirac-normalized wavefunction. 
Of course, the the usual normalization and Dirac-normalization are not the same and should not be thought of as just two different flavors of essentially the same thing. The fact that the position eigenstates are Dirac-normalized has an imprint on (and is essential for) their relation to the usual normalizable states. A usual normalizable state, when expressed as a linear combination of other such states, the linear combination takes the form of a summation. Whereas, when a usual normalizable state is expressed as a linear combination of position eigenstates, the linear combination takes the form of an integration. Here the measure of integration invariably comes with dimensionality of its own and the dimensionality of the position eigenstates becomes essential in making the dimensionality of the linear combination the same as the dimensionality of the usual normalizable state. Notice that all of this business would fall apart if the position eigenstates hadn't been Dirac-normalized--because since the complete set of the position eigenstates is continuous, a generic normalizble state would have to be expressed as an integration of different position eigenstates (and not as a usual summation of them) and thus, the dimensionality of the measure of integration would make the equations inconsistent without the "unusual" dimensionality of the position eigenstates (which they owe to their Dirac-normalization). 
Finally, the direct probabilistic interpretation of a non-normalizable state is as such forbidden. They are nonetheless important as normalizable states can be represented as linear combination of such non-normalizable states. See: https://ocw.mit.edu/courses/physics/8-05-quantum-physics-ii-fall-2013/lecture-notes/MIT8_05F13_Chap_01.pdf (Section $2$). One can often invent certain clever ways of thinking about probabilistic interpretations of non-normalizable states--often via thinking of the non-normalizable state as a limit of a normalizable state. For example, the Dirac-delta can be seen as a limit of a Gaussian. Similarly, as such non-normalizable momentum eigenstates are often thought of as the limit of normalizable discrete set of momentum eigenstates on a lattice, etc. 
