It is well known that the sum of the odd powers of the $n$ first natural numbers is divisible by $\frac{n(n+1)}{2}$. See here for a simple proof, and here or more recently here for related questions.
I have observed (numerically) that it actually seems to be divisible by $\frac{n^2(n+1)^2}{4}\frac{1}{\mathrm{rad}\big(\frac{n(n+1)}{2}\big)}$, where $\mathrm{rad}\big(\frac{n(n+1)}{2}\big)$ is the product of the distinct prime factors of $\frac{n(n+1)}{2}$.
I am trying to derive an elementary proof involving modular arithmetics (Euler's theorem and/or tricks like in the proof here, ...) and I would appreciate any help.
EDIT here is a tentative proof that I have found, which is not as simple as I expected, as it requires Kummer theorem for the $q$-adic order of the binomial coefficient, Faulhaber formula for the sum of powers, involving Bernoulli numbers and the Von Staudt-Clausen theorem on the prime divisors of their denominators. I would like a verification of this proof, checking for any gap, and suggestions for a simpler proof.
We suppose $k\ge 1$. We first observe that it suffices to show that $S_{2k+1}(n):=\sum_{j=1}^nj^{2k+1}$ is divisible respectively by $\frac{n^2}{\mathrm{rad}n}$, when $n$ is odd or by $\frac{\frac{n^2}{4}}{\mathrm{rad}\frac{n}{2}}$, when $n$ is even. Indeed $n$ and $n+1$ are coprime and $S_{2k+1}(n)= S_{2k+1}(n+1)-(n+1)^{2k+1}$ is then also divisible by $\frac{\frac{(n+1)^2}{4}}{\mathrm{rad}(\frac{n+1}{2})}$ or$\frac{(n+1)^2}{\mathrm{rad}(n+1)}$ respectively, since $(n+1)^{2k+1}$ is divisible by $(n+1)^{3}$, hence by $\frac{\frac{(n+1)^2}{4}}{\mathrm{rad}(\frac{n+1}{2})}$ or$\frac{(n+1)^2}{\mathrm{rad}(n+1)}$ respectively.
Let $q$ be a prime number and $v_q(n)$ be the $q$-adic order of $n$. So, we want to show that when $q$ is odd and $v_q(n)=\alpha \ge 1$, then $v_q(S_{2k+1}(n))\ge 2\alpha -1$ and that when $v_2(n)=\alpha \ge 1$, then $v_2(S_{2k+1}(n))\ge 2\alpha -2$.
We write the Faulhaber formula involving the Bernoulli numbers $B_h$: $$S_{2k+1}(n)=\frac{n^{2k+2}}{2k+2}+\frac{n^{2k+1}}{2}+\sum_{h=1}^k\frac{{2k+2\choose 2h}}{2k+2}B_{2h}n^{2k+2-2h}.$$ We first suppose that $q$ is odd and does not divide $k+1$.
We clearly see that $\frac{n^{2k+2}}{2k+2} \equiv 0 \bmod q^{4\alpha}$, $\frac{n^{2k+1}}{2} \equiv 0 \bmod q^{3\alpha}$ and $\frac{n^{2k+2-2h}}{2k+2} \equiv 0 \bmod q^{2\alpha}$ for $h\le k$. Then since by the Von Staudt-Clausen Theorem, $v_q(B_{2h}) \ge -1$, we have $v_q(S_{2k+1}(n))\ge 2\alpha -1$.
Then we suppose that $q$ odd divides $k+1$ and $v_q(k+1) =\beta \ge 1$.
We have $v_q(\frac{n^{2k+2}}{2k+2})= 2k\alpha +2 \alpha - \beta $. Since $k+1 \ge q^\beta$, $v_q(\frac{n^{2k+2}}{2k+2}) \ge 2\alpha q^\beta -\beta$, but $q^\beta \ge \beta $, then $v_q(\frac{n^{2k+2}}{2k+2}) \ge \beta (2\alpha-1) \ge 2\alpha-1$ and we still have $v_q( \frac{n^{2k+1}}{2}) \ge 3\alpha$.
Now by Kummer theorem, $v_q({2k+2 \choose 2h})$ is the number of carries when adding $2h$ and $2k+2-2h$ in base $q$. In the base $q$ representation, the $\beta$ first digits of $2k+2$ are zero. When $q$ does not divide $2h$, it does not divides $2k+2-2h$ either, and $2h$ and $2k+2-2h$ cannot have a zero digit at the same position in their $\beta$ first digits. Then there are at least $\beta$ carries to obtain the $\beta$ first zero digits when adding $2h$ and $2k+2-2h$ and in this case where $q$ does not divide $2h$, we have $v_q({2k+2 \choose 2h}) \ge \beta$. More generally, we suppose that $v_q(2h)=v_q(h)= \gamma$, then $v_q(2k+2-2h)=v_q(k+1-h)= \mathrm{min}(\gamma,\beta)$ and there are at least $\beta-\mathrm{min}(\gamma,\beta) $ carries to obtain the $\beta$ first zero digits when adding $2h$ and $2k+2-2h$ and in the general case, we have $v_q({2k+2 \choose 2h}) \ge \beta-\mathrm{min}(\gamma,\beta)$. Then
$$v_q\Big(\frac{{2k+2 \choose 2h}}{2k+2}\Big) \ge -\mathrm{min}(\gamma,\beta) $$ Now $v_q\big(n^{2k+2-2h}\big) = 2\alpha(k+1-h) $ and $ k+1-h \ge q^{v_q(k+1-h)}= q^{\mathrm{min}(\gamma,\beta)}$. Then $$v_q\Big(\frac{{2k+2 \choose 2h}}{2k+2}n^{2k+2-2h}\Big) \ge 2\alpha\cdot q^{\mathrm{min}(\gamma,\beta)} -\mathrm{min}(\gamma,\beta)$$ but $0\le x\le q^x -1$ then $$v_q\Big(\frac{{2k+2 \choose 2h}}{2k+2}n^{2k+2-2h}\Big) \ge (2\alpha-1)\cdot q^{\mathrm{min}(\gamma,\beta)} +1 \ge 2\alpha.$$ Then again by the Von Staudt-Clausen Theorem, $v_q(B_{2h}) \ge -1$ and we have $v_q(S_{2k+1}(n))\ge 2\alpha -1$.
We now suppose that $q=2$ and does not divide $k+1$.
We clearly see that $\frac{n^{2k+2}}{2k+2} \equiv 0 \bmod 2^{4\alpha-1}$, $\frac{n^{2k+1}}{2} \equiv 0 \bmod 2^{3\alpha-1}$ and $\frac{n^{2k+2-2h}}{2k+2} \equiv 0 \bmod 2^{2\alpha-1}$ for $h\le k$. Then since by the Von Staudt-Clausen Theorem, $v_2(B_{2h}) \ge -1$, we have $v_2(S_{2k+1}(n))\ge 2\alpha -2$.
Then we suppose that $2$ divides $k+1$ and $v_2(k+1) =\beta \ge 1$.
We have $v_2(\frac{n^{2k+2}}{2k+2})= 2k\alpha +2 \alpha - \beta-1 $. Since $k+1 \ge 2^\beta$, $v_2(\frac{n^{2k+2}}{2k+2}) \ge 2\alpha 2^\beta -\beta-1$, but $2^\beta \ge \beta $, then $v_2(\frac{n^{2k+2}}{2k+2}) \ge \beta (2\alpha-1)-1 \ge 2\alpha-1$ and we still have $v_2( \frac{n^{2k+1}}{2}) \ge 3\alpha-1$.
Again by Kummer theorem, $v_2({2k+2 \choose 2h})$ is the number of carries when adding $2h$ and $2k+2-2h$ in base $2$. In the base $2$ representation, the $\beta+1$ first digits of $2k+2$ are zero. We suppose that $v_2(2h)=v_2(h)+1= \gamma+1$, then $v_2(2k+2-2h)=v_2(k+1-h)+1= \mathrm{min}(\gamma,\beta)+1$ and there are at least $\beta+1-\mathrm{min}(\gamma,\beta)-1 $ carries to obtain the $\beta+1$ first zero digits when adding $2h$ and $2k+2-2h$ and we have $v_2({2k+2 \choose 2h}) \ge \beta-\mathrm{min}(\gamma,\beta)$. Then
$$v_2\Big(\frac{{2k+2 \choose 2h}}{2k+2}\Big) \ge -\mathrm{min}(\gamma,\beta)-1 $$ Now $v_2\big(n^{2k+2-2h}\big) = 2\alpha(k+1-h) $ and $ k+1-h \ge 2^{v_2(k+1-h)}= 2^{\mathrm{min}(\gamma,\beta)}$. Then $$v_2\Big(\frac{{2k+2 \choose 2h}}{2k+2}n^{2k+2-2h}\Big) \ge 2\alpha\cdot 2^{\mathrm{min}(\gamma,\beta)} -\mathrm{min}(\gamma,\beta)-1$$ but $0\le x\le 2^x -1$ then $$v_2\Big(\frac{{2k+2 \choose 2h}}{2k+2}n^{2k+2-2h}\Big) \ge (2\alpha-1)\cdot 2^{\mathrm{min}(\gamma,\beta)} \ge 2\alpha-1.$$ Then again by the Von Staudt-Clausen Theorem, $v_2(B_{2h}) \ge -1$ and we have $v_2(S_{2k+1}(n))\ge 2\alpha -2$.
