Let

${H}_{n}=1+\frac{1}{2}+\cdots +\frac{1}{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

*A. Eswarathasan* and

*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 [

*R. L. Graham*,

*D. E. Knuth*,

*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\left({p}^{2}{(loglogp)}^{2+\epsilon}\right)$ and that

$|{J}_{p}|\ge {p}^{2}{(loglogp)}^{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}$.