Why is the energy-momentum tensor $T_{\mu \nu}$ the source of the gravitational field instead of the momentum vector $P_{\mu}$? In Newton's gravity theory, the mass $m$ is the source of the gravitational field.
In Special Relativity mass can change depending on the observer, so we construct a Lorentz four-vector $P_{\mu}$ that is conserved under Lorentz transformations.
General Relativity is, by hypothesis, locally-Minkowski, so mass still behaves as a component of a four-vector and it would not make sense if only the mass sourced the gravitational field, as then gravitation would be different for different inertial observers, violating the principle of relativity of Einstein.
However, I cannot understand from a physical point of view why we should go any further than taking $P_{\mu}$ as the source of the field.
I know that the all the properties of a manifold are reflected in the metric, a 2-tensor with components $g_{\mu \nu}$, and in particular all the curvature properties are collected in the Riemann 4-tensor $R^\rho_{\sigma \mu \nu}$, obtained from second derivatives of the metric. I also know that the only symmetric (tosion-free) conserved tensor that can be written that involves only second derivatives of the metric is Einstein's tensor, $G_{\mu \nu}$.
I guess that because the Riemann tensor has 4 indexes, we cannot obtain an object with a single index by contracting it with other Lorentz-invariant tensors, but then why should the equations involve only this 2-tensor object instead of the full 4-tensor Riemann? And what is the physical rationale that justifies using $T_{\mu \nu}$ as the source? Why not something like $P_{\mu} P_{\nu}$?
 A: Because for a continuous distribution of matter (or energy), such as a fluid, what you really want is the density of momentum, and this density is given by the energy-momentum tensor.
If you have a fluid (or, again, any distribution of energy), you can define at each point the velocity $\mathbf{u}$, which in spacetime generalizes to the four-velocity $u^\mu$. No problem so far. But if you want the momentum at each point, that's the mass times the velocity. And we run into a problem, because the mass of a single point is zero! Mass is not a quantity that can be defined at each point; it's not a field.
That's fine, because we know that we really should have used the mass density. And so in Newtonian mechanics we can write the momentum density as $\rho \mathbf{u}$, and the energy density as $\frac12 \rho u^2$. But in relativity this doesn't really work, because density is not a scalar: it changes after a Lorentz boost. Some math reveals that the energy density transforms as the zero-zero component of a tensor (picking up two factors of $1/\sqrt{1-v^2}$), which is the energy-momentum tensor.
Or to put it another way: the energy-momentum density, whatever it is, must have two indices. You need one index to tell me what your time axis is; and then I can take the 3D space orthogonal to that time axis, and give the momentum density in that space, accounting for the other index. But I cannot define a density if you don't tell me how you choose to split spacetime into space and time.
A: Here is one line of reasoning: Recall that a source $J_k$ is (minus) the derivative of the effective action $\Gamma[\phi_{\rm cl}]$ wrt. the classical field $\phi^k_{\rm cl}$. By analogy, one may consider the Hilbert/metric SEM tensor $$T^{\mu\nu}~:=\pm\frac{2}{\sqrt{|g|}}\frac{\delta S}{\delta g_{\mu\nu}}$$ as the source of the gravitational field $g_{\mu\nu}$, cf. OP's title question.
A: Actually we do use something a bit like $P_\mu P_\nu $. If you sum the expression  $P_\mu v_\nu $ for all the particles of matter in a region, you get $ T_{\mu\nu} $, at least for the matter element. You still need the e.m. component from photons, but you can see that it must be included because of conservation of energy-momentum.
Einstein had realised much earlier that curvature must be involved, by thinking of the twin paradox which lead him to express the equivalence principle. One cannot use Riemann because that would mean spacetime was flat in regions with no mass. He tried using Ricci, but it did not work. Only when he found that the Einstein tensor obeys the contracted Bianci identity, and consequently automatically conserves energy-momentum did he find the correct equation.
A: The starting point of the rationale should be the Newton's classical gravitation field equation. This is not Newton's famous law of gravitation, but the Poisson-equation with the mass density on the RHS:
$$\Delta \phi = 4\pi \rho$$
This equation is much more general and allows to compute the gravitational potential for any kind of mass distribution expressed by the mass distribution density $\rho(\vec{r})$, in particular for $\rho(\vec{r})= m\delta^3(\vec{r}-\vec{r}')$ for a mass $m$ at $\vec{r}'$. 
The EFE (Einstein's Field equations) should in any case contain the Poisson-equation in Newtonian approximation. 
If we use the EFE in the following form (let's put speed of light $c=1$):
$$R^i_k = 8\pi G (T^i_k - \frac{1}{2}\delta^i_k T^l_l)$$
We limit the analysis to the $^0_0$ component:
$$R^0_0 = 8\pi G (T^0_0 - \frac{1}{2} T^l_l)= 4\pi G(T^0_0-T^1_1-T^2_2-T^3_3)\approx 4\pi G\rho  $$
using $T_{ik}=diag(\epsilon, p,p,p)$ where $\epsilon \equiv \rho$ if $c=1$ is the energy density. This is the energy-momentum tensor of an incompressible liquid at rest. 
Moreover we assume that in most cases pressure $p$ can be neglected as source of gravitation in this approximation
(in the Newtonian theory the pressure at least does not serves as source of gravitation).
So we already reproduce the source term (RHS) of the Poisson-equation. Trying to use the energy-momentum 4-vector or binomials of it would, however, fail to construct such an expression.
The analysis of the LHS is a bit more involved. The best approach is via the geodesic deviation equation. It is well-known that in free fall the felt acceleration towards the massive source can be completely eliminated, whereas extended probe bodies  still feel tidal forces that cannot be eliminated. The formal description is done via the equation of geodesic deviation containing components of the Riemann-tensor:
$$\frac{d^2 n^i}{ds^2} =-R^i_{0j0} n^j$$
which describes the acceleration of a  normal vector $n^i$ between two adjacent geodesics. In GR this tidal force is associated with the curvature of space, expressed in the equation by the Riemann-tensor.
Actually something comparable can also be done in Newtonian theory: 
For two adjacent mass points on the orbits: $x^i(t)$ and $x^i(t)+n^i(t)$ we get 2 equations of motion:
$$\ddot{x^i}(t)= - \frac{\partial\phi}{\partial x^i}\mid_{x(t)}$$
and 
$$\ddot{x^i}(t)+\ddot{n^i}(t) = - \frac{\partial\phi}{\partial x^i}\mid_{x(t)+n(t)}$$
Taking the difference of both equations yields: 
$$\ddot{n^i}(t)= - \frac{\partial^2\phi}{\partial x^i\partial x^j}n^j(t)$$
So we can establish the following correspondence:
$$ R^i_{0j0} \leftrightarrow \partial_i\partial_j \phi$$
If we contract the indices $i$ and $j$ we get:
$$ R_{00} \leftrightarrow \Delta\phi $$
So if replace $ R_{00}$ by $\Delta\phi $ in the $^0_0$ component of the EFE we finally get: 
$$\Delta\phi \approx 4\pi G\rho$$
Again, if instead we had tried it with the uncontracted Riemann-tensor on the LHS, we would have ended up with terms like $\partial_i\partial_j \phi$, which 
don't appear on the LHS of the Poisson-equation. If we were not able to reproduce Newtionian physics in simple circumstances, something would be evidently wrong with the EFEs. 
One could of course wonder if there are more subtle ways to reproduce the Poisson-equation from more subtle gravitational field equations since we need a quantum-compatible theory (replacing GR). This is subject of actual research. 
