Let $\sigma_0(n)$ be the divisor counting function: $$\sigma_0(n) = \sum_{d \vert n} 1.$$ I ran some numerical experiments that showed when $p$ is prime, the function $\sigma_0(n)$ is equidistributed mod $p$. That is, for any residue class $a \mod p$, $$\lim_{X \to \infty} \dfrac{ \vert \{ n<X: \sigma_0(n) \equiv a \mod p \} \vert }{X} = \dfrac{1}{p}.$$
Is this fact correct? If so, could anyone sketch a proof / provide a reference for a proof?
$\newcommand{\Y}{\mathfrak{X}_p(X)}$I haven't checked all the details on the application of Selberg–Delange below, so it's possible something I am saying is nonsense, but even after accounting for the correction noted by Noam Elkies in the comments, it appears to be the case that $\sigma_0(n)$ is not, in general, equidistributed among the non-zero congruence classes mod $p$ unless $2$ is a primitive root modulo $p$ (see A001122 on OEIS). In fact, we have that $$\lim_{X \to \infty}\biggl(\frac{1}{\Y}\sum_{n\leqslant X} 1_{\sigma_0(n) \equiv a \mod p}\biggr) = \frac{1}{p-1}\biggl(1+\frac{1}{\delta_p}\sum_{\substack{\chi \neq \chi_0\\\chi(2) = 1}} \overline{\chi}(a)G_{\chi}(1)\biggr),$$ where the sum runs over all nontrivial Dirichlet character mod $p$ with $\chi(2) = 1$, $$\Y = \sum_{n\leqslant X} 1_{p \nmid \sigma_0(n)}, \qquad \delta_p = \lim_{p \to \infty} \frac{\Y}{X},$$ are the counting function and density of the set of $n$ with $p \nmid \sigma_0(n)$ and $G_\chi$ is defined below. Note if $n$ is squarefree, then $\sigma_0(n) = 2^{\omega(n)}$, where $\omega(n)$ is the prime-divisor counting function, and hence $p \nmid \sigma_0(n)$ for $p > 2$. Since the squarefrees have positive density, this implies that $\Y \gg X$ and $\delta_p > 0$.
When $2$ is a primitive root modulo $p$, the only character with $\chi(2) = 1$ is the trivial one, so this would imply equidistribution in that case. But, for example, if you take $p = 7$ or $p = 17$, and consult the table of values of Dirichlet characters mod $p$, then you'll find that the only characters which have $\chi(2) = 1$ are the trivial character and the quadratic character; from this it follows that $\sigma_0(n)$ is biased towards being a quadratic residue over being a nonresidue. The number of characters that will occur in this sum for a given prime $p$ is $1$ less than A001917 on OEIS.
Roughly speaking, this is because $\sigma_0(p) = 2$, and hence the characters which have $\chi(2) = 1$ give a main term contribution. Here's a sketch:
Guided by Weyl's criterion for $(\mathbb{Z}/p\mathbb{Z})^\times$, suppose $(a,p)=1$, and note that orthogonality of Dirichlet characters reads $$1_{n \equiv a \mod p} = \frac{1}{p-1}\sum_{\chi} \overline{\chi}(a) \chi(n),$$ where the sum is over all Dirichlet characters mod $p$. Thus, $$\frac{1}{\Y}\sum_{n\leqslant X} 1_{\sigma_0(n) \equiv a \mod p} = \frac{1}{p-1} + \frac{1}{p-1}\sum_{\chi\neq \chi_0} \overline{\chi}(a)\left( \frac{1}{\Y} \sum_{n\leqslant X} \chi(\sigma_0(n)) \right), \tag{$\star$}\label{star}$$ where we have separated the contribution of the trivial character. It thus suffices to study $$M_\chi(X) = \sum_{n\leqslant X} \chi(\sigma_0(n)), $$ which is a mean-value of multiplicative function. The characters which satisfy $M_\chi(X) = o(X)$ do not contribute a main term to \eqref{star} since $\Y \gg X$, while the characters with $M_\chi(X) \asymp X$ do.
Standard multiplicative number theory techniques apply here. To put this into effect, define the Dirichlet series $$F_\chi(s) = \sum_{n\geqslant 1} \frac{\chi(\sigma_0(n))}{n^s},$$ so that an application of Perron's formula gives that $$\sum_{n\leqslant X} \chi(\sigma_0(n)) = \frac{1}{2\pi i} \int_{2-i\infty}^{2+i\infty} F_\chi(s) X^s \frac{ds}{s}.$$ Investigating the Euler product of $F_\chi$, we find $$F_\chi(s) = \prod_p \biggl(\sum_{k=0}^\infty \frac{\chi(\sigma_0(p^k))}{p^{ks}}\biggr) = \prod_p \biggl(\sum_{k=0}^\infty \frac{\chi(k+1)}{p^{ks}}\biggr) = \zeta(s)^{\chi(2)} G_\chi(s),$$ where $G_\chi(s)$ is convergent in $\Re(s) > 1/2$. An application of the Selberg-Delange method (see Chapter 5 of Tenenbaum's "Introduction to analytic and probabilistic number theory", in particular Theorem 5.2 with $z = \chi(2)$, $N = 1$ and $F = F_\chi$) should then give that $$ \sum_{n\leqslant X} \chi(\sigma_0(n)) = \frac{X (\log X)^{\chi(2) - 1}}{\Gamma(\chi(2))} \biggl(G_\chi(1) + O\Big(\frac{1}{\log X}\Big)\biggr),$$ from which it is clear that if $\chi(2) \neq 1$, then the expression on the right is $o(X)$ while, if $\chi(2) = 1$, then we get something of size $\asymp X$. Putting this back into \eqref{star}, and concentrating on the main terms proves the claim.
