The reasoning of the definition of $S$-matrix The definition of the $S$-matrix is given by
$$S=\lim_{t_{f}\rightarrow\infty}\lim_{t_{i}\rightarrow-\infty}U(t_{f},t_{i}).$$
Where $U(t_{f},t_{i})$ is the evolution operator, given by the $$U(t_{f},t_{i})=Te^{-i\int_{t_{i}}^{t_{f}}V(\tau)d\tau},$$ which $V(t)$ is the interaction Hamiltonian, $T$ is the time-ordering.
One thing I do not quite understand is the reason for taking $t_{f},t_{i}$ to infinity. What is the purpose of that? The definition makes the integral $\int V(t)dt$ hard to converge due to the upper and lower limit being unbounded. Does there exists solid physical reasoning behind this?
 A: As lurscher already said correctly, the idea is to choose an initial and a final time that is (in time) away from the actual time-dependent interaction in order to make sure that the interaction does not act at the initial and at the final time (well, this actually not completely true: the supposed free particles at $t_i$ and $t_f$ are not really free, in order to treat this correctly the concept of a dressed particle has to be introduced which under the effect of self-interaction has a different mass etc.. But this is out of scope here).
I would even say putting the initial and end time to infinite values makes the evaluation of the integral easier. Because in a further step one can construct Fourier integrals out of it which can be in the most cases evaluated --- and if not QFT disposes of a couple methods of regularization of integrals that will make the integral finite.
Example:
QED in first order of perturbation theory:
$V_I(t)\equiv H_I(t)$ represents the interaction part of the Hamiltonian.
So the S-matrix is up to first order:
$$ S = id - i\int_{t_i\rightarrow-\infty}^{t_f\rightarrow\infty} \hat{H_I} dt$$
As long as there is no derivative in the expression the interaction Hamiltonian is just a volume integral over the negative of the interaction part of the  Lagrange density.
$$ \hat{H_I} = -\int d^3x {\cal{L}} \quad$$
Now in QED the interaction operator is just the volume  integral  of the  electron current $j^\mu$ and the electromagnetic 4-potential  operator $\hat{A}_\mu$
$${\cal{L}} =-e\hat{j}^\mu \hat{A}_\mu$$
so plugging this into the precedent expression gives:
$$S = id - i \int_{t_i\rightarrow-\infty}^{t_f\rightarrow \infty} dt \int d^3x e \hat{j}^\mu  \hat{A}_\mu   =   id - i \int d^4x e \hat{j}^\mu  \hat{A}_\mu  $$
In the last 4-dim. integral it is supposed that the integration borders in time and space are all put to $(\pm)$infinity.
Now we only need to plug in the electromagnetic potential operator into the last expression ($c$ and $c^\dagger$ are the annihilation and creation operators of photons):
The 4-potential operator can be developed in annihilation and creation operators:
$$\hat{A}_\mu = \sum_{\epsilon=1,2}\sum_p [c_{pe} \epsilon_\mu e^{-ipx}  + c^\dagger_{pe} \epsilon^\ast_\mu e^{ipx}   ]$$
$\epsilon_\mu$ are components of the polarisation vector of the photon.
Then we can compute S-matrix elements between an electron state $i$ and an electron state $f$ plus 1 photon like:
$$\langle 1_k f|S| 0i\rangle  =\langle 1_k f| 0i\rangle - i \int d^4x e \langle f|j^\mu|i\rangle \sum_{\epsilon=1,2}\sum_p \langle 1_k| c^\dagger_{pe}| 0\rangle \epsilon^\ast_\mu e^{ipx} =  - i\int d^4x \sum_{\epsilon=1,2} e \langle f|j^\mu(x)|i\rangle e^{ikx} \epsilon^\ast_\mu$$
where $\tilde{j}^\mu_{if}$ is the 4-dim. Fourier transform of the electromagnetic current:
$$\tilde{j}^\mu_{if}(k)  =   \int d^4x \langle f|j^\mu(x)| i\rangle e^{ikx}$$
Imagine the borders of the time integral were not infinite, the integral  in time direction could not be (at least not easily) computed as it would not be a Fourier integral. For more details of the calculation see my post: Why is $\mathcal{M}(k)$ given by this? (Ward Identity derivation in Peskin & Schroeder)
EDIT
Using infinite integration border technique provides even more physics to uncover:
If the initial $\psi_i$ and end states $\psi_f$ in the current $j^\mu$ describe free particles we know that they behave like $\psi \sim e^{-ipx} \sim e^{-iEt} $
$$\langle f |j^\mu(x)| i \rangle = \bar{\psi}^f \gamma^\mu \psi_i \sim e^{iE_ft} e^{-iE_i t}$$
So:
$$\tilde{j}^\mu_{if}(k)  =   \int d^4x \langle f|j^\mu(x)| i\rangle e^{ikx} \sim \int_{-\infty}^\infty e^{iE_ft} e^{-iE_i t} e^{i\omega_\mathbf{k} t}dt \sim \delta(E_f -E_i +\omega_\mathbf{k})$$
You are right if $E_f = E_i + \omega_\mathbf{k}$ the integral diverges. But actually we found that the considered interaction fulfills energy conservation. These kind of delta functions appear systematically in the S-matrix, but disappear once the S-matrix element gets  integrated over the phase space of the out-going particles.  This delta-function appears also in simpler computations as for instance of transition probabilities via Fermi's golden rule. So using infinite integration borders provides us energy conservation.
