Let $n = p_1^{a_1}p_2^{a_2}\cdots p_m^{a_m}$ be the prime factorization of $n$ and let $f(n) = a_1 a_2 \cdots a_m$. Using a heuristic argument I am able to show that the mean value of $\lim_{n \to \infty} \frac{1}{n}\sum_{k = 1}^n f(k)$ is close to
$$ \bigg(2 - \frac{1}{\zeta(2)}\bigg)\bigg(N_c + \frac{1}{\zeta(2)} - \frac{2}{\zeta(3)}\bigg) - \frac{1}{\zeta(2)} + \frac{2}{\zeta(3)} \approx 1.95979 $$
where $N_c \approx 1.7052$ is the Niven's constant. Experimental data seems to be consistent with this heuristics. For $n = 10^6$, the mean is $1.9406$ while for $n = 1.42 \times 10^{10}$, the mean increased slightly to $1.94357$.
Question: What is the limiting value of the mean? Is there a closed form for the mean of $f(n)$?
Let $h$ be a multiplicative function defined by $h(p^k)=k$ for all prime powers $p^k$. Then you are trying to calculate the mean value of $h$, and we have that:
This can be done using the general method in this answer: On the mean value of a multiplicative function: Prove that $\sum\limits_{n\leq x} \frac{n}{\phi(n)} =O(x) $
The idea is that since $h$ is "very close" to $1$, $f=\mu*h$ will be very close to zero, where $\mu$ is the Mobius mu function, and $*$ represents Dirichlet convolution. In particular, computing the Mobius inversion, we have $f(p)=0$, $f(p^k)=1$ for $k\geq 2$, and $(1*f)=h$.
From the method in the above answer, we have
Thus it follows that
$$\sum_{n\leq x} h(n) = x\sum_{d}\frac{f(d)}{d}+O(\sqrt{x}),$$ and hence the mean value is $$\sum_{d}\frac{f(d)}{d} = \prod_{p}\left(1+\frac{1}{p(p-1)}\right)$$
where the final form is an Euler-product. The Euler product can be rearranged as $$\prod_{p}\left(1+\frac{1}{p(p-1)}\right)=\prod_{p}\left(\frac{p^{2}-p+1}{p(p-1)}\right)$$ $$=\prod_{p}\left(\frac{p^{2}-p+1}{p(p-1)}\right)=\prod_{p}\left(\frac{p^{6}-1}{p(p^{2}-1)(p^{3}-1)}\right)=\frac{\zeta(3)\zeta(2)}{\zeta(6)}.$$
See also this list of similar questions using the same technique: