What is sum of the Bernoulli numbers?

Using divergent summation (Nörlund sum) one can find it numerically that it is surely $\zeta(2)-1 \approx 0.644934$.
I've asked this question a couple of years ago in the usenet (sc.math) in a more general form, where the bernoulli-numbers are cofactored with sequences of binomial-coefficients. One reader (Robert Israel) showed an analytical derivation using the psi-function, which gave as a result, that the results were systematically $\zeta(2)-1,\zeta(3)-1,\zeta(4)-1,...$ with that binomial-coefficients as cofactors. (See the link at the end below)

I've an old article (with many inexact and sloppy formulation due to my limited knowledge then) containing a compilation of similar sums on my website where in the appendix I've reproduced that proof.
On page 5 you can find the definition of a matrix, which contains the bernoulli-numbers and indicated cofactors; this matrix contains the coefficients of the integrals of the Bernoulli-polynomials (or, equivalently said, the coefficients of the Faulhaber-polynomials for the sums-of-like-powers) if you read it rowwise instead of columnwise.

Here is a link to the discussion in the google-news-archive

From the the definition of the Bernoulli numbers $B_n$, we know that, $ \frac{x}{e^x-1}=\sum_{n=0}^{\infty}\frac{B_nx^n}{n!} $ ........(1)

Multiply oth the sides of this equation by $x^ke^{-x}$ , so we get $ \frac{x^{k+1}e^{-x}}{e^x-1}=\sum_{n=0}^{\infty}\frac{B_nx^{n+k}e^{-x}}{n!} $

Integrate both the sides of this equation from zero to infinity,so we get $ \int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx=\int_{0}^{\infty}\sum_{n=0}^{\infty}\frac{B_nx^{n+k}e^{-x}}{n!}dx $

We note here that, definition given as per Eq.(1) is valid only for '$ 0 \leq x < 2 \pi $' and therefore the integral taken on both the sides w.r.t. '$x$' in the range of $0$ to $\infty$ is a valid step only at LHS but since infinite series at RHS is convergent only under above constraints of x , hence such integral on RHS will give rise to a divergent sum. But, I think, in the sense of analytic continuation only, the divergent sum on RHS can be said to have the finite value equal to value of the integral at LHS. But, direct calculation will always say that RHS is a divergent sum. Other similar example i can give here is, the sum of factorials and alternate sum of factorials of all integers. These sums are actually divergent, but a finite value can be assigned to it. This is a famous sum given by Euler and by Ramanujan also.

Now, proceeding further,

$ \int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx=\sum_{n=0}^{\infty}\frac{B_n}{n!}\int_{0}^{\infty}x^{n+k}e^{-x}dx $

$ \int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx=\sum_{n=0}^{\infty}\frac{B_n}{n!}\times (n+k)! $

$ \int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx=\sum_{n=0}^{\infty}\frac{B_n(n+k)!}{n!} $

So, finally we got $ \sum_{n=0}^{\infty}\frac{B_n(n+k)!}{n!}=\int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx $

So, let $I=\int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx$ This can be written equivalently as, $I=\int_{0}^{\infty}x^{k+1}e^{-2x}(\frac{1}{1-e^{-x}})dx$

Since, $e^{-x}<1 $ for $0<x<\infty $ , hence we can write, $I=\int_{0}^{\infty}x^{k+1}e^{-2x}(1+e^{-x}+e^{-2x}+e^{-3x}+e^{-4x}+......)dx$

$I=\int_{0}^{\infty}x^{k+1}e^{-2x}dx+\int_{0}^{\infty}x^{k+1}e^{-3x}dx+\int_{0}^{\infty}x^{k+1}e^{-4x}dx+\int_{0}^{\infty}x^{k+1}e^{-5x}dx+......$

Now, consider a general term of above integral as, $I_r=\int_{0}^{\infty}x^{k+1}e^{-rx}dx$

Substitute, $rx=t$ , so we get $ rdx=dt$ i.e. $dx=\frac{1}{r}dt$ and $x=\frac{t}{r}$ The limits of integration will remain the same, hence we finally get $I_r=\frac{1}{r^{k+2}}\int_{0}^{\infty}t^{k+1}e^{-t}dt$ which is $I_r=\frac{1}{r^{k+2}}\times (k+1)!$ So, finally we got, $I_r=\frac{(k+1)!}{r^{(k+2)}}$ Hence, from the form of above integral $I$ , we can write $I=\frac{(k+1)!}{2^{k+2}}+\frac{(k+1)!}{3^{k+2}}+\frac{(k+1)!}{4^{k+2}}+.... $

i.e. $I=(k+1)!(\zeta(k+2)-1) $ So, we finally got $ \sum_{n=0}^{\infty}\frac{B_n(n+k)!}{n!}=\int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx=(k+1)!(\zeta(k+2)-1) $

Dividing everything by 'k!' , we get,

So, we finally got $ \sum_{n=0}^{\infty}\frac{B_n(n+k)!}{k!n!}=\frac{1}{k!}\int_{0}^{\infty}\frac{x^{k+1}e^{-x}}{e^x-1}dx=(k+1)(\zeta(k+2)-1) $

As per definition of Binomial coefficient, we know that $\frac{(n+k)!}{k!n!}$ can be replaced by Binomial coefficient ${n+k}\choose{k}$

So, when $k=0$ , we get $\sum_{n=0}^{\infty} B_n=\int_{0}^{\infty}\frac{x^{1}e^{-x}}{e^x-1}dx=\zeta(2)-1 $

I have plotted $y=\frac{x^{k+1}e^{-x}}{e^x-1}$ for different values of 'k' I have also checked and confirmed the validity of the result $ \int_{0}^{\infty}\frac{x^{1}e^{-x}}{e^x-1}dx=\zeta(2)-1 $ which holds 100 percent true as the graph suggest that this integral do converge to finite value

This result is analogous in some sense to analytic continuation of a function - just as in case of Zeta function values at negative integers. If anyone could suggest me some good e-books about detailed analysis of such concepts , i will be very much thankful.