Are asymptotic states in scattering experiments really momentum eigenstates? In a typical collider experiment, two particles, generally in approximate momentum eigenstates at $t=-\infty$, are collided with each other and we measure the probability of finding particular outgoing
momentum eigenstates at $t=\infty$.
Firstly, what does it mean for the particles to be in approximate momentum eigenstates? Does it mean that there is a small spread in the distribution of the momentum of the particles, due to the particles being in a superposition of momentum eigenstates, so that the particles are rather in approximate momentum eigenstates, so to speak?
Secondly, how can we be certain that the experiment produces only momentum eigenstates at $t=\infty$ and not some superposition of momentum eigenstates at $t=\infty$?
 A: The experiment certainly does produce a very general complex superposition of momentum eigenstates. The spread is not "small" in any way – virtually all allowed (by conservation laws etc.) final states are represented in the superposition for any initial state. We detect particles of particular momenta in the final states because the detectors (e.g. at the LHC) perform a measurement, just like everywhere in quantum mechanics, and this measurement collapses the final wave function (a very general superposition, indeed) to an eigenstate of all the particles' momenta (and before that, an eigenstate of the particle number and other operators that determine the list of the particles). The Born rule applied to the probability amplitudes quantifies the probability or probability density that the collapse will end up with one set of values of the final momenta or another.
The calculation of the S-matrix elements, $S_{fi}$, is a calculation of all the coefficients $S_{fi}$ that appear in this superposition of states, in front of the states $f$ with particles of different momenta in the final state, assuming that the initial state is $i$. That's why the S-matrix elements are functions of the (initial and final) momenta.
Also, the initial state is not a true momentum eigenstate in the real world. It can never be because the momentum eigenstates are not normalizable to one. That doesn't matter. The point is that the S-matrix $S_{fi}$ is computed for a basis of states $i$ of the initial Hilbert space and a basis of states $f$ of the final Hilbert space.
The probability amplitude that a general superposition in the initial state $u$ evolves to a general superposition in the final state, the vector $v$, may be computed by a simple matrix multiplication of all the coefficients, $v\cdot S\cdot u$. Only the $S$-matrix middle factor in this product is "hard" – it has to be computed by summing all the Feynman diagrams etc.
But the matrix elements $S_{fi}$ for momentum eigenstates contain all the relevant dynamical information about the behavior of the particles or other physical systems. For example, some totally universal formula allow us to calculate the cross sections as some universal kinematic coefficients (that depend on momenta and energies) times $|S_{fi}|^2$. Only the matrix elements $S_{fi}$ depend on the coupling constants and the theory-specific features of the interactions, not the universal kinematic coefficients.
When the initial momentum has a big uncertainty, the cross section will be some combination or weighted average of the cross sections calculated for different values of the momenta.
Concerning the final states, the probabilities or cross sections are written as "densities" for particular values of the final momenta. In the case of cross sections, these densities are known as the "differential cross sections" (well, the "differential cross section" is just the density with respect to the angular parts of the final momentum vectors; the lengths i.e. radial parts are determined by the conservation laws). The total cross section (basically a probability up to the coefficients determining the precision and intensity of the beams etc.) that something happens may be computed as an integral of the "differential cross sections" over the (angular part of the) momenta of the final particles – over the region in the momentum space we're interested in (where a certain condition is fulfilled – when we calculate the probability that the condition will be obeyed by the scattering).
I just mentioned the final states – $|S_{fi}|^2$ times some simple factors determine the differential cross section which may be integrated over angles to get the overall cross section. The kinematic treatment of the initial state is completely different – basically due to the arrow of time. In practice, the cross sections are never integrated over the initial states (initial particles' momenta). Instead, the fact that the cross section "could be integrated" means that the formula proportional to $|S_{fi}|^2$ doesn't directly quantify the probabilities but rather the cross section. 
The multiplicative difference between a probability and a cross section exactly takes care of the "momentum spread" of the initial state. Note that the smaller the momentum spread of the initial state is, the larger the position spread is, by the uncertainty principle. So a small initial momentum uncertainty unavoidably makes the probability of an interaction/collision small, too. At the LHC, the protons are confined to a tube whose width is a fraction of a millimeter or so. It means that $\Delta x_{\rm transverse}$ is rather small which implies that there's some nonzero $\Delta p_{\rm transverse}$. The value of the latter is still vastly smaller than $|\vec{p}|$ of the particles so this uncertainty of the momentum may be neglected because no one can measure the dependence of the cross sections on the momenta this accurately.
