Regularizing divergent series and Bernoulli numbers

Here is the actual problem I need a proof for: (I made an error writing down the equation initially)

$$B_{k+1} = \frac{(-1)^k}{2^{k+2}-2} \sum_{q=0}^{k} \binom{k+1}{q} 2^q B_q$$

Below is my motivation why I am looking at the problem.

Kindly pardon the long post.

I was looking at some diverging series on how to regularize them to get a finite answer. For instance, $$1+1+\cdots = - \frac12$$ can be obtained through the analytic continuation based on geometric series (or) the zeta function evaluated at the origin. Similarly, we have that $$ \begin{array}{lr} 1^1 + 2^1 + \cdots = & -\frac1{12}\\ 1^2 + 2^2 + \cdots = & 0\\ 1^3 + 2^3 + \cdots = & \frac1{120}\\ 1^4 + 2^4 + \cdots = & 0 \end{array} $$ And in general, the zeta-regularization (or) analytic continuation based on geometric series gives us $$ \begin{align} 1^k + 2^k + \cdots & = \zeta(-k) \end{align} $$ However, they don't satisfy some nice relations/ recurrence satisfied by their finite counter parts. The first and foremost one which struck me was the one below. $$ \begin{align} 1^3 + 2^3 + 3^3 + \cdots + n^3 = \left( 1 + 2 + 3 + \cdots + n\right)^2 \end{align} $$ whereas $$ \zeta(-3) \neq \zeta(-1)^2 $$ Though there is no reason for $\zeta(-3)$ to equal $\zeta(-1)^2$, I was wondering if only we could assign a different value to the diverging series $1^3 + 2^3 + 3^3 + \cdots$ and $1 + 2 + 3 + \cdots$ so that we can still have $$ \begin{align} 1^3 + 2^3 + 3^3 + \cdots = (1+2+3+\cdots)^2 \end{align} $$ To do so, we turn out attention to their finite counterparts, $S(n;k) = 1^k + 2^k + \cdots n^k$. As mentioned earlier, we have some nice relations among $S(n;k)$ for different $k$'s. For instance, we have $$ \begin{array}{rrc} S(n;1) & = & \frac{S(n;0)(S(n;0)+1)}{2}\\ S(n;2) & = & \frac{S(n;1)(2S(n;0) + 1)}{3}\\ S(n;3) & = & S(n;1)^2\\ 1^3 + 2^3 + \cdots + n^3 & = & (1+2+ \cdots + n)^2\\ S(n;4) & = & \frac{S(n;2) (6S(n;1)-1)}{5} \end{array} $$ In general, we have $S(n;k) = P_k(S(n;0))$ where $P_k(a)$ is a polynomial in $a$ of degree $k+1$ and are called the Faulhaber polynomials and given by $$P_k(n) = \frac{n^{k+1}}{k+1} \sum_{j=0}^{k} \left[ \binom{k+1}{j} \left( - \frac1n \right)^{j} B_j \right]$$ where $B_k$ are the Bernoulli numbers. For instance, $$ \begin{array}{rll} P_0(a) & = a\\ P_1(a) & = \frac{a^2+a}{2} & = \frac{P_0(a) \left( P_0(a)+1 \right)}{2}\\ P_2(a) & = \frac{a(a+1)(2a+1)}{6} & = \frac{\left( 2P_0(a)+1 \right)P_1(a)}{3}\\ P_3(a) & = \frac{a^2(a+1)^2}{4} & = \left( P_1(a) \right)^2 \end{array} $$ We would like to define $$1^k + 2^k + 3^k + \cdots = \lim_{n \rightarrow \infty} S(n;k) = s_k$$ in a way consistent with the recurrence satisfied by the finite sums i.e. we would want $$s_k = \lim_{n \rightarrow \infty} S(n;k) = \lim_{n \rightarrow \infty} P_k(S(n;0)) = P_k \left( \lim_{n \rightarrow \infty} S(n;0) \right) = P_k(s_0)$$ As mentioned before, unfortunately, the analytic continuation based on geometric series (or) zeta-regularization doesn't satisfy this. For instance, $$\zeta(-1) = -\frac1{12} \neq -\frac18 = \frac{-\frac12 (-\frac12 + 1)}{2} = \frac{\zeta(0) (\zeta(0) + 1)}{2} = P_1 \left( \zeta(0) \right)$$ Similarly, $$\zeta(-3) = \frac1{120} \neq \left( -\frac1{12} \right)^2 = \zeta(-1)^2$$ So I was wondering if we could find another regularization still respecting the recurrence relation satisfied by their finite counterparts subject to $\displaystyle s_0 = \lim_{n \rightarrow \infty} S(n;0) = -\frac12$. Let us look at what we would want our $s_k$'s to be if we want the recurrences to be satisfied. $$ \begin{align} s_1 = & P_1(s_0) = & \frac{s_0(s_0 + 1)}{2} = & -\frac18\\ s_2 = & P_2(s_0) = & \frac{(2s_0 + 1)s_1}{3} = & 0\\ s_3 = & P_3(s_0) = & \frac{s_0^2(s_0 + 1)^2}{2^2} = & \frac1{64} \end{align} $$ One nice thing to notice is that $P_k(-1/2) = 0$ whenever $k$ is even which coincides with $\zeta(-k)$ whenever $k$ is even. Hence, $s_{k} = P_k(s_0) = P_k(-1/2) = 0$, whenever $k$ is even. Proceeding on the same lines for the odd numbers, we obtain the table shown below. $$ \begin{array}{lcrlrrl} 1^1 + 2^1 + \cdots& \text{ }& P_1(-1/2) = & -\frac18& \text{ }& \zeta(-1) = & -\frac1{12}\\ 1^3 + 2^3 + \cdots& \text{ }& P_3(-1/2) = & \frac1{64}& \text{ }& \zeta(-3) = & \frac1{120}\\ 1^5 + 2^5 + \cdots& \text{ }& P_5(-1/2) = & -\frac1{128}& \text{ }& \zeta(-5) = & -\frac1{252}\\ 1^7 + 2^7 + \cdots& \text{ }& P_7(-1/2) = & \frac{17}{2048}& \text{ }& \zeta(-7) = & \frac1{240}\\ 1^9 + 2^9 + \cdots& \text{ }& P_9(-1/2) = & -\frac{31}{2048}& \text{ }& \zeta(-9) = & -\frac1{132}\\ 1^{11} + 2^{11} + \cdots& \text{ }& P_{11}(-1/2) = & \frac{691}{16384}& \text{ }& \zeta(-11) = & \frac{691}{32760}\\ 1^{13} + 2^{13} + \cdots& \text{ }& P_{13}(-1/2) = & -\frac{5461}{32768}& \text{ }& \zeta(-13) = & -\frac{1}{12}\\ 1^{15} + 2^{15} + \cdots& \text{ }& P_{15}(-1/2) = & \frac{929569}{1048576}& \text{ }& \zeta(-15) = & \frac{3617}{8160} \end{array} $$ and so on. However, there is a nice relationship which I observed empirically but am not able to prove $$\zeta(-k) = \frac{P_{k}(-1/2)}{2 - \frac1{2^{k}}}$$ From some results on $\zeta$ functions and Bernoulli numbers, we have$$\displaystyle \zeta(-k) = - \frac{B_{k+1}}{k+1}$$ and $$\displaystyle P_{k}(-1/2) = \frac1{k+1} \frac{(-1)^{k+1}}{2^{k+2}-2} \sum_{q=0}^{k} \binom{k+1}{q} 2^q B_q$$

Hence, it boils down to proving the following relationship $$B_{k+1} = \frac{(-1)^k}{2^{k+2}-2} \sum_{q=0}^{k} \binom{k+1}{q} 2^q B_q$$

Hence, now we not only have $S(n;k) = P_k(S(n;0))$ but also $s_k = P_k(s_0)$ where $$s_k = \lim_{n \rightarrow \infty} S(n;k) = \left( 2 - 2^{-k} \right) \zeta(-k)$$

Hence, if the above is true, instead of assigning the value of the diverging series $1^k + 2^k + 3^k + \cdots$ to $\zeta(-k)$, to carry over the nice recurrence relations satisfied by their finite counterparts, would it be nice to assign a value of $(2-2^{-k}) \zeta(-k)$ to the diverging series $1^k + 2^k + 3^k + \cdots$ so that for instance, $$1^3 + 2^3 + 3^3 + \cdots = (1+2 +3 + \cdots)^2$$ would still be true.

Apart from what I have argued above, is there any other reason why one would like to associate the value $(2-2^{-k}) \zeta(-k)$ to the diverging series $1^k + 2^k + 3^k + \cdots ?$

In this article, the authors give the recurrence relationship, I am looking for (though they seem to employ a different sign convention for the Bernoulli numbers). In fact, they have a general version of the recurrence relationship, $$B_{k+1} = \frac{(-1)^{k}}{n^{k+2}-n} \sum_{q=0}^{k} \left[ n^q \binom{k+1}{q} B_q \sum_{j=1}^{n-1} j^{k+1-q} \right]$$ $n=2$ gives the desired relationship, I am looking for.