Here's a massive generalization, with references.
Consider a sequence $\{a_n \}_{n\ge 1}$ of integers. We'll say it satisfies the necklace congruences if $a_{n} \equiv a_{n/p} \mod n$ whenever $p \mid n$ ($p$ prime).
Here are some equivalent formulations:
- $a_{p^{k+1} m} \equiv a_{p^k m} \mod {p^{k+1}}$ for all $p,m,k$
- $\sum_{d \mid n} a_d \mu(n/d) \equiv 0 \mod n$ for all $n$
- The Mobius function can be replaced by any arithmetic function $f$ satisfying $\sum_{d \mid n} f(d) =0 \mod n, f(1) = \pm 1$.
So any such sequence gives rise to an integer sequence $\lambda_{\tilde{a}}(n) := \frac{1}{n}\sum_{d \mid n} a_d \mu(n/d)$ which satisfies $a_n = \sum_{d \mid n} d \lambda_{\tilde{a}}(n)$ by Mobius inversion. One can then define $\sigma_{\tilde{a}}(n):=\sum_{d \mid n, 1<d<n} d\lambda_{\tilde{a}}(d)$ and obtain the following generalization of Fermat's little theorem:
$a_n \equiv a_1 + \sigma_{\tilde{a}}(n) \mod n$
Now, let's go back to necklace congruences. A less obvious equivalence is the following:
Theorem: $\{a_n\}_{n \ge 1}$ satisfies the necklace congruences iff $\zeta_{a}(x):=\exp(\sum_{n\ge1} \frac{a_n}{n}x^n)$ has integer coefficients.
Proof: Write, $\zeta_a$ formally as $\prod_{n\ge 1} (1-x^n)^{-\frac{b_n}{n}}$. The $b_n$'s are uniquely determined, and in fact (by taking logarithmic derivative) it can be seen that they are $b_n = \sum_{d \mid n} a_d \mu(n/d)$. Now one direction is immediate and the other requires an inductive argument. $\blacksquare$
(Reference: Exercise 5.2 (and its solution) in Richard Stanley's book "Enumerative Combinatorics, vol. 2")
As Sergei remarked, for any integer square matrix $A$, the sequence $a_n := \text{Tr}(A^n)$ satisfies the necklace congruences. This is immediate from the last theorem, as the corresponding "zeta" function $\zeta_a$ is just:
$\exp( \sum_{n \ge 1} \frac{\text{Tr}((Ax)^n)}{n}) = \exp (\text{Tr} (- \ln (I -Ax)))=\det(I-Ax)^{-1} \in \mathbb{Z}[x]$
In particular, by taking a 1 on 1 matrix, we recover your observation.
Another generalization is $a_n = [x^n]f^n(x)$, where $f \in \mathbb{Z}[[x]]$. By taking the linear polynomial $f(x)=1+ax$ we recover your example, by taking $f(x)=(1+x)^m$ we get $a_n = \binom{nm}{n}$.
The really interesting feature is that this notion has a local version. Specifically, theorems of Dieudonne-Dwork and Hazewinkel give the following result:
Theorem: Given $\{ a_n \}_{n\ge 1} \subseteq \mathbb{Q}_{p}$ ($p$-adic numbers), the following are equivalent:
- $\zeta_{a}(x) := \exp(\sum_{n\ge 1} \frac{a_n }{n}x^n)\in\mathbb{Z}_{p}[x]$
- $\exp( \sum_{n \ge 1} p \frac{a_n - a_{n/p}}{n} x^n) \in 1+px\mathbb{Z}_{p}[x]$ (I defnite $a_{n/p}=0$ if $p\nmid n$)
- $\sum_{n \ge 1} \frac{a_n - a_{n/p}}{n} x^n) \in \mathbb{Z}_{p}[x]$
- $p \mid n \implies a_n \equiv a_{n/p} \mod {n\mathbb{Z}_{p}}$
Applying this simultaneously for all $p$, we recover the previous theorem.
A reference for this theorem is "A Course in p-adic Analysis" by Alain M. Robert (the theorems are stated there in a much more general context).
$\newcommand{floor}[1]{\left\lfloor #1 \right\rfloor}$
I doubt that this answer is still useful to you, since this question is one year old. Anyway, I'll leave the answer here, and maybe it will help someone (or not).
Your congruence is, modulo details, a special case of a congruence due to Voronoi. Let $k$ be an even positive integer, and write $B_k=\dfrac{N_k}{D_k}$ in lowest terms with $D_k>0$. Let $a$ and $n$ be coprime positive integers. Voronoi congruence states [2, Proposition 9.5.20]
that
\begin{equation}
\tag{1}\label{voronoin}
(a^k-1)N_k\equiv D_kka^{k-1}\sum_{x=1}^{n-1} x^{k-1}\floor{\frac{ax}{n}} \quad(\mathrm{mod}\ n),
\end{equation}
where $\floor{x}$ is the usual floor function.
Note that both sides are integers.
Let $n=p$ a prime number, and suppose $D_k$ is not divisible by $p$. Thus $p\ge5$ since $D_k$ is always divisible by $6$, by Clausen-von Staudt which states that [2, Corollary 9.5.15]
$$D_k=\prod_{(\ell-1)\mid k}\ell,$$
the product being over all primes $\ell$ such that $\ell-1$ divides $k$. Then \eqref{voronoin}
can be written as
\begin{equation}
\tag{2}\label{voronoip}
(a^k-1)B_k\equiv ka^{k-1}\sum_{x=1}^{p-1} x^{k-1}\floor{\frac{ax}{p}} \quad(\mathrm{mod}\ p).
\end{equation}
Let $c$ be a positive integer (coprime to $p$)
such that $ac\equiv 1 \ (\mathrm{mod}\ p)$, and for simplicity write $a=c^{-1}$, being an inverse modulo $p$. (Obviously $c$ can be chosen such that $0<c<p$ as in your statement.)
Multiplying both sides in \eqref{voronoip} by $c^k$ we obtain
\begin{equation}
\tag{3}\label{voronoic}
(1-c^k)B_k\equiv k\sum_{x=1}^{p-1} x^{k-1}c\floor{\frac{c^{-1}x}{p}} \quad(\mathrm{mod}\ p).
\end{equation}
Now, $cc^{-1}=1+pN$ for some integer $N$, and we use this as follows to relate the summands with $h_c(x)$.
For $1\le x\le p-1$, so that $x=\left(x \, \mathrm{mod}\, p\right)$, we have
\begin{align*}
c\floor{\frac{c^{-1}x}{p}}
&=\frac{cc^{-1}x}{p}-\frac{c\left(c^{-1}x \; \mathrm{mod}\, p\right)}{p}
=\frac{x}{p}+Nx-\frac{c\left(c^{-1}x \; \mathrm{mod}\, p\right)}{p}\\
&=h_c(x)+Nx-\frac{c-1}{2},
\end{align*}
so that \eqref{voronoic} becomes
\begin{equation}
\tag{4}\label{voronoitemp}
(1-c^k)B_k\equiv k\sum_{x=1}^{p-1} x^{k-1}h_c(x)+
kN\sum_{x=1}^{p-1}x^{k}
-k\frac{c-1}{2}\sum_{x=1}^{p-1}x^{k-1}
\quad(\mathrm{mod}\ p).
\end{equation}
To get rid of the last two sums, recall that for all $q$ not divisible by $p-1$ one has [1, Lemma 2.5.1]
\begin{equation}
\notag%\tag{4}\label{modp3}
\sum_{x=1}^{p-1}x^{q}\equiv 0
\quad(\mathrm{mod}\ p).
\end{equation}
Now we use your hypothesis $2\le k\le p-3$; then both $k-1$ and $k$ are not divisible by $p-1$, and the above implies that the last two sums in \eqref{voronoitemp} are zero. Also under this hypothesis, Clausen-von Staudt implies that $p$ does not divide $D_k$, which was used in \eqref{voronoip}.
If moreover $c^k \not\equiv 1 \pmod p$ we finally arrive to your congruence
\begin{equation}
\notag%\tag{4}\label{modp3}
B_k\equiv
\frac{k}{1-c^k}\sum_{x=1}^{p-1} x^{k-1}h_c(x)\quad(\mathrm{mod}\ p).
\end{equation}
[1] Cohen, Henri Number theory. Vol. I. Tools and Diophantine equations. Graduate Texts in Mathematics, 239. Springer, New York, 2007.
[2] Cohen, Henri Number theory. Vol. II. Analytic and modern tools. Graduate Texts in Mathematics, 240. Springer, New York, 2007.
(The original article by Voronoi can be accesed here. His congruence appears as Corollary IV in page 146.)
Best Answer
A proof is essentially given in Section 5.1 of Notes on primitive lambda-roots by P. J. Cameron and D. A. Preece.