The constant term is difficult to construct. I'll present the terms to $O(1/n)$ which includes the constant. This problem was asked in MSE 2891159, in which case I answered one of three; this problem is that problem's $\sigma_c(n).$ The problem can be setup exactly as that previous problem, by differentiating with respect to $s$ the following, and taking limits when needed:
$$ \sum_{k=1}^n k^{-s}=\zeta{(s)}-\frac{1}{(s-1)n^{s-1}}+\frac{1}{2n^s}-
\frac{1}{12}\frac{s}{n^{s+1}} +.... $$
Higher terms are neglected because what is present is sufficient to get to $O(1/n)$ terms. The formula follows from the Euler-McLaurin summation. As in the other problem, break the sum you want into primitives consisting of sums over $\log{k}/k^m$ or $\log^2{k}/k^m.$ That is,
$$\sum_{k=1}^n H_k \log{k} = \sum_{k=1}^n \Big(\gamma + \log{k} + \frac{1}{2k} -\frac{1}{12k^2} \Big) \log{k}$$
where the asymptotic formula for $H_k = \gamma + \psi(k+1)$ has been used and a sufficient number of terms have been taken, $except$ for those need to derive the constant term (we'll come back to that later). Letting $L=\log{n}$ we have the following:
$$ v_0=\sum_{k=1}^n \log^2{k}=n\big(L^2+2(1-L)) + \frac{L^2}{2}+\frac{L}{6n}+\frac{\gamma^2}{2}-\frac{\pi^2}{24}-\frac{1}{2}\log^2{(2\pi)}+\gamma_1$$
$$ \quad v_1=\sum_{k=1}^n \log{k}=n(L-1)+\frac{L}{2}-\zeta'(0)+\frac{1}{12n}$$
$$v_2= \sum_{k=1}^n \frac{\log{k}}{k} =\frac{L^2}{2}+\gamma_1+\frac{L}{2n}$$
$$ v_3=\sum_{k=1}^n \frac{\log{k}}{k^2} =-\zeta'(2)-\frac{L+1}{n}$$
Let $\tilde{v_k}=(v_k$ with constant term set to 0). Then
$$\sigma_c(n)=\sum_{k=1}^n H_k \log{k} =\tilde{v_0}+\gamma\, \tilde{v_1} + \frac{1}{2} \tilde{v_2}-\frac{1}{12} \tilde{v_3} + C $$
$$=n\big(L^2+(\gamma-2)(L-1)\big) + \frac{3L^2}{4}+\gamma\frac{L}{2}+\frac{1}{2n}(L+\frac{\gamma+1}{6}) + C $$
where $C$ is the unknown constant term. A correct way to determine C is simply by
$$ C=\lim_{n \to \infty} \Big(\sum_{k=1}^n H_k \log{k} - \Big(n\big(L^2+(\gamma-2)(L-1)\big) + \frac{3L^2}{4}+\gamma\frac{L}{2} \Big)\, \Big)$$
Proceed further only if you are comfortable with formal mathematics. Now I'm going to find $C$ in a different manner. It is well-known that
$$ \gamma = \lim_{n \to \infty}\Big( \sum_{k=1}^n\frac{1}{k} - \log{n} \Big) =
\int_0^\infty \Big(\frac{1}{e^x-1} - \frac{e^{-x}}{x} \Big)dx.$$ The first expression is similar in spirit to how we defined $C,$ but the second equation is useful in that many digits of $\gamma$ can be extracted quite easily with your favorite numerical integration routines.
In derving $v_0$ through $v_3$ I stopped because this is sufficient to get terms to $O(1/n).$ If we kept going, there would be a $v_5 = \sum_{k=1}^n \log{k}/k^4= -\zeta'(4) + o(1/n)$, etc. The natural thing to do is add all the $-\zeta'(2n)$ terms up with the appropriate weights to get a new constant that, when added to the ones in the $v_k$, constitute $C.$ This was done in MSE 2891159 and very fortunately the series converged. No such luck here. The weights are from the asymptotic expansion of the harmonic numbers,
$$H_k = \gamma + \log{k} + \frac{1}{2k}-\sum_{m=1}^\infty \frac{B_{2m}}{2m} k^{-2m} $$ so the constant we want to assign meaning to is
$$ \kappa \,\,\dot{=} \sum_{m=1}^\infty \frac{B_{2m}}{2m} \zeta'(2m) $$
where the dotted equals means 'representation' instead of 'equals.' The idea is akin to saying $$\sum_{k=1}^\infty k \,\, \dot{=} -1/12$$ which can be given rigor in terms of zeta regularization. First differentiate the well-known Euler-like integral for the zeta function to find
$$\zeta'(2s) =\frac{1}{\Gamma(2s)} \int_0^\infty \frac{t^{2s-1}}{e^t-1} \Big( \log{t} - \psi(2s) \Big) .$$
Insert into the definition for $\kappa$ and interchange $\int$ and $\sum$
$$ (K) \quad \kappa \,\,\dot{=} \int_0^\infty \frac{dt/t}{e^t-1}\sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \big( \log{t} - \psi(2m) \big). $$ The inner sum converges (its been Borel transformed), we just need to find an expression that's not a power series. We need some formulas. It is well-known that, and it will be used repeatedly,
$$\sum_{m=1}^\infty \frac{B_m}{m!} t^m = \frac{t}{e^t-1} -1, \quad \sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m}= \frac{t}{e^t-1} -1+t/2 $$
An integration of the previous equation leads to
$$\sum_{m=1}^\infty\frac{B_{2m}}{(2m)(2m)!} t^{2m}= -\log\big(\frac{t}{e^t-1}\big) -t/2 $$ Let's begin on the 'psi' term:
$$ \Psi(t):=\sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \,\psi(2m)=
\sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \,\big(\psi(2m)+\frac{1}{2m} - \frac{1}{2m} \big)=$$
$$= \sum_{m=1}^\infty \frac{B_{2m}}{(2m)!} t^{2m} \,\psi(2m+1)+\Big(t/2+\log{\big(\frac{t}{e^t-1}\big)} \Big) =$$
$$= \sum_{m=1}^\infty \frac{B_{m}}{m!} t^m \,\psi(m+1)+
\Big( \frac{1}{2}(1-\gamma)t \Big) +
\Big(t/2+\log{\big(\frac{t}{e^t-1}\big)} \Big) .$$
Again use $H_m = \gamma + \psi(m+1)$ to find
$$\Psi(t)=\sum_{m=1}^\infty \frac{B_m}{m!} t^m \,H_m -
\gamma\Big(\frac{t}{e^t-1} - 1\Big)+ \big(1-\frac{\gamma}{2}\big)t + \log{\big(\frac{t}{e^t-1}\big) }$$
The harmonic number is used because there is the integral relation
$$H_m = -m \int_0^1 dx x^{m-1} \log{(1-x)} $$ Insert this, switch $\int$ and $\sum$, sum up the series in closed form, and $finally$
$$\Psi(t)=-\int_0^t \frac{\log{(1-u/t)}}{e^u-1}\Big(1-\frac{u\,e^u}{e^u-1}\Big)\,du-
\gamma\Big(\frac{t}{e^t-1} - 1\Big)+ \big(1-\frac{\gamma}{2}\big)t + \log{\big(\frac{t}{e^t-1}\big) }$$
In eq. $(K)$ the sum before the $\log{t}$ is one of the known formulas. Therefore a have a (double) integral relationship that is perfectly well-behave for $\kappa$
$$\kappa=\int_0^\infty \frac{dt/t}{e^t-1}\Big(\, \log{t}\big(\frac{t}{e^t-1}-1+t/2\big) - \Psi(t)\Big) = -0.077596...$$
When using $\kappa$ the answer can be stated as
$$ C=\kappa + \frac{\gamma^2}{2}-\frac{\pi^2}{24}+\frac{1}{2}\gamma \, \log(2\pi) - \frac{1}{2} \log^2{(2\pi)} + \frac{3}{2} \gamma_1 $$
Six digits of agreement are obtained by using the asymptotic formula and comparing to the brute force summation for $n=$ 20.
Rewriting your sum as :
\begin{align}
\tag{1}S(n)&:=\frac {n!}{n^n}\sum_{k=1}^{n} \frac{n^{n-k}}{k\,(n-k)!}\\
&=\frac {n!}{n^n}\sum_{i=0}^{n-1} \frac{n^{i}}{(n-i)\,i!}\\
&=\frac {n!}{n^n}\sum_{i=0}^{n-1} \frac{i^{i}}{(n-i)\,i!}\left(\frac {n}{i}\right)^i\\
&\approx\sqrt{2\pi\,n}\sum_{i=1}^{n-1} \frac{e^{\,i-n}}{(n-i)\,\sqrt{2\pi i}}\left(\frac {n}{i}\right)^i\\
\tag{2}&\approx\sum_{i=1}^{n-1} \frac{e^{\,i-n}}{(n-i)}\left(\frac {n}{i}\right)^{i+1/2}\,=:A(n)\\
\end{align}
Stirling's approximation with the error term $\displaystyle\; n!\approx\sqrt{2 \pi\,n}\, \left(\frac{n}{e} \right)^n \left( 1 + O \left( \frac{1}{n} \right) \right)$ was used.
Observe that the last term of $A(n)$ is nearly $\dfrac 11$, the previous nearly $\dfrac 12$ and then $\dfrac 13$ and so on so that we obtain 'nearly' a sum of terms $\dfrac 1k$ with an error coefficient of order $\left(1+ O \left( \frac{1}{n-k} \right)\right)\left(1+ O \left( \frac{1}{n} \right)\right)$.
This allows to obtain $\;S(n)-A(n)=O \left( \frac{\log n}{n} \right)$ and thus to study $A(n)$ instead of $S(n)$.
We can't conclude that $\,A(n)\sim H(n)\,$ because the values of the sum is decreasing much faster than $\dfrac 1k$ after "a while". More exactly for $\,k:=n-i\,$ and $\,n\gg 1\,$ the first $k$ terms become :
\begin{align}
b_k&:= \frac{e^{\,-k}}{k}\left(\frac {n}{n-k}\right)^{n-k+1/2}\\
&= \frac{e^{\,-k}}{k}\left(1-\frac kn\right)^{\large{-n\left(1-\frac{k-1/2}n\right)}},\quad\text{expanding }\;\left(1-\frac kn\right)^{-n}\;\text{to second order :}\\
&\approx \frac{e^{\,-k}}{k}e^{\large{\left(k+\frac{k^2}{2n}\right)\left(1-\frac{k-1/2}n\right)}},\quad\text{neglecting the $\frac 1{n^2}$ terms after simplification :}\\
&\approx \frac{e^{\large-{\frac{k\,(k-1)}{2\,n}}}}{k}\\
\end{align}
Notice that the exponential will become smaller than $\dfrac 12$ for $\,k\,$ somewhere between $\,\left[\sqrt{n}\right]\,$ and $\,\left[\sqrt{2\,n}\right]\,$ so that the main term of the asymptotic should rather be $\,\dfrac {\log(n)}2\,$ than $\,\log(n)$.
With some black magic used to get the additional corrective term we will write our new approximation as a series (since $b_k$ is very small for $\,k\ge n\gg 1$) :
$$\tag{3}B(n):=\sum_{k=1}^{\infty} \frac{e^{\large-{\frac{k\,(k-1)}{2\,n}}}}{k}-\frac 16\sqrt{\frac{\pi}{2\,n}}$$
The interest of the approximations $(2)$ and $(3)$ is that numerical evaluation is much faster with these two sums (as you may see in the table at the end). Further $(3)$ allows to get the asymptotic term with any precision (I got $\dfrac{\gamma+\log(2)}2$) and thus my "final conjecture" with $\gamma$ the Euler–Mascheroni constant :
$$\tag{4} S(n)= \frac{\log(2\,n)+\gamma}2+O\left(\frac 1{\sqrt{n}}\right)$$
Table of values for $\ \displaystyle S(n)-\frac 12 \log(n)\ $ and the corresponding approximations $(2)$ and $(3)$ :
\begin{array} {c|c|c|c}
n&\text{Exact}\;S(n)-\log(n)/2&A(n)-\log(n)/2&B(n)-\log(n)/2\\
\hline
1000& 0.648336893251805& 0.648340255305227& 0.648310001946645\\
10000& 0.639353578731186& 0.639353683741043& 0.639350829257599\\
100000& 0.636501976201853 & 0.636501979510217& 0.636501699321884\\
1000000&-& 0.635599138656474& 0.635599110802617\\
10000000&-& 0.635313528088334& 0.635313525308151\\
100000000&-& 0.635223199313196& 0.635223199035342\\
1000000000&-& 0.635194633766179& 0.635194633738398\\
\infty&-&-& 0.635181422730739\\
\end{array}
(the limits as $\,n\to \infty\,$ should all be equal!)
Best Answer
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.