It is important to recognize the intuition being used in this problem, before we come to the actual proof.
In fact, a lot of probabilistic number theory is based around some very simple intuition and reasoning involving distinct primes.
First, let me put a setting in place. We consider the set $\{1,2,...,n\}$ as a sample space with the uniform probability measure.
On this, we can define the following random variable : for $m \in \mathbb N$, $I_m(x)=0$ if $m \nmid x$ and $I_m(x)=1$ if $m | x$. The question is : if $p_1\neq p_2$ are distinct primes, then what is the relationship between $I_{p_1}$ and $I_{p_2}$?
The answer is actually quite simple : for large enough $n$, we have $E[I_{p_1}] \approx \frac{1}{p_1}$ , $E[I_{p_2}] \approx \frac 1{p_2}$ and $$E[I_{p_1}I_{p_2}] = E[I_{p_1p_2}] \approx \frac 1{p_1p_2} \approx E[I_{p_1}]E[I_{p_2}]$$
and therefore, $I_{p_1}$ and $I_{p_2}$ are (we've approximated, since the probabilities are not exact) weakly dependent for large enough $n$. (The exact notion of weak independence is established by a concrete bound for $|E[I_{p_1}I_{p_2}] - E[I_{p_1}]E[I_{p_2}]|$, which is quite easy to calculate).
This particular idea is fundamental. Let's quickly see why.
Note : The idea quickly generalizes to a family of such variables i.e. for distinct primes $p_1,...,p_n$, the family $I_{p_1},...,I_{p_n}$ is "weakly" dependent in that none of the random variables affect the other to a great degree for large $n$, with a concrete result easily providable.
Consider the quantity $\Omega(n)$, the number of distinct prime divisors of $n$. Now, the question is : given a prime $p$,does $p$ divide $n$? Now, we look over all primes $p$. In other words, fixing a $m$, for $n>m$ the value of the random variable $\sum_{p \text{ prime}} I_p$ is $\Omega(m)$!
But now, because of what we've written earlier, $\sum_{p \text{ prime}} I_p$ is the sum of a weakly dependent collection of random variables.
Now, we all know that if the random variables were independent, then the CLT would apply. However, the CLT does apply in weaker conditions where the dependence is present but really low.
But let's imagine independence did apply. If it did,then let $S_n = \sum_{k=1}^{n} I_{p_k}$. We have $E[S_n] = \sum_{k=1}^n \frac 1{p_k}$ (doesn't use independence), and it's well known that this behaves like $\log \log p_k$ to the first order (a proof using the prime number theorem to a first order correction is one way of doing this). What about $Var(S_n)$? Well, IF the random variables were weakly dependent, THEN the variance is the sum of variances, and we've got indicator random variables so the variance is quite clearly $Var(S_n) = \sum_{k=1}^n \frac {p_k-1}{(p_k)^2}$, which once again (basically because $p_k-1 \approx p_k$ will cancel the $p_k$ in the denominator) will look like $\log \log p_k$.
Therefore, we know that if CLT applied then $\frac{S_n - E(S_n)}{\sqrt{Var(S_n)}} \to N(0,1)$. Applying that here, and noting that $\Omega(n) = S_n$, we have
$$
\frac{\Omega(n) - \log \log n}{\sqrt{\log \log n}} \to N(0,1)
$$
and now , all that is left to see is that some variant of the CLT is available that works for weakly dependent RVs.
What about $d(n)$ then? Here, we need a way of counting the multiplicity of a prime. What kind of random variable counts multiplicity of a prime dividing a uniformly randomly chosen number? Well, to see this let's fix a prime $p$ and a large integer $k$.
It's easy to see that in the interval $\{1,...,p^k\}$, exactly $1$ element has multiplicity $k$, $p-1$ have multiplicity $k-1$, $p(p-1)$ have multiplicity $k-2$ , ,... , $p^{k-1}(p-1)$ have multiplicity $0$. Writing this as a probability vector after dividing by $p^k$ gives :
$$
\left[1 - \frac 1p , \frac 1p\left(1 - \frac 1p\right), \frac 1{p^2}\left(1-\frac 1p\right) , ... , \frac{1}{p^k}\right]
$$
which, barring the $k$-multiplicity part, clearly resembles a geometric distribution with success probability $1 - \frac 1p$.
And therefore, define geometric random variables $G_p$ for each prime $p$, given by $G_p(k)$ as the multiplicity of $p$ dividing $k$.
Now, it is quite clear that $G_p$ and $G_{p'}$ are weakly dependent for $p \neq p'$ (Explicit bounds are elementary). But here's the interesting point : for large primes $p$, it's extremely unlikely that $p$ will even divide a randomly chosen number, and it's even more unlikely that the multiplicity is bigger than $1$. Therefore, for all purposes, we must imagine that for large $p$, the value of $G_p$ is almost always $0$ or $1$. Explicit bounds will be available, once again!
With that, let's think about $\log d(n)$. Because $d$ is weakly multiplicative and $d(p^k) = k+1$, it is quite clear that $\log d(n) = \sum_{i} \log (G_{p_i}+1)$.
Now, let's repeat the CLT rigmarole. $d(n)$ is like the sum of weakly dependent random variables, so let's find their expectation and variance.
But HOW? We've got a logarithm of a Geometric random variable, and let's say that calculations get complicated quite quickly!
So we go back to our intuition : recall what I said about the geometric random variable nearly always being $0$ or $1$? Well, this means that we should be looking to compare the lower order behaviours of $\log(1+G_{p_i})$ with the random variable $\log(1+I_{p_i})$, where we note that $G_{p_i} = I_{p_i}$ with really high probability for large primes $p_i$.
Therefore, with this intuition, the expectations and variance will hopefully match with that of $\log(1+I_{p_i})$ up to the first order! The expectation of $\log (1+I_{p_i})$ is $\frac{\log 2}{p_i}$. The variance is $(\log 2)^2\frac 1{p_i}\left(1 - \frac 1{p_i}\right)$. Taking the sum over the primes, we get approximately $\log 2 (\log \log n)$ and noting that $\sum \frac 1{p_i^2}<\infty$ we get approximately $(\log 2)^2 \log \log n$ for the variance.
Once again, the heuristics apply and if the dependencies and the bounds work out, we should get $\frac{\log d(n) - \log 2\log \log n}{(\log 2)\sqrt{\log \log n}} \to N(0,1)$.
Of course, that's the heuristics. We learnt a lot from here , the main points being :
Independence of prime division counting indicator random variables.
Imagining the CLT held and seeing what happens.
Powers via geometric random variables.
The log - geometric approximation by the aforementioned indicators.
It is a good instructive exercise to see what one can do with similar additive functionals, like the sum of divisors function.
So what about the main proof?
At this point, I felt like I would be able to elementarily derive moment bounds which are required to solve this question from the approximations I made above. However, to be honest I just got so lost in the calculation! So I decided to make use of machinery.
Machinery, which is actually very heavy, and found as Theorem 3.1 of
Billingsley, Patrick, The probability theory of additive arithmetic functions, Ann. Probab. 2, 749-791 (1974). ZBL0327.10055.
which I state verbatim (note that it's a vast generalization of what we actually need!).
Suppose that $f_n$ are completely additive functions i.e. $f_n(ab)=f_n(a)+f_n(b)$ for $a,b$ coprime, and define
$$
A_n = \sum_{p \leq n} \frac{f_n(p)}{p} \\
B_n = \sqrt{\sum_{p \leq n}\frac{f_n^2(p)}{p}}
$$
suppose , furthermore that
- $\lim_{m \to \infty} \frac{f_n(m)}{B_n} \to 0$.
- $\max_{p \leq n} \frac{|f_n(p)|}{B_n} \to 0$.
Then
$$
\frac{f_n-A_n}{B_n} \to N(0,1)
$$
However, the proof isn't staggeringly difficult from what we have discussed! Indeed, IF random variables supported on different primes were independent, then $A_n$ is the exact mean and $B_n$ is the approximate variance (once you make $1 - p^{-1} \approx 1$ for large primes $p$) of each $f_n$. The proof is literally nothing but making sure that the dependence is weak enough for the moments to reflect , up to their first order, the products of the individual moments.
So what one does is considers an alternate family of random variables $g_n$, which are actually independent. How? quite simple : we first make each $f_n$ into a sum-of-weakly-dependent random variables, then we replace each random variable by independent instances of the "rough" random variable that they "should" ideally be (so replace $I_p$ by an actual independent Bernoulli random variable with parameter $1-\frac 1p$, and likewise the Geometric case). Now, we get a sum of independent random variables, for which we apply a traditional CLT (turns out Lyapunov's CLT works). So the CLT implies that the moments for the independent sum converge.
Now, we just make sure that the moments for the independent and the dependent case differ by something that vanishes in the limit, and that gives the result.
Let's apply this to our situation : we have $f_n \equiv f = \log d$, the divisor function. Let's calculate $A_n$ and $B_n$ for this.
Indeed, we first have
$$
A_n = \sum_{p \leq n} \frac{\log d(p)}{p} = \log 2\sum_{p \leq n} \frac 1p = \log 2 \log \log n + O(1/n) \\
B_n^2 = (\log 2)^2\sum_{p \leq n} \frac 1p = (\log_2)^2 \log \log n + O(1/n)
$$
Fix an $m$, then $d(m)$ is constant while $B_n \to \infty$ so the first condition is quite clearly satisfied. The second maximum is just , well, $d(p)=2$ so it's $\frac{2}{B_n}$ which goes to $0$ as $n \to \infty$. Hence, the conditions are satisfied and we are done!
Let me know if :
