Can you control the oscillation of f(x) as x increases? If you can show that the ratio of f(x) to your 'simplified' form is 'slowly varying' then your asymptotics will probably work out.
A simple example of what you cannot afford is a log-periodic oscillation; this is because the limits of oscillation of the function and Laplace transform need not agree. The simplest example of a log-periodic oscillation is a complex exponential:
$\int_{0}^{\infty }e^{-st}t^{\alpha +i\beta }dt=\left[ \frac{\Gamma \left( \alpha +i\beta +1\right) }{\Gamma \left( \alpha +1\right) }e^{i\beta \log t}\right] \frac{\Gamma \left( \alpha +1\right) }{s}s^{-\alpha }$
In a sense you can view the imaginary part of the exponential as a 'wobbly constant' which changes more and more slowly. The point is the amplitude of the wobble in the transform depends on β but not for the function.
If you do have this sort of problem (it happens all the time in analysis of algorithms and chaotic dynamics) then you can for example resort to the 'gamma function method' of DeBruijn.
The same thing holds true for the moment question. If you look up a counterexample for the moment problem, (e.g. Feller volume II p. 227) you see the ubiquitous log-periodic oscillation.
Not surprisingly the log-periodic oscillation also shows up in convergence questions of Fourier series, but there it is not oscillating more and more slowly, but faster and faster.
After taking a quick look at the paper, I agree with Robert Israel.
Almost no detailed justification is given (perhaps not surprising for a Physics journal...), and it seems to be not so easy to justify.
Take the special case $N=1$ for simplicity; so up to constants, we consider (note the use of (A3) in the paper, this constant is subtracted from the original $\tilde{I}$)
$$\tilde{I}(s) = \frac{(1-\exp(-Ks/(s+r))}{s/(s+r)} - (1-\exp(-K))
$$
for $K$, $r$ constants. Now consider what happens as $|s| \to \infty$, then we get [assuming, dangerously, that my calculations are error free]
$$\tilde{I}(s) \sim -Kr e^{-K}/(s+r) \sim A/s$$
as $|s| \to \infty$, for some constant $A$. The contour is originally a vertical line, but $1/s$ is not in $L^1$ so the integral is not a standard Lebesgue integral (i.e., is not absolutely convergent). However, it is is $L^2$, i.e. is square integrable, so you can use the $L^2$ theory of the Fourier transform to give a meaning.
Now, if you deform the original vertical line contour over { Re(s) = c } by Cauchy's Theorem (push part of it to the left), you get the integral of $\tilde{I}(s) e^{st}$ over the circle $C$ plus an error term, being the sum of integrals over the contours
C1 = { $L + iy : |y| < R $ },
C2 = { $x \pm i R : L < x < c $ },
C3 = { $c + iy : |y|>R$ }
where $c$ is fixed; now let $L \to -\infty$ and $R \to +\infty$ appropriately to get the error term going to zero.
[More detail: as $R \to \infty$ the integral over C3 goes to 0 uniformly in $t$, by the $L^2$ properties of the Fourier transform; also the integral over C2 goes to zero because $|1/s| < 1/R$ on C2.
The integral over C1 goes to zero as $L \to -\infty$ because $| e^{st} / s | < e^{Lt} $ on C1, which clearly $\to 0$ rapidly, at least for $t>0$. ]
There are still 2 things to be justified in this special case: (i) what about $t<0$, and (ii) how do we know the analytic function $\tilde{I}(s)$ really can be represented as the Laplace transform of something?
I think (ii) follows from a general result in my last paper (about Laplace transform representation theorems, which is in Documenta Mathematica 2010, or on my website); but (i) I am not sure about.
[Strangely enough, the contours C1, C2, C3 are exactly the ones I used in my paper, although I got them from an earlier paper by C.Batty and M.D.Blake; and I suspect the original use of these contours dates back at least 50 years, since similar stuff has been used in Tauberian theory and analytic number theory for a long time].
So, even in this special case with just one singularity, there are a lot of extra details to fill in if you want to make it properly rigorous. Since the original paper is from a Physics journal, I would guess it's not the best place to look for rigour(!) Approach with caution...
Best Answer
The answer depends on the class of functions $\phi(t):(0,\infty)\to\mathbb R$ where you want to define the Laplace transform. A standard assumption is that $$e^{-ct}\phi(t)\in L^2(0,\infty)\tag{1}\label{1} $$ for some $c\in \mathbb R$. In this case the Laplace transform $$f(s)=\int_{0}^{\infty}e^{-st}\phi(t)dt\tag{2}\label{2}$$ can be extended analytically to the right half-plane $\{s=\sigma+i\tau:\ \sigma>c\}$. Moreover, it is easy to check that $$\sup\limits_{\sigma>c}\int_{-\infty}^{\infty}|f(\sigma+i\tau)|^2d\tau<\infty.\tag{3}\label{3}$$ Now rewriting \eqref{2} as $$f(\sigma+i\tau)=\int_{0}^{\infty}e^{-it\tau}e^{-\sigma\tau}\phi(t)dt,$$ we observe that $f$ is just the Fourier transform of the function $e^{^{-\sigma t}}\phi(t)$ (trivially extended by $0$ to $t\leq 0$) belonging to $L^2(\mathbb R)$ for $\sigma=c$ and to $L^1(\mathbb R)\cap L^2(\mathbb R)$ for $\sigma>c$. Taking the inverse Fourier transform, we get that $$ e^{-\sigma t}\phi(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{it\tau}f(\sigma+i\tau)d\tau,\qquad t>0,$$ and $$0=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{it\tau}f(\sigma+i\tau)d\tau,\qquad t<0,$$ or, equivalently, $$\lim\limits_{T\to\infty}\frac{1}{2\pi i}\int_{\sigma-iT}^{\sigma+iT}e^{st}f(s)ds=\begin{cases} \phi(t), & t>0 \\\ \\\ 0, & t\leq 0. \end{cases} $$
One can show also that the Parseval identity
$$\frac{1}{2\pi}\int_{-\infty}^{\infty}|f(\sigma+i\tau)|^2d\tau=\int_{0}^{\infty}e^{-2\sigma t}|\phi(t)|^2dt$$ holds, so there is a complete analogy with the standard Fourier transform.
Executive summary. A function $f$ is the Laplace transform of some function $\phi$ satisfying condition \eqref{1}, if and only if it can be extended analytically to the right half-plane $\{s=\sigma+i\tau:\ \sigma>c\}$ and condition \eqref{3} holds. This class of functions is known as the Hardy space on a (right) half-plane.