Just a hint as a possible way to go, through the Beta Function.
$$
\eqalign{
& f_{\,N,\,n} = \sum\limits_{k = 0}^{\min (N,n) - 1}
{{{k!\left( {n - 1 - k} \right)!\left( {N - 1 - k} \right)!} \over {\left( {N + n - 1 - k} \right)!}}} = \cr
& = n\sum\limits_{k = 0}^{\min (N,n) - 1}
{{{\left( {k + 1 - 1} \right)!\left( {n - k - 1} \right)!\left( {n - 1} \right)!\left( {N - k - 1} \right)!}
\over {\left( {n + 1 - 1} \right)!\left( {N + n - k - 1} \right)!}}} = \cr
& = n\sum\limits_{k = 0}^{\min (N,n) - 1}
{{\rm B}\left( {k + 1,n - k} \right){\rm B}\left( {n,N - k} \right)} \cr}
$$
Now, the integral representation of Beta looks to get clumsy.
On the other hand, Beta is rapidly decaying at growing of either of its arguments, so
it might be possible to try by truncating the sum.
Another way is that $F_{\, N}(x)$ actually reads
$$
\eqalign{
& F_{\,N} (x) = \sum\limits_{n = 0}^\infty {\sum\limits_{k = 0}^{\min (N,n) - 1}
{{{k!\left( {n - 1 - k} \right)!\left( {N - 1 - k} \right)!}
\over {\left( {N + n - 1 - k} \right)!}}} \;x^{\,n} } = \cr
& = \sum\limits_n {\sum\limits_k
{{{\left[ {0 \le k < n} \right]\left[ {0 \le k < N} \right]k!\left( {n - 1 - k} \right)!
\left( {N - 1 - k} \right)!} \over {\left( {N + n - 1 - k} \right)!}}x^{\,n} } } = \cr
& = \sum\limits_k {\left[ {0 \le k < N} \right]k!\left( {N - 1 - k} \right)!
\sum\limits_n {{{\left[ {k < n} \right]\left( {n - 1 - k} \right)!}
\over {\left( {N + n - 1 - k} \right)!}}x^{\,n} } } = \cr
& = \sum\limits_{0 \le k < N} {k!\left( {N - 1 - k} \right)!
\sum\limits_{k < n} {{{1^{\,\overline {\,n - 1 - k\,} } }
\over {1^{\,\overline {\,N + n - 1 - k\,} } }}x^{\,n} } } = \cr
& = \sum\limits_{k = 0}^{N - 1} {k!\left( {N - 1 - k} \right)!
\sum\limits_{k + 1 \le n} {{{1^{\,\overline {\,n - 1 - k\,} } }
\over {1^{\,\overline {\,n - 1 - k\,} } \left( {n - k} \right)^{\,\overline {\,N\,} } }}x^{\,n} } } = \cr
& = \sum\limits_{k = 0}^{N - 1} {k!\left( {N - 1 - k} \right)!
\sum\limits_{k + 1 \le n} {{{x^{\,n} } \over {\left( {n - k} \right)^{\,\overline {\,N\,} } }}} } = \cr
& = \left( {N - 1} \right)!\left( {\sum\limits_{k = 0}^{N - 1} {{{x^{\,k} }
\over {\left( \matrix{ N - 1 \cr k \cr} \right)}}} } \right)
\left( {\sum\limits_{1 \le j} {{{x^{\,j} } \over {j^{\,\overline {\,N\,} } }}} } \right) \cr}
$$
where:
I am sure you can continue by yourself.
Let's consider a bit more general case
$$I(x,a)=\int_0^1\cos(xt^a)\tan tdt=\Re\int_0^1e^{ixt^a}\tan tdt\overset{xt^a=z}{=}\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\tan\left(\frac{z}{x}\right)^\frac1az^{\frac1a-1}dz$$
Decomposing $\tan$ into the series
$$=\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\left(\big(\frac{z}{x}\big)^{1/a}+\frac{1}{3}\big(\frac{z}{x}\big)^{3/a}+...\right)z^{\frac 1a-1}dz\tag{1}$$
Let's consider the first term:
$$I_1=\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\Big(\frac{z}{x}\Big)^{1/a}z^{\frac 1a-1}dz=\frac{1}{ax^{2/a}}\Re\int_0^xe^{iz}z^{\frac2a-1}dz\tag{2}$$
Next, we integrate along a quarter-circle in the complex plane: $0\to x\to ix\to i0$, adding the arch of radius $x$. We do not have poles inside our contour; therefore,
$$\oint_C=0\,\Rightarrow\,I_1+\frac{1}{ax^{2/a}}\Re\,e^{\frac{\pi i}{a}}\int_x^0e^{-z}z^{\frac2a-1}dz+I_{C_x}=0$$
where $I_{C_x}$ is the integral along the arch of the radius $x$.
Extending integration of the second term to $\infty$ (and, therefore, dropping exponentially small corrections),
$$I_1\sim\frac{\cos\frac{\pi}{a}\Gamma\left(\frac2a\right)}{ax^{2/a}}-I_{C_x}$$
To estimate $I_{C_x}$, we use the Jordan' lemma. Using $\sin\phi>\frac2\pi\phi$ for $\phi\in[0;\frac\pi2]$,
$$|I_{C_x}|=\frac{1}{ax^{2/a}}\left|\int_0^{\pi/2}e^{ix\cos\phi-x\sin\phi}x^\frac2ae^{i\phi\frac2a}id\phi\right|<\frac{1}{a}\int_0^{\pi/2}e^{-\frac2\pi x\phi}d\phi<\frac{\pi}{2ax}$$
We see that $|I_{C_x}|\ll I_1$ at $x\to\infty$ and can be dropped , if $x^\frac2a\ll x$, or $\,\boxed{a>2}$.
In the same way we can evaluate next terms of $\tan$ decomposition and get, for example,
$$I(x,a)=\int_0^1\cos(xt^a)\tan tdt\sim\frac{\cos\frac{\pi}{a}\Gamma\left(\frac2a\right)}{a\,x^{2/a}}+\frac{\cos\frac{2\pi}{a}\Gamma\left(\frac4a\right)}{3a\,x^{4/a}},\,\,a>4\tag{3}$$
For $a=3$ we are only allowed to keep the first term:
$$\int^1_0 \cos(xt^3)\tan(t)dt\sim\frac{\cos\frac{\pi}{3}\Gamma\left(\frac23\right)}{3x^{2/3}}=\frac16\Gamma\big(\frac23\big)\,x^{-\frac23}\tag{4}$$
$\bf{Explicite \,form \,of\, full \,asymptotic}$
In his post @AxelT evaluated first oscillating asymptotic term, which gives very good approximation even at $x\sim1$. It would be also interesting to get full asymptotic.
Using integration in the complex plane, we got above
$$I(x,a)=\frac{1}{ax^{1/a}}\Re\int_0^xe^{iz}\tan\left(\frac{z}{x}\right)^\frac1az^{\frac1a-1}dz=\frac{1}{ax^{1/a}}\Re \,e^\frac{\pi i}{2a}\int_0^xe^{-z}\tan\left(\frac{z^\frac1a}{x^\frac1a}e^\frac{\pi i}{2a}\right)z^{\frac1a-1}dz$$
$$-\,I_{C_x}\tag{5}$$
where $I_{C_x}$ is the integral along the arch of the radius $x$:
$$I_{C_x}=\frac{1}{ax^{1/a}}\Re\int_0^{\pi/2}e^{ixe^{i\phi}}e^{\frac{i\phi}{a}}\tan\big(e^{\frac{i\phi}{a}}\big)x^\frac1ae^{i\phi}id\phi\overset{s=e^{i\phi}}{=}\frac1a\int_1^{e^{\frac{\pi i}{2}}}e^{ixs}\tan(s^\frac1a)s^\frac1ads\tag{6}$$
Integrating (6) several times by part and dropping exponentially small terms ($\sim e^{-\frac{\pi x}{2}}$)
$$I_{C_x}\sim\frac{1}{a}\Re\left(-\frac{e^{ix}}{ix}f^{(0)}(s,a)|_{s=1}+\frac{e^{ix}}{(ix)^2}f^{(1)}(s,a)|_{s=1}-\frac{e^{ix}}{(ix)^3}f^{(2)}(s,a)|_{s=1}-...\right)$$
where we denoted
$$f^{(k)}(s,a)|_{s=1}:=\frac{\partial^k}{\partial s^k}s^\frac1a\tan(s^\frac1a)\Big|_{s=1}$$
Rearranging
$$-I_{C_x}\sim\frac{\sin x}{ax}\left(f^{(0)}(s,a)-\frac{f^{(2)}(s,a)}{x^2}+...\right)\bigg|_{s=1}+\frac{\cos x}{ax^2}\left(f^{(1)}(s,a)-\frac{f^{(3)}(s,a)}{x^2}+...\right)\bigg|_{s=1}$$
$$=\frac1a\sum_{n=1}^\infty\frac{\sin\big(\frac{\pi(n-1)}{2}+x\big)}{x^n}\left(s^\frac1a\tan s^\frac1a\right)^{(n-1)}\,\bigg|_{s=1}\tag{7}$$
Using the series of $\tan$
$$\tan t=\sum_{n=1}^\infty(-1)^{n-1}\frac{4^n(4^n-1)}{(2n)!}B_{2n}t^{2n-1}$$
and integrating term by term (extending integration to $\infty$)
$$\frac{1}{ax^{1/a}}\Re \,e^\frac{\pi i}{2a}\int_0^xe^{-z}\tan\left(\frac{z^\frac1a}{x^\frac1a}e^\frac{\pi i}{2a}\right)z^{\frac1a-1}dz\sim\frac1a\sum_{n=1}^\infty\frac{(-1)^{n-1}4^n(4^n-1)B_{2n}}{(2n)!}\frac{\cos\frac{\pi n}{a}\Gamma\big(\frac{2n}{a}\big)}{x^\frac{2n}{a}}\tag{8}$$
Putting (7) and (8) into (5), we get the desired asymptotics:
$$\color{blue}{I(x,a)\sim\frac1a\sum_{n=1}^\infty\frac{(-1)^{n-1}4^n(4^n-1)B_{2n}}{(2n)!}\frac{\cos\frac{\pi n}{a}\Gamma\big(\frac{2n}{a}\big)}{x^\frac{2n}{a}}}$$
$$\color{blue}{+\,\frac1a\sum_{n=1}^\infty\frac{\sin\big(\frac{\pi(n-1)}{2}+x\big)}{x^n}\left(s^\frac1a\tan s^\frac1a\right)^{(n-1)}\,\bigg|_{s=1}}$$
Here we have the only limitation $\color{blue}{a>0}$.
If we want to take several first terms, the solution - the sequence of the terms we have to keep - depends, of course, on $a$. For example, for $a=3$ we get the series
$$I(x,3)\sim\frac{\cos\frac\pi3\Gamma\big(\frac23\big)}{3x^\frac23}+\frac{\sin x\tan1}{3x}+\frac{\cos\frac{2\pi}{3}\Gamma\big(\frac43\big)}{9x^\frac43}+\frac{2\cos\pi\,\Gamma(2)}{45x^2}+\frac{\cos x}{9x^2}\left(\tan1+\frac{1}{\cos^21}\right)$$
$$+\,O\left(\frac1{x^\frac83}\right)$$
Best Answer
As noticed, a simple Laplace method cannot be used here, as 2 scales are involved. A uniform asymptotic expansion should be found. Alternatively, in this case, we can recognize a modified Bessel function. Indeed, changing $x=\frac{h}{M}\sinh t$, the integral can be written as \begin{align} I&=\int_{-\infty}^\infty \exp(-\sqrt{h^2+M^2x^2}) \,dx\\ &=\frac{h}{M}\int_{-\infty}^\infty \exp(-h\cosh t)\cosh t \,dt \end{align} which is proportional to an integral representation of a modified Bessel function (DLMF): \begin{equation} K_{\nu}\left(z\right)=\int_{0}^{\infty}e^{-z\cosh t}\cosh\left(\nu t\right)% \mathrm{d}t \end{equation} with $\nu=1,z=h$, \begin{equation} I=\frac{2h}{M}K_1(h) \end{equation} Using the series expansion of the Bessel function near $h=0$ (DLMF), \begin{align} I&\sim \frac{2h}{M}\left[ \frac{1}{h}+\frac{h}{2}\left(\ln\frac{h}{2} +\gamma-\frac{1}{2}\right)+\ldots\right]\\ &\sim\frac{2}{M}+\frac{h^2}{M}\left(\ln\frac{h}{2} +\gamma-\frac{1}{2}\right)+\ldots \end{align} EDIT:Another method using the Mellin transform technique
Changing $x=th/M$ and then $u=\sqrt{1+t^2}-1$, the problem is equivalent to find the small $h$ behavior of \begin{align} I&=2\frac{h}{M}\int_0^\infty\exp(-h\sqrt{1+t^2})\,dt\\ &=2\frac{h}{M}e^{-h}\int_0^\infty\frac{u+1}{\sqrt{u(u+2)}}\exp(-hu)\,du \end{align} We have thus to find a Laplace transforms with a small parameter. A classical method which uses the Mellin transform technique is given in (DLMF). Intermediate results are given below, with the help of a CAS. Defining \begin{equation} H(u)=\frac{u+1}{\sqrt{u(u+2)}} \end{equation} he following behaviors hold:\ \begin{array}{lll} H(u)&\sim 1+\frac{1}{2u^2}+O\left( u^{-3} \right)& \text{ for }u\to\infty\\ &\sim O(u^{-1/2})& \text{ for } u\to 0 \end{array} and the Mellin transform is \begin{equation} \mathcal{M}\left[H(u) \right](z)=\pi^{-1/2}2^{z-1}(z-1)\Gamma(-z)\Gamma\left( z-\frac{1}{2} \right) \end{equation} For the function \begin{equation} F(z)=-h^{z-1}\Gamma(1-z)\mathcal{M}\left[H(u) \right](z) \end{equation} the residues for $z=0,1,2,3$ are \begin{align} \left. \operatorname{res}F(z)\right|_{z=0}&=\frac{1}{h}\\ \left. \operatorname{res}F(z)\right|_{z=1}&=1\\ \left. \operatorname{res}F(z)\right|_{z=2}&=\frac{h}{2}\left[ \ln\frac{h}{2}+\gamma+\frac{1}{2}\right]\\ \left. \operatorname{res}F(z)\right|_{z=3}&=\frac{h^2}{12}\left[6 \ln\frac{h}{2}+6\gamma-1\right] \end{align} With $e^{-h}= 1-h+h^2/2+O(h^3)$, \begin{align} I&\sim \frac{2h}{M}(1-h+\frac{h^2}{2})\left[ \frac{1}{h}+\frac{h}{2}\left( \ln\frac{h}{2}+\gamma+\frac{1}{2}\right)+\frac{h^2}{12}\left(6 \ln\frac{h}{2}+6\gamma-1\right)\right]\\ &\sim \frac{2}{M}+\frac{h^2}{M}\left( \ln\frac{h}{2}+\gamma-\frac{1}{2} \right) \end{align}