The results in the referenced paper do not directly give the generating function of $\operatorname{Tr} [(D+uu^T)^k]$, but the underlying technique can be used to obtain it.
Actually, it is convenient to start with a little bit more general quantity - the (matrix-valued) generating function of $A_k\equiv(D+uu^T)^k$. First, let's expand $(D+uu^T)^k$ using binomial formula, where in each term we can have $0\leq m\leq k$ choices of $uu^T$ positioned at $1\leq t_1 < t_2 < \dots < t_m \leq k$
$$
\begin{align}
A_k &= \sum_{m=0}^k \; \sum_{1\leq t_1 < \dots < t_m \leq k} D^{t_1-1}uu^TD^{t_2-t_1-1}uu^T \ldots uu^T D^{k-t_m-1} \\
\tag{1}
\label{binom}
&= D^k + \sum_{1 \leq t_1 \leq k} D^{t_1-1}uu^T A_{k-t_1}
\end{align}
$$
Here in the second equality we separated the term with $m=0$ and all the other terms with $m\geq1$. Then, we used that all terms with fixed $t_1$ share the same first $t_1$ factors $D^{t_1-1}uu^T$, and all the configurations of the rest $k-t_1$ factors are exactly the same as in $A_{k-t_1}$.
Next, let's denote generating functions of $A_k$ and $D^k$ as $\widetilde{A}(z)=\sum_{k=0}^\infty z^k A_k$ and $\widetilde{D}(z)=\sum_{k=0}^\infty z^k D^k$. Now, thanks to convolution nature of the sum $\sum_{1 \leq t_1 \leq k} D^{t_1-1}uu^T A_{k-t_1}$, we can use our representation \eqref{binom} of $A_k$ to get a algebraic equation for $\widetilde{A}(z)$
$$
\begin{align}
\widetilde{A}(z)&=\widetilde{D}(z)+\sum_{k=0}^\infty z^k\sum_{1\leq t_1 \leq k} D^{t_1-1}uu^TA_{k-t_1} \\
&=\widetilde{D}(z) + z \left(\sum_{t_1=1}^\infty z^{t_1-1}D^{t_1-1}\right) uu^T\left(\sum_{k=t_1}^{\infty}z^{k-t_1}A_{k-t_1}\right) \\
\tag{2}
\label{Atilde_eq}
&= \widetilde{D}(z)+z\widetilde{D}(z) uu^T \widetilde{A}(z)
\end{align}
$$
Finally, we use Woodberry matrix inversion formula to solve \eqref{Atilde_eq} w.r.t. $\widetilde{A}(z)$
$$
\tag{3}
\label{Atilde_solution}
\widetilde{A}(z) = \left(I - z \widetilde{D}(z) uu^T \right)^{-1} \widetilde{D}(z)=\widetilde{D}(z)+\frac{z}{1-z u^T \widetilde{D}(z) u}\widetilde{D}(z)uu^T\widetilde{D}(z)
$$
From this point, it is easy to use \eqref{Atilde_solution} to get simpler results in a particular application. For example, in paper authors are interested in $u^T (D+uu^T)^k v$ with generating function $u^T \widetilde{A}(z)v$ for a certain vector $v$. Multiplying \eqref{Atilde_solution} by $u^T$ and $v$, this scalar generating function reduces to two other scalar generating functions $u^T \widetilde{D}(z) u$ and $u^T \widetilde{D}(z) v$ as
$$
u^T \widetilde{A}(z)v = \frac{u^T \widetilde{D}(z)v}{1-zu^T \widetilde{D}(z)u}
$$
Going back to the original question of finding generating function $\operatorname{Tr}[\widetilde{A}(z)]$ of $\operatorname{Tr} [(D+uu^T)^k]$, we just need to take trace of \eqref{Atilde_solution}. Then, the result is expressed in terms of 3 other scalar generating functions $\operatorname{Tr}[\widetilde{D}(z)]$, $u^T \widetilde{D}(z)u$ and $u^T \big(\widetilde{D}(z)\big)^2u$ as
$$
\operatorname{Tr}[\widetilde{A}(z)] = \operatorname{Tr}[\widetilde{D}(z)] + \frac{z}{1-zu^T\widetilde{D}(z)u} u^T \big(\widetilde{D}(z)\big)^2 u
$$
Now, the last question of finding the large $k$ asymptotic of $\operatorname{Tr} [(D+uu^T)^k]$ can be answered, in certain cases, using generating function. For example, iт the case of a power-law asymptotic $\operatorname{Tr} [(D+uu^T)^k]\sim k^{-s}$, we can obtain exponent $s$, and also constant if needed, by looking at the singularity of $\operatorname{Tr} [\widetilde{A}(z)]$ at $z=1$. This can be done in the following way. Assume a sequence $a_k$ has asymptotic $a_k\sim k^{-s}$ and $\widetilde{a}(z)$ is its generating function. Then, if $0<s<1$ the sequence of partial sums diverges as $\sum_{l=0}^k a_l \sim k^{1-s}$, which implies the singularity of generating function $\widetilde{a}(1-\varepsilon)\sim \varepsilon^{s-1}$. In our problem, we go in the opposite direction and find the singularity of generating function to infer the asymptotic of the sequence. If $1\leq s<2$, the powe-law singularity in generating function will disappear but persist in its first derivative $\widetilde{a}'(z)=\sum_k ka_k z^{k-1}$ as $\widetilde{a}'(1-\varepsilon)\sim \varepsilon^{s-2}$, which gives $\sum_{l=0}^k l a_l \sim k^{2-s}$. And for $s\geq 2$ we can do the same by taking more derivatives.
