A question about generalized harmonic numbers modulo $p$

Question

Let $p \equiv 1 \pmod{3}$ be a prime and denote $H_{n,m} = \sum_{k = 1}^n 1/k^m$ as the $n,m$-th generalized harmonic number. I'm interested in computing $H_{(p-1)/3,\, 2}$ and $H_{(p-1)/6,\,2}$ modulo $p$. From this paper I know \begin{align*} H_{(p-1)/6,2} \equiv -\frac{B_{2p-3}(1/6)}{2p-3}\pmod p \end{align*} and \begin{align*} H_{(p-1)/3,2} \equiv -\frac{B_{2p-3}(1/3)}{2p-3}\pmod p \end{align*} where $B_n(x) = \sum_{k = 0}^n {n \choose k}x^{n-k}B_k$ is the $n$-th Bernoulli polynomial.

It is well known, I think, that \begin{align*} B_n(1/6) = (1 - 3^{1-n})(1 - 2^{1-n})\frac{B_n}{2\cdot 3^{n-1}} \end{align*} and \begin{align*} B_n(1/3)=(1 - 3^{1-n})\frac{B_n}{2\cdot 3^{n-1}} \end{align*} but only for even $n$. From this we can eventually express $H_{(p-1)/3} = H_{(p-1)/3, 1}$ and $H_{(p-1)/6} = H_{(p-1)/6,1}$ modulo $p$ in terms of the Fermat quotients $q_p(2)$ and $q_p(3)$ where $q_p(a) = (a^{p-1} - 1)/p$ for $a$ co-prime to $p$.

It is mentioned in the paper above that there is some kind of expression for for $B_n(1/3)$ and $B_n(1/6)$ when $n$ is odd in terms of "$I$ numbers", but I can't find the source, or anything else that is relevant.

Is it possible to write $H_{(p-1)/3,2}$ and $H_{(p-1)/6,2}$ in terms of Fermat quotients modulo $p$ analogous to $H_{(p-1)/3}$ and $H_{(p-1)/6}$? What are these "$I$ numbers?"?

I've asked a similar question here but I thought I would ask here.

Answer 1

Glaisher's I-numbers are described in J. W. L. Glaisher, On a set of coefficients analogous to the Eulerian numbers, Proc. London Math. Soc., 31 (1899), 216-235.

Answer 2

Here is what I've learned: \begin{align*} 2\sum_{n = 0}^{\infty} B_{2n+1}\left(\frac{1}{3}\right) \frac{x^{2n+1}}{(2n+1)!} & = \frac{xe^{\frac{1}{3}x}}{e^x -1} + \frac{xe^{-\frac{1}{3}x}}{e^{-x} -1} = \frac{xe^{\frac{1}{3}x}}{e^x -1} - \frac{xe^{\frac{2}{3}x}}{e^{x} -1} \\ &= \frac{x(e^{\frac{1}{3}x}-e^{\frac{2}{3}x})}{e^x -1} = \frac{xe^{\frac{1}{3}x}(1-e^{\frac{1}{3}x})}{(e^{\frac{1}{3}x} -1)(e^{\frac{2}{3}x}+e^{\frac{1}{3}x}+1)} \\ & = \frac{-xe^{\frac{1}{3}x}}{e^{\frac{2}{3}x}+e^{\frac{1}{3}x}+1} = \frac{-x}{e^{\frac{1}{3}x}+e^{-\frac{1}{3}x}+1} \end{align*}

The $n$-th $I$ number is defined by \begin{align*} \frac{3/2}{e^x + e^{-x} + 1} = \sum_{n = 0}^{\infty} I_n \frac{x^n}{n!} \end{align*}

Notice that \begin{align*} 2 \sum_{n = 0}^{\infty} I_{2n+1} \frac{x^{2n+1}}{(2n+1)!} &= \sum_{n = 0}^{\infty} I_n \frac{x^n}{n!} - \sum_{n = 0}^{\infty} I_n \frac{(-x)^n}{n!} \\ & = 0. \end{align*} Therefore, $I_n = 0$ whenever $n$ is odd, and $ \sum_{n = 0}^{\infty} I_n \frac{x^n}{n!} = \sum_{n = 0}^{\infty} I_{2n} \frac{x^{2n}}{(2n)!}$. Now,

\begin{align*} \frac{-x}{e^{\frac{1}{3}x}+e^{-\frac{1}{3}x}+1} &= -\frac{2}{3}x \frac{3/2}{e^{\frac{1}{3}x}+e^{-\frac{1}{3}x}+1} \\ & = -\frac{2}{3}x \sum_{n = 0}^{\infty} I_{2n} \frac{\left( \frac{1}{3}x\right)^{2n}}{(2n)!} \\ & = 2\sum_{n = 0}^{\infty} -\frac{(2n+1)I_{2n}}{3^{2n+1}} \frac{x^{2n+1}}{(2n+1)!}. \end{align*} Therefore, \begin{align*} \sum_{n = 0}^{\infty} -\frac{(2n+1)I_{2n}}{3^{2n+1}} \frac{x^{2n+1}}{(2n+1)!} = \sum_{n = 0}^{\infty} B_{2n+1}\left(\frac{1}{3}\right) \frac{x^{2n+1}}{(2n+1)!} \end{align*} which implies, \begin{align*} B_{2n+1}\left(\frac{1}{3}\right) = -\frac{(2n+1)I_{2n}}{3^{2n+1}}. \end{align*}