Super Congruences Involving Multiple Harmonic Sums and Bernoulli Numbers

Let $m$, $r$ and $n$ be positive integers. We denote by ${\bf k}\vdash n$ any tuple of odd positive integers ${\bf k}=(k_1,\dots,k_t)$ such that $k_1+\dots+k_t=n$ and $k_j\ge 3$ for all $j$. In this paper we prove that for every sufficiently large prime $p$ $$\sum_{\substack{l_1+l_2+\cdots+l_n=mp^r p\nmid l_1 l_2 \cdots l_n }} \frac1{l_1l_2\cdots l_n} \equiv p^{r-1} \sum_{{\bf k}\vdash n} C_{m,{\bf k}} B_{p-{\bf k}} \pmod{p^r}$$ where $B_{p-{\bf k}}=B_{p-k_1}B_{p-k_2}\cdots B_{p-k_t}$ are products of Bernoulli numbers and the coefficients $C_{m,{\bf k}}$ are polynomials of $m$ independent of $p$ and $r$. This generalizes previous results by many different authors and confirms a conjecture by the authors and their collaborators.