Let $H_n = 1 + {1 \over 2} + \cdots + {1 \over n}$ be the $n$-th partial sum of the harmonic series. For a given prime $p$, denote by $J_p$ the set of $n$ for which $p$ divides the numerator of $H_n$. In 1991 {\it A. Eswarathasan} and {\it E. Levine} [Discrete Math. 91, 249-257 (1991;

Zbl 0764.11018)] have determined $J_p$ for $p \in \{2,3,5,7\}$ and made the conjecture that $|J_p |$ is finite for all $p$. However, they didn’t prove even that $|J_{11} |$ is finite. The author remarks that this fact appears as a problem in [{\it R. L. Graham}, {\it D. E. Knuth}, {\it O. Patashnik}, Concrete Mathematics (Addision-Wesley, 1989;

Zbl 0668.00003)], and shows that this set contains exactly 638 integers, the largest of which is a number of 31 decimal digits. He determines $J_p$ for all $p < 550$, with three exceptions: 83, 127, 397. In this eye-opening paper the author strengthens the above conjecture on $|J_p |$, by using the theory of branching processes. This is based on a new $p$-adically convergent formula for $H_{pn} - H_n/p$. A probabilistic model predicts that $|J_p |= O (p^2 (\log \log p)^{2 + \varepsilon})$ and that $|J_p |\ge p^2 (\log \log p)^2$ for infinitely many $p$. Another interesting conjecture, supported by a probabilistic argument, is that the density of primes $p$ with $|J_p |= 3$ is $1/e$. This is confirmed experimentally for all $p \le 10^5$.