Sums of integer powers similar to Prouhet–Tarry–Escott problem Recall Prouhet–Tarry–Escott problem. Its solution shows that certain sums of powers of integers can be made to vanish simultaneously if their signs are chosen to follow the Thue–Morse sequence, e.g.
$$\sum_{k=0}^{2^4-1}(-1)^{\sigma_2(k)}\,k^n=0\quad\text{for}\;n<4,\quad\text{i.e.}\\
\small
\begin{array}{l}
 0^1+3^1+5^1+6^1+9^1+10^1+12^1+15^1=1^1+2^1+4^1+7^1+8^1+11^1+13^1+14^1 \\
 0^2+3^2+5^2+6^2+9^2+10^2+12^2+15^2=1^2+2^2+4^2+7^2+8^2+11^2+13^2+14^2 \\
 0^3+3^3+5^3+6^3+9^3+10^3+12^3+15^3=1^3+2^3+4^3+7^3+8^3+11^3+13^3+14^3, \\
\end{array}$$
where $\sigma_2(k)$ is the sum of digits of $k$ in base-2.
I am interested in sums of powers of a similar form, but with larger exponents (such that the sum does not vanish anymore). I was able to make a few conjectures:
$$\begin{align}
&\sum _{k=0}^{2^n-1} (-1)^{\sigma _2(k)}\, k^{n+2}\stackrel{\color{gray}?}=\frac{2^{\binom{n}{2}}}{36}\, (-1)^n\, \left(2^n-1\right)\, \left(5\times2^n-4\right)\, (n+2)!\\
\\
&\sum _{k=0}^{2^n-1} (-1)^{\sigma _2(k)}\,k^{n+3}\stackrel{\color{gray}?}=\frac{2^{\binom{n}{2}}}{72}\,(-1)^n \,\left(2^n-1\right)^2\,\left(2\times2^n-1\right)\,(n+3)!\\
\\
&\small\sum _{k=0}^{2^n-1} (-1)^{\sigma_2(k)}\,k^{n+4}\stackrel{\color{gray}?}=\frac{2^{\binom{n}{2}}}{32400}\,(-1)^n\,\left(2^n-1\right)\,\left(143\!\times\! 8^n-307\!\times\!4^n+193\!\times\! 2^n-32\right)\,(n+4)!\end{align}$$
This pattern seems to continue within increasingly complicated polynomials of $2^n$ on the right-hand side. Can you suggest an approach to prove these, and to find a more general formula for an arbitrary exponent $k^{n+p}$?
 A: The exponential generating function for $S_{n,m}=\displaystyle\sum_{k=0}^{2^n-1}(-1)^{\sigma_2(k)}k^m$ (w.r.t. $m$) is $$\sum_{m=0}^\infty S_{n,m}\frac{x^m}{m!}=\sum_{k=0}^{2^n-1}(-1)^{\sigma_2(k)}e^{kx}=\prod_{k=0}^{n-1}(1-e^{2^k x})=2^{n(n-1)/2}(-x)^n\exp\sum_{k=0}^{n-1}f(2^k x),$$ where $$f(x)=\log\frac{e^x-1}{x}=\sum_{m=1}^\infty\frac{B_m^+}{m}\frac{x^m}{m!}$$ using the Bernoulli numbers (look at $f'(x)$ to see the expansion). Thus,
\begin{align*}
S_{n,n+p}&=(-1)^n 2^{n(n-1)/2}(n+p)!\ A_p(2^n),
\\A_p(a)&=[x^p]\exp\sum_{m=1}^\infty\frac{a^m-1}{2^m-1}\frac{B_m^+}{m}\frac{x^m}{m!}.
\end{align*}
This gives an algorithmic way to compute $S_{n,n+p}$ for any fixed $p$. Continuing your table,
\begin{align*}
  A_0(a)&=1
\\2A_1(a)&=a-1
\\36A_2(a)&=(a-1)(5a-4)
\\72A_3(a)&=(a-1)^2(2a-1)
\\32400A_4(a)&=(a-1)(143a^3-307a^2+193a-32)
\\64800A_5(a)&=(a-1)^2(2a-1)(19a^2-24a+2)
\\85730400A_6(a)&=(a-1)(5765a^5-19372a^4+22670a^3-10405a^2+1355a+32)
\\342921600A_7(a)&=(a-1)^2(2a-1)(1166a^4-2850a^3+1715a^2+60a-1)
\end{align*}
Observe that $A_p(a)$ is always divisible by $a-1$, and even by $(a-1)^2(2a-1)$ if $p>1$ is odd (because of $B_m^+=0$ for odd $m>1$, and [the g.f. of] $A_p(1/2)$ in closed form).
