Eswarathasan–Levine conjecture on infinitely many harmonic primes

Let $H_n=1+\frac12+\cdots+\frac1n$ be the $n$th harmonic number. For a prime $p$, define

$$J_p:=\{n\geq 1:\nu_p(H_n)\geq 1\},$$

where $\nu_p$ is the $p$-adic valuation. For every $p\geq 5$, the set $J_p$ contains $p-1$, $p^2-p$, and $p^2-1$. A harmonic prime is a prime $p$ such that $|J_p|=3$.

Eswarathasan–Levine conjecture. There are infinitely many harmonic primes.

The conjecture was suggested by Eswarathasan and Levine. Boyd's probabilistic model gives the quantitative refinement that harmonic primes should have density $e^{-1}$ among all primes; neither this refinement nor the infinitude assertion is established in the source.