Finding a closed form for a limit of a sequence Consider a sequence $$u_n=\sum_{r=1}^{n} \frac{{n\choose r} f^{(r)}(1)}{(r-1)!} $$ where $f$ is an infinitely differentiable real valued function on $\mathbb
{R}$.

Question: Given that the limit exists, find a closed form for $$\lim_{n\to \infty} u_n$$

My try: Using binomial theorem $$(x+y)^n= x^n+{n\choose1}x^{n-1}y+{n\choose2}x^{n-2}y^2+{n\choose3}x^{n-3}y^3+\dots+{n\choose r}x^{n-r}y^r +\dots+y^n $$
I cannot see how to get a closed form for $\lim_{n\to \infty} u_n$.
Any help will be appreciated. Thank you.
 A: New Answer. Define $h(x) = f\left(\frac{1}{1-x}\right)$. Then we have:

Claim. $\displaystyle u_n := \sum_{k=1}^{n} \binom{n}{k} \frac{f^{(k)}(1)}{(k-1)!} = \frac{h^{(n)}(0)}{(n-1)!} $.
(The proof is postponed to the end.)

Now, assume that $f(z)$ is analytic for $\operatorname{Re}(z) > \frac{1}{2}$ (so that $h(z)$ is analytic for $|z| < 1$) and that either (i) $(u_n)$ converges in $\mathbb{C}$, or (ii) $(u_n) \subseteq \mathbb{R} $ and it has a limit in $[-\infty, \infty]$. Then the limit of $(u_n)$ is also the Abel mean of $(u_n)$:
$$ \bbox[color:navy;padding:8px;border:1px navy solid;]{\lim_{n\to\infty} u_n
= \lim_{r \to 1^-} \frac{\sum_{n=1}^{\infty} u_n r^n}{\sum_{n=0}^{\infty} r^n}
= \lim_{r \to 1^-} \frac{h'(r)}{1/(1-r)}
= \lim_{x\to\infty} \frac{f'(x)}{x}.} \tag{$\diamond$} $$
For example, if $f(x) = \log(1+\log x)$ and assuming that the corresponding $(u_n)$ converges, then the limit must coincide with
$$ \lim_{x\to\infty} \frac{f'(x)}{x} = \lim_{x\to\infty} \frac{1}{x^2(1+\log x)} = 0. $$
However, in the case $f(s) = \xi(s) = \frac{1}{2} s(s-1) \pi^{-s/2} \Gamma\left(\frac {s}{2}\right) \zeta(s)$ is the Riemann Xi function, the limit $(\diamond)$ diverges to $+\infty$ since $\Gamma(s/2)$ grows super-exponentially fast as $s \to \infty$. Therefore $(u_n)$ does not converge as well. This can also be glimpsed from the exact computation:

Proof of Claim. When $f$ is analytic near $1$, then this formula can be easily proved by noting that
\begin{align*}
u_n
&= \sum_{k=1}^{n} \binom{n}{k} \left[ \frac{1}{2\pi i}  \int_{|z-1| = \varepsilon} \frac{f'(z)}{(z-1)^k} \, \mathrm{d}z \right] \\
&= \frac{1}{2\pi i} \int_{|z-1| = \varepsilon} \left(\frac{z}{z-1}\right)^n f'(z) \, \mathrm{d}z \\
&= \frac{1}{2\pi i} \int_{|w| = \varepsilon'} \frac{h'(w)}{w^n} \mathrm{d}w. \tag{$w = \tfrac{z-1}{z}$}
\end{align*}
In general, if $f$ is only assumed to be indefinitely differentiable near $1$, write $p(x) = \frac{1}{1-x}$ and note that $p^{(k)}(0) = k!$. Then by the Faà di Bruno's formula,
\begin{align*}
\frac{h^{(n)}(0)}{(n-1)!}
&= \frac{1}{(n-1)!} (f \circ p)^{(n)}(0) \\
&= \frac{n!}{(n-1)!} \sum_{k=1}^{n} f^{(k)}(p(0)) \sum_{\substack{a_1 + 2a_2 + \cdots + na_n = n \\ a_1 + a_2 + \cdots + a_n = k}} \prod_{i=1}^{n} \frac{1}{a_i!} \biggl( \frac{p^{(i)}(0)}{i!} \biggr)^{a_i} \\
&= n \sum_{k=1}^{n} f^{(k)}(1) \sum_{\substack{a_1 + 2a_2 + \cdots + na_n = n \\ a_1 + a_2 + \cdots + a_n = k}} \prod_{i=1}^{n} \frac{1}{a_i!}
\end{align*}
However, for each given $(a_1, \ldots, a_n)$ with $\sum_{i=1}^{n} ia_i = n$ and $\sum_{i=1}^{n} a_i = k$,
$$ k! \prod_{i=1}^{n} \frac{1}{a_i!} = \binom{k}{a_1, \ldots, a_n} $$
is the number of ways of arranging $a_1$ $\bullet$'s, $a_2$ $\bullet\bullet$'s, ..., $a_n$ $\overbrace{\bullet\bullet\cdots\bullet}^{n}$'s. So, the sum of this multinomial coefficient for all such $(a_i)$'s gives the number of ways of partitioning a chain of $n$ $\bullet$'s into $k$ parts, each having at least one $\bullet$. This is precisely $\binom{n-1}{k-1}$, and so, it follows that
\begin{align*}
\frac{h^{(n)}(0)}{(n-1)!}
= n \sum_{k=1}^{n} f^{(k)}(1) \frac{1}{k!} \binom{n-1}{k-1}
= \sum_{k=1}^{n} \binom{n}{k} \frac{f^{(k)}(p(0))}{(k-1)!}
= u_n,
\end{align*}
proving the desired claim.

Old Answer. As OP mentioned in the comment, we will work with $f(x)=\log(1+\log x)$. Then for any small $\varepsilon > 0$, Cauchy integration formula and the binomial theorem together yield
\begin{align*}
u_n
= \sum_{k=1}^{n} \binom{n}{k} \left[ \frac{1}{2\pi i}  \int_{|z-1| = \varepsilon} \frac{f'(z)}{(z-1)^k} \, \mathrm{d}z \right]
= \frac{1}{2\pi i} \int_{|z-1| = \varepsilon} g(z) \, \mathrm{d}z,
\end{align*}
where
$$ g(z) = \left(\frac{z}{z-1}\right)^n \frac{1}{z(1+\log z)}. $$
Now, "inflating" the contour $|z-1| = \varepsilon$ continuously while wrapping any singularities of $g(z)$ encountered during this process, and noting that the integrand is $\mathcal{O}\left(\frac{1}{|z|\log|z|}\right)$ as $|z| \to \infty$, it follows that
\begin{align*}
u_n
&= - \mathop{\mathrm{Res}}_{z=e^{-1}} g(z) + \frac{1}{2\pi i} \int_{-\infty+0^+i}^{0^+i} g(z) \, \mathrm{d}z - \frac{1}{2\pi i} \int_{-\infty+0^-i}^{0^-i} g(z) \, \mathrm{d}z \\
&= \frac{(-1)^{n-1}}{(e - 1)^n} + \int_{0}^{\infty} \left( \frac{x}{x+1} \right)^n \frac{1}{x \left( (1 + \log x)^2 + \pi^2 \right)} \, \mathrm{d}x \\
&= \frac{(-1)^{n-1}}{(e - 1)^n} + \int_{-\infty}^{\infty} \left( \frac{e^t}{e^t + 1}\right)^n \frac{1}{(1 + t)^2 + \pi^2} \, \mathrm{d}t. \tag{*}
\end{align*}
For example, the picture below compares the exact values of $u_n$ computed by the original formula and the numerical values of $u_n$ obtained from the integral representation above:

Now, since the integrand of $\text{(*)}$ is bounded by the integrable function $\frac{1}{(t+1)^2 + \pi^2}$, we can apply the dominated convergence theorem to get
$$ \lim_{n\to\infty} u_n
= \int_{-\infty}^{\infty} \lim_{n\to\infty} \left( \frac{e^t}{e^t + 1}\right)^n \frac{1}{(1 + t)^2 + \pi^2} \, \mathrm{d}t
= 0. $$
