What is the sense of introducing generating functional to the summands of expansion of S-matrix? Let's have generating functional $Z(J)$:
$$
Z(J) = \langle 0|\hat {T}e^{i \int d^{4}x (L_{Int}(\varphi (x)) + J(x) \varphi (x))}|0 \rangle , \qquad (1)
$$
where $J(x)$ is the functional argument (source), $\hat {T}$ is the chronological operator, $\varphi (x)$ - some field. 
I want to understand the reasons for its introduction for the summands of expansion of S-matrix. As I read in the books, it helps to consider only the vacuum expectation values​​, forgetting about in- and out-states. But in $(1)$ appear summands like $\int \frac{J(p)dp}{p^2 - m^2 + i0}$ instead of the contributions from external lines. It may refer to the internal lines. So what to do with them and are there some other reasons to introducing $(1)$ except written by me?
 A: The primary utility in introducing the generating functional is in using it to compute correlation functions of the given quantum field theory.
Let's restrict the discussion to that of a theory of a single, real scalar field on Minkowski space, and let $x_1, \dots, x_n$ denote spacetime points. Of central importance are time-ordered vacuum expectation values of field operators evaluated at such points;
\begin{align}
  \langle0|T[\phi(x_1)\cdots\phi(x_n)]|0\rangle.
\end{align}
It can be shown that these objects can be obtained from the generating functional by taking functional derivatives with respect to the $J(x_i)$ as follows:
\begin{align}
   \langle0|T[\phi(x_1)\cdots\phi(x_n)]|0\rangle = \frac{1}{Z[0]}\left(-i\frac{\delta}{\delta J(x_1)}\right)\cdots \left(-i\frac{\delta}{\delta J(x_n)}\right)Z[J]\Bigg|_{J=0}.
\end{align}
This standard fact is proven in many books on QFT.  It's often proven using the path integral approach which makes it pretty transparent why it's true.  The crux of the argument is that every time you take a functional derivative with respect to the source $J(x_i)$, it pulls down a factor of the field $\phi(x_i)$.  Dividing by $Z[0]$ is an important normalization relating to vacuum bubbles, and setting $J=0$ after computing the appropriate functional derivatives eliminates terms with more than $n$ factors of the field and renders the final result source-independent as it should be.
A: If you can calculate vacuum-to-vacuum transition amplitudes, you can calculate S-matrix elements, because the two are related by the LSZ reduction formula. The LSZ will in any case chop off the propagators for external lines that the generating functional inserts, so you will end up only needing to compute amputated diagrams.
