Virtual particles and S-matrix One of methods of introducing of virtual particles is using perturbation theory. We say that scattering matrice amplitude $M_{in \to out}$ contains of $\delta(P_{out} - P_{in})$, which realizes energy-impulse conservation and relativistic connection $E^{2} - p^{2} = m^{2}$. When we leave only a few terms of the expansion this leads to unconservation energy at the interaction stage. But in final we obtain that energy-impulse is conserved. 
So my question is following: how exactly the dropping of terms of higher order by the coupling constant leads to energy-impulse non-conservation and how exactly do we obtain that in final energy-impulse is conserved?
Also there is a question about delta-function. Whether it is administered solely by requiring the conservation of energy-momentum? Or it arises naturally from the very formulation of the problem?
 A: It is better to begin with amplitudes calculus in position space, because things become more clear (usual "Feynman diagrams" correspond to momentum space).
Suppose, for instance, a massless scalar field theory, with a $\phi^4$ interaction. In this theory, each vertex has 4 "legs", which may be real particles or field perturbations/field correlations (also unfortunately called "virtual particles")
Now, look at the corresponding process : $2$ in-particles "beginning" at $x_1, x_2$, joining at a first vertex in $z$, then two field perturbations beginning at $z$ and ending at a new vertex $w$, plus $2$ out-particles beginning at $w$ and ending at points $x_3$ and $x_4$.
We are interested in the amplitude of the process $A(x^1,x^2,x^3,x^4)$. The amplitude for a field perturbation beginning at $x$ and ending at $x'$ is $D(x-x')$, the propagator. The amplitude for a real particle (the in- and out- particles) is $\delta(x-x')$
So, we have : 
$A(x^1,x^2,x^3,x^4) \sim \int dz ~dw ~\delta(x^1 - z)\delta(x^2 - z) [D(z-w)]^2 \delta( w - x^3 )\delta(w - x^4) \tag{1}$
(here all the integrations and $\delta$ functions are $4$-dimensional)
Now, look at the Fourier transform : $\tilde A(p_1,p_2,p_3,p_4)\sim \int dx^1dx^2dx^3dx^4 ~e^{-i(p_1x^1 +p_2x^2+p_3x^3+p_4x^4)}A(x^1,x^2,x^3,x^4)$
We suppose here that the in and out particles are real particles, so that $p_i^2=0$ $1 \leq i\leq 4$ (for simplicity we do not put explicitely the $\delta(p_i^2)$ functions, this is just a global term)
You get : 
$\tilde A(p_1,p_2,p_3,p_4) \sim \int dz dw\int e^{-i((p_1 +p_2) z -(p_3+p_4)w)}[D(z-w)]^2 \tag{2}$
Let $u=z+w, v=z-w$, so that $z =\frac{u+v}{2}, w =\frac{u-v}{2}$, we have : 
$\tilde A(p_1,p_2,p_3,p_4) \sim \int e^{-i((p_1 +p_2) +(p_3+p_4))\large \frac{v}{2}}[D(v)]^2 \int du  e^{-i((p_1 +p_2) -(p_3+p_4))\large \frac{u}{2}} \tag{3}$
The integration on $u$ gives a $\sim \delta((p_1 +p_2) -(p_3+p_4))$, so you have :
$\tilde A(p_1,p_2,p_3,p_4) \sim \delta(p_1 +p_2- p_3-p_4)\int dv ~e^{-i(p_1 +p_2)v}[D(v)]^2 \tag{4}$
Let $\tilde D(p)$ be the Fourier transform of the propagator $D(x)$, we have finally : 
$\tilde A(p_1,p_2,p_3,p_4) \sim \delta(p_1 +p_2- p_3-p_4)\int dp \tilde D(p) \tilde D(p_1+p_2-p) \tag{5}$
Now, we are in momentum space, the domain of Feynman diagrams.
We see the difference between real particles and fields perturbations. Real particles, (in and out) have a on-shell momenta $p_i^2=0$, and there is a global conservation of momentum-energy.
There is also a conservation of momentum/energy between real particles and field perturbations, this may be seen, in the last expression in the term $\tilde D(p)\tilde D(p_1+p_2-p)$. At the first (momentum) vertex, the total momentum is $p_1+p_2$, and there is a repartition of the momentum into the two legs (fields perturbations). But you may see that the integration is on all $p$, so there is no mass-shell conditions for the fields perturbations. Field perturbations are not particles. 
[EDIT]
If you are not satisfyed with formula $(1)$,it could be obtained from the LSZ reduction formula, which states, roughly :
$$A(x^1,x^2,x^3,x^4) \sim \square_{x_1}\square_{x_2}\square_{x_3}\square_{x_4} \int dz ~dw ~D(x^1 - z)D(x^2 - z) [D(z-w)]^2 \\
D( w - x^3 )D(w - x^4) \tag{6}$$
The connection with the formula $(1)$ is done with $\square_{x_1}D(x^1 - z) = \delta(x^1 - z)$
