The "two words" answer is that the $S$-matrix is Poincare covariant. The Poincare covariance in particular requires (independently on the number of particles in initial and final states, as well as on the detils of the theory) the proportionality of the $S$-matrix to the delta-function $\delta(p_{1}+...+p_{n}-p_{1'}-...-p_{n'})$, where $\{p_{i}\}$ are momenta of initial state particles and $\{p_{i'}\}$ are momenta of final state particles. This can be qualitatively understood by taking into account that the subgroup of the Poincare group - the translation group - leads to conservation of the stress-energy tensor, which in particular requires 4-momentum conservation. That's why for any amplitude within Feynman diagram approach we use the delta-function.
The text below just contains the derivation of this statement.
Suppose the S-matrix:
$$
\tag 1 S_{\alpha \beta} = \langle \alpha, \text{out}|\beta, \text{in}\rangle
$$
The QFT you mentioned in the question is based on the Poincare symmetry. In particular this means that we need to postulate that in the space of any state $|\alpha, \text{out}\rangle$ and $|\beta,\text{in}\rangle$ is realized the unitary transformation of the Poincare group equivalent to a representation in the Fock space. This means in particular that (and same for $|\beta,\text{out}\rangle$)
$$
\tag 2 |\alpha,\text{in}\rangle = |\mathbf p_{1},\sigma_{1},...,\mathbf p_{n},\sigma_{n}\rangle,
$$
where $\{p_{i}, \sigma_{i}\}$ defines the one-particle state with given 4-momentum $p_{i}$ corresponding to fixed energy-momentum orbit (say, massive or massless particles) and helicity $\sigma_{i}$ (assuming the particle has spin $s_{i}$). The Poincare group transformation law for the state $(2)$ is
$$
\tag 3 |\alpha,\text{out}\rangle \to \sqrt{\frac{(\Lambda p_{1})^{0}}{p_{1}^{0}}...\frac{(\Lambda p_{n})^{0}}{p_{n}^{0}}}e^{ia_{\mu}\Lambda^{\mu}_{\ \nu}(p_{1}+...p_{n})^{\nu}}\times
$$
$$
\times \sum_{\tilde{\sigma}_{1},\tilde{\sigma}_{2},...}D_{\tilde{\sigma}_{1}\tilde\sigma_{1}}^{s_{1}}...D_{\tilde{\sigma}_{n}\sigma_{n}}^{s_{n}}|(\Lambda p_{1}),\tilde{\sigma}_{1},...,(\Lambda p_{n}),\tilde{\sigma}_{n}\rangle
$$
Here $a_{\mu}$ and $\Lambda_{\mu}^{\ \nu}$ are defined through the Poincare group transformation acting on the arbitrary 4-vector $x_{\mu}$,
$$
\tag 4 x_{\mu} \to \Lambda_{\mu}^{\ \nu}x_{\nu} + a_{\mu},
$$
and $D_{\sigma'\sigma}$ is the unitary Lorentz group transformation corresponding to the fixed orbit of the 4-vector. The same is true for $|\beta,\text{in}\rangle$.
Since the transformation $(3)$ is unitary, this means that $S$-matrix $(1)$ is covariant, i.e.,
$$
S_{\alpha\beta} \to S_{\alpha\beta}'= S_{\alpha\beta}
$$
In particular for the transformation $(4)$ corresponding to $\Lambda = \mathbf{1}$ and arbitrary $a_{\mu}$ one obtains
$$
S_{\alpha\beta}' = e^{ia_{\mu}(p_{1}+p_{2}+...p_{n}-p_{1'}-p_{2'}-...-p_{n'})}S_{\alpha\beta} = S_{\alpha\beta}
$$
For non-trivial scattering, this equality requires zero $S_{\alpha\beta}$ unless
$$
p_{1}+p_{2}+...+p_{n}-p_{1'}-p_{2'}-...-p_{n'} = 0
$$
But this means nothing but the proportionality of the $S$-matrix $S_{\alpha \beta}$ to the delta-function $\delta(p_{1}+...p_{n}-p_{1}'-...-p_{n}')$,
$$
S_{\alpha\beta} \equiv T_{\alpha\beta}\delta(p_{1}+...+p_{n} - p_{1'}-...-p_{n'})
$$
P.S. You can find more on the Poincare covariance of the $S$-matrix (as well as the derivation given above) in Weinberg's QFT, Vol. 1. I also think that the derivation similar to this can be found in Schwartz's book.
