On Eswarathasan--Levine and Boyd's conjectures for harmonic numbers

We provide numerical evidence towards three conjectures on harmonic numbers by Eswarathasan–Levine and Boyd. Let $J_p$ denote the set of integers $n\ge 1$ such that the harmonic number $H_n$ is divisible by a prime $p$. The conjectures state that: $(i)$ $J_p$ is always finite and of the order $O(p^2(\log \log p{)}^{2+\epsilon })$; $(ii)$ the set of primes for which $J_p$ is minimal (called harmonic primes) has density $e^{-1}$ among all primes; $(iii)$ no harmonic number is divisible by $p^4$. We prove $(i)$ and $(iii)$ for all $p\le 16843$ with at most one exception, and enumerate harmonic primes up to~$50\cdot {10}^5$, finding a proportion close to the expected density. Our work extends previous computations by Boyd by a factor of about $30$ and $50$, respectively.