Let me point out that the cited paper is a landmark brilliant event
containing profound insights that surpass the MSE question answer
format.
As an example we do the calculation of the expected cycle length
$\mu(u_0)$ and the tail length $\lambda(u_0)$ using the labelled tree
function. We refer the reader to the paper as regards the
asymptotics.
We will provide a closed form of the exponential generating function
of the quantities that are involved.
The species of labelled trees has the specification
$$\def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}\mathcal{T} =
\mathcal{Z} \times \textsc{SET}(\mathcal{T})$$
which gives the functional equation
$$T(z) = z \exp T(z).$$
We have that $$T(z) = \sum_{n\ge 1} n^{n-1} \frac{z^n}{n!}.$$
Cycle length.
To compute the expected cycle length note that these random mappings
are sets of cycles of trees having combinatorial specification
$$\textsc{SET}
\left(\sum_{q\ge 1} \textsc{CYC}_{=q}(\mathcal{T}(\mathcal{Z}))\right).$$
Now a tree on $n$ nodes that goes into a cycle of size $q$ contributes
$nq$ to the total count so we get the species
$$\textsc{SET}
\left(\sum_{q\ge 1}
\textsc{CYC}_{=q}(\mathcal{T}(\mathcal{V}^q\mathcal{Z}))\right).$$
It follows that the bivariate generating function of mappings by node
count and cycle size is
$$G(z, v) = \exp\left(\sum_{q\ge 1} \frac{T(v^q z)^q}{q}\right).$$
This yields for the generating function $H(z)$ of the expected cycle size
$$H(z) = \left.\frac{\partial}{\partial v} G(z,v)\right|_{v=1}.$$
This works out to
$$H(z) = \left. \exp\left(\sum_{q\ge 1} \frac{T(v^q z)^q}{q}\right)
\left(\sum_{q\ge 1}
\frac{qT(v^q z)^{q-1} T'(v^q z)q v^{q-1} z}{q}\right)
\right|_{v=1}
\\ = \exp\left(\sum_{q\ge 1} \frac{T(z)^q}{q}\right)
\left(\sum_{q\ge 1}
T(z)^{q-1} T'(z)q z\right)
\\ = \frac{z T'(z)}{(1-T(z))^2} \exp\log\frac{1}{1-T(z)}
\\ = \frac{z T'(z)}{(1-T(z))^3}.$$
We get from the functional equation
$$T'(z) = \exp(T(z)) + z \exp(T(z)) T'(z)$$
so that
$$T'(z) = \frac{\exp(T(z))}{1-z\exp(T(z))}
= \frac{T(z)/z}{1-zT(z)/z}
= \frac{1}{z} \frac{T(z)}{1-T(z)}$$
and therefore
$$H(z) = \frac{T(z)}{(1-T(z))^4}.$$
Extracting coefficients via Lagrange inversion we have
$$Q_n = n! [z^n] H(z) =
n! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}}
\frac{T(z)}{(1-T(z))^4} dz.$$
Put $T(z)=w$ so that $z=w/\exp(w) = w\exp(-w)$ and
$dz = \exp(-w) - w\exp(-w)$
to get
$$n! \frac{1}{2\pi i}
\int_{|w|=\gamma}
\frac{\exp(w(n+1))}{w^{n+1}}
\frac{w}{(1-w)^4}
(\exp(-w) - w\exp(-w)) \; dw
\\ = n! \frac{1}{2\pi i}
\int_{|w|=\gamma}
\frac{\exp(wn)}{w^{n+1}}
\frac{w}{(1-w)^4} (1 - w) \; dw
\\ = n! \frac{1}{2\pi i}
\int_{|w|=\gamma}
\frac{\exp(wn)}{w^{n}}
\frac{1}{(1-w)^3} \; dw.$$
This gives the closed form
$$Q_n = \frac{1}{2}
n! \sum_{q=0}^{n-1} \frac{n^q}{q!}
(n+1-q) (n-q).$$
The average in question is
$$\frac{Q_n}{n^n \times n}.$$
The sequence $Q_n$ starts with
$$1, 10, 117, 1648, 27425, 528336, 11581885,
\\ 284878336, 7772592897, 233010784000,\ldots$$
The Maple code for verification by total enumeration of this sequence
was as follows.
Q :=
proc(n)
option remember;
local ind, d, gf, pos, q, x, seen, traj, cycinit;
if n = 1 then return v fi;
gf := 0;
for ind from n^n to 2*n^n-1 do
d := convert(ind, base, n);
d := map(l->l+1, [seq(d[q], q=1..n)]);
q := 0;
for pos to n do
seen := {}; x := pos; traj := [];
while not(x in seen) do
traj := [op(traj), x];
seen := seen union {x};
x := d[x];
od;
cycinit := 1;
while traj[cycinit] <> x do
cycinit := cycinit + 1;
od;
q := q + nops(traj)-cycinit+1;
od;
gf := gf+v^q;
od;
gf;
end;
EX :=
proc(n)
local T;
T := solve(TT=z*exp(TT), TT);
n!*coeftayl(T/(1-T)^4, z=0, n);
end;
EX2 :=
proc(n)
n! * residue(exp(w*n)/w^n*1/(1-w)^3, w=0);
end;
EX3 :=
proc(n)
1/2*n! * add(n^q/q!*(n+1-q)*(n-q), q=0..n-1);
end;
Addendum Mon Sep 22 2019. A better way to extract the coefficients
from $H(z)$ is not to compute the derivative of $T(z)$ and use
$$Q_n = n! [z^n] H(z) =
n! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n+1}}
\frac{z T'(z)}{(1-T(z))^3} \; dz
\\ = n! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n}}
\frac{T'(z)}{(1-T(z))^3} \; dz.$$
Put $T(z)=w$ so that $z=w/\exp(w) = w\exp(-w)$ and
$dw = T'(z) \; dz$ to get
$$n! \frac{1}{2\pi i}
\int_{|w|=\gamma} \frac{\exp(nw)}{w^{n}}
\frac{1}{(1-w)^3} \; dw.$$
We then continue as before. The standard Lagrange inversion does not
help here, we get
$$Q_n = (n-1)! [z^{n-1}] H'(z) \\ =
(n-1)! \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{n}}
\left( \frac{1}{(1-T(z))^4}
+ \frac{4T(z)}{(1-T(z))^5} \right) T'(z) \; dz
\\ = (n-1)! \frac{1}{2\pi i}
\int_{|w|=\gamma} \frac{\exp(nw)}{w^{n}}
\left( \frac{1}{(1-w)^4}
+ \frac{4w}{(1-w)^5} \right) \; dw.$$
This yields
$$(n-1)! \sum_{q=0}^{n-1} \frac{n^q}{q!} {n-1-q+3\choose 3}
+ 4(n-1)! \sum_{q=0}^{n-2} \frac{n^q}{q!} {n-2-q+4\choose 4}
\\ = (n-1)! \sum_{q=0}^{n-1} \frac{n^q}{q!} {n-q+2\choose 3}
+ 4(n-1)! \sum_{q=0}^{n-2} \frac{n^q}{q!} {n-q+2\choose 4}.$$
With
$${n-q+2\choose 3} + 4 {n-q+2\choose 4}
= {n-q+2\choose 3} + (n-q-1) {n-q+2\choose 3}$$
we find
$$(n-1)! \sum_{q=0}^{n-1} \frac{n^q}{q!}
(n-q) {n-q+2\choose 3}.$$
Tail length.
This requires a modification to the tree function in order to count
the total tail length for all nodes in the tree.
We obtain the functional equation
$$T(z, v) = z\sum_{q\ge 0} \frac{T(vz,v)^q}{q!}
= z\exp(T(vz, v)).$$
This is because when we attach a set of trees to the new root all
nodes in the tree have their tail length incremented by one.
We then have the a simple species for average tail lengths:
$$\textsc{SET} \left(\textsc{CYC}(\mathcal{T})\right).$$
It follows that the bivariate generating function of mappings by node
count and tail length is
$$G(z, v) = \exp\log\frac{1}{1-T(z,v)}
= \frac{1}{1-T(z,v)}.$$
This yields for the generating function $H(z)$ of the expected tail
length (same as before)
$$H(z) = \left.\frac{\partial}{\partial v} G(z,v)\right|_{v=1}.$$
This works out to
$$H(z) = \left.\frac{1}{(1-T(z,v))^2}
\frac{\partial}{\partial v} T(z,v)
\right|_{v=1}
= \frac{1}{(1-T(z,v))^2}
\left.z\exp T(vz, v)
\left(z \frac{\partial}{\partial z} T(z,v)+
\frac{\partial}{\partial v} T(z,v)
\right)\right|_{v=1}
\\ = \frac{T(z)}{(1-T(z))^2}
(z T'(z) + H(z) (1-T(z))^2)
\\ = \frac{z T(z) T'(z)}{(1-T(z))^2} + H(z) T(z).$$
This yields
$$H(z) = \frac{z T(z) T'(z)}{(1-T(z))^3}
= \frac{T(z)^2}{(1-T(z))^4}.$$
The Lagrange inversion computation can be carried out as before with
an extra $w$ in the integrand and it produces
$$Q_n = \frac{1}{2}
n! \sum_{q=0}^{n-2} \frac{n^q}{q!}
(n-q) (n-1-q).$$
The sequence $Q_n$ starts with
$$0, 2, 36, 624, 11800, 248400, 5817084, 150660608,
\\ 4285808496, 133010784000,\ldots$$
The Maple code for verification by total enumeration of this sequence
was as follows.
Q :=
proc(n)
option remember;
local ind, d, gf, pos, q, x, seen, traj, cycinit;
if n = 1 then return 1 fi;
gf := 0;
for ind from n^n to 2*n^n-1 do
d := convert(ind, base, n);
d := map(l->l+1, [seq(d[q], q=1..n)]);
q := 0;
for pos to n do
seen := {}; x := pos; traj := [];
while not(x in seen) do
traj := [op(traj), x];
seen := seen union {x};
x := d[x];
od;
cycinit := 1;
while traj[cycinit] <> x do
cycinit := cycinit + 1;
od;
q := q + cycinit-1;
od;
gf := gf+v^q;
od;
gf;
end;
EX :=
proc(n)
local T;
T := solve(TT=z*exp(TT), TT);
n!*coeftayl(T^2/(1-T)^4, z=0, n);
end;
EX2 :=
proc(n)
n! * residue(exp(w*n)/w^(n-1)*1/(1-w)^3, w=0);
end;
EX3 :=
proc(n)
1/2*n! * add(n^q/q!*(n-q)*(n-1-q), q=0..n-2);
end;
The labelled tree function recently appeared at this
MSE link
and at this
MSE link II
which is closely related to the present subject.