# Limit approach to finding $1+2+3+4+\ldots$

When exploring the divergent series consisting of the sum of all natural numbers

$$\sum_{k=1}^\infty k=1+2+3+4+\ldots$$

I came across the following identity involving a one-sided limit:

$$\lim_{x\to0^-}\sum_{k=1}^\infty k\exp(kx)\cos(kx)=-\frac{1}{12}$$

Using zeta function regularization yields the same value:

$$\sum_{k=1}^\infty k^{-s}=\zeta(s)$$ $$\zeta(-1)=-\frac{1}{12}$$

In general, I found that for $$y\neq0$$ and $$n=1,5,9,13,\ldots$$ (i.e. $$4m+1$$ where $$m\in\mathbb{N}$$),

$$\lim_{x\to0^-}\sum_{k=1}^\infty k^n\exp(kxy)\cos(kxy)=\zeta(-n)$$

However, a similar limit did not exist for other powers of $$k$$, e.g.

$$\lim_{x\to 0^-}\sum_{k=1}^\infty k^2\exp(kx)\cos(kx)$$ $$\lim_{x\to 0^-}\sum_{k=1}^\infty k^3\exp(kx)\cos(kx)$$

The regularized values of the corresponding series are $$\zeta(-2)=0$$ and $$\zeta(-3)=\frac{1}{120}$$.

Given this information, I have the following questions:

1. What is the connection between the limit approach and the zeta function approach?
2. Why does the limit expression only seem to converge for $$n=4m+1$$?
3. Can the limit approach be used to find the sum for other powers of $$k$$? If so, how?

Here is a plot I made with Mathematica using the command

Plot[Evaluate[Sum[k*Exp[x*k]*Cos[x*k], {k, 1, Infinity}]], {x, -16, 16}, PlotRange -> {-.25, .25}, AspectRatio -> 1]


Notice how it approaches $$-1/12\approx-0.08333$$ as $$x$$ approaches $$0$$.

Further information:

$$\sum_{k=1}^\infty k^n\exp(kx)=\mathrm{Li}_{-n}(\exp(x))=\frac{n!}{(-x)^{n+1}}+\zeta(-n)+O(x)$$

Edit:

Based on Micah's answer, we seek an $$f$$ such that $$f(s,0) = 1$$ and \begin{align} \int_0^\infty x^s f(s,x) \,\mathrm{d}x &= 0 \end{align}

Let \begin{align} f(s,x) &= \mathrm{e}^{-x}(1+ax) \end{align}

Then $$f(s,0) = 1$$. Assuming $$\mathrm{Re}(s) > -1$$, \begin{align} 0 &= \int_0^\infty x^s f(s,x) \,\mathrm{d}x \\ &= \int_0^\infty x^s \mathrm{e}^{-x}(1+ax) \,\mathrm{d}x \\ &= \int_0^\infty x^s \mathrm{e}^{-x} \,\mathrm{d}x + a \int_0^\infty x^{s+1} \mathrm{e}^{-x} \\ &= \Gamma(s+1) + a (s+1) \Gamma(s+1) \\ &= (1 + a(s+1)) \Gamma(s+1) \\ \Longrightarrow a &= -\frac{1}{s+1} \\ \Longrightarrow f(s,x) &= \mathrm{e}^{-x}\left( 1 - \frac{x}{s+1} \right) \end{align}

Thus \begin{align} \zeta(-s) &= \lim_{\varepsilon \rightarrow 0^+} \lim_{m \rightarrow \infty} \sum_{n=1}^m n^s f(s, n \varepsilon) \\ &\overset{\star}{=} \lim_{m \rightarrow \infty} \lim_{\varepsilon \rightarrow 0^+} \sum_{n=1}^m n^s f(s, n \varepsilon) \\ &= \lim_{m \rightarrow \infty} \sum_{n=1}^m n^s f(s, 0) \\ &= \lim_{m \rightarrow \infty} \sum_{n=1}^m n^s \\ &= \sum_{n=1}^\infty n^s \end{align}

where the star indicates the non-rigorous step of exchanging limits. In fact, this regulator seems to work for all $$s \neq -1$$ (see some plots for $$s < -1$$ below). Hence we can also regularize the harmonic series as follows: \begin{align} \gamma &= \frac{1}{2} \lim_{\varepsilon \rightarrow 0^+} (\zeta(1+\varepsilon) + \zeta(1-\varepsilon)) \\ &= \frac{1}{2} \lim_{\varepsilon \rightarrow 0^+} \left(\lim_{m \rightarrow \infty} \sum_{n=1}^m n^{-1-\varepsilon} f(-1-\varepsilon, n \varepsilon) + \lim_{m \rightarrow \infty} \sum_{n=1}^m n^{-1+\varepsilon} f(-1+\varepsilon, n \varepsilon) \right) \\ &= \frac{1}{2} \lim_{\varepsilon \rightarrow 0^+} \lim_{m \rightarrow \infty} \left(\sum_{n=1}^m n^{-1-\varepsilon} f(-1-\varepsilon, n \varepsilon) + \sum_{n=1}^m n^{-1+\varepsilon} f(-1+\varepsilon, n \varepsilon) \right) \\ &\overset{\star}{=} \frac{1}{2} \lim_{m \rightarrow \infty} \lim_{\varepsilon \rightarrow 0^+} \left(\sum_{n=1}^m n^{-1-\varepsilon} f(-1-\varepsilon, n \varepsilon) + \sum_{n=1}^m n^{-1+\varepsilon} f(-1+\varepsilon, n \varepsilon) \right) \\ &= \frac{1}{2} \lim_{m \rightarrow \infty} \left(\sum_{n=1}^m n^{-1} f(-1, 0) + \sum_{n=1}^m n^{-1} f(-1, 0) \right) \\ &= \frac{1}{2} \lim_{m \rightarrow \infty} \left(\sum_{n=1}^m n^{-1} + \sum_{n=1}^m n^{-1} \right) \\ &= \lim_{m \rightarrow \infty} \sum_{n=1}^m n^{-1} \\ &= \sum_{n=1}^\infty n^{-1} \\ \end{align}

Here is the Mathematica code and its plots:

f[s_, x_] := Exp[-x] (1 + a x)
g[s_, t_] :=
Evaluate@Simplify[
f[s, t] /. Solve[Integrate[x^s f[s, x], {x, 0, Infinity}] == 0, a],
Assumptions -> Re[s] > -1]
g[s, t]
Table[{s,
Plot[{Zeta[-s],
Sum[n^s g[s, n \[Epsilon]], {n, 1, 1000}]}, {\[Epsilon], 0, 1},
Evaluated -> True]}, {s, -4, 4, 1/2}] // TableForm // Quiet
Plot[{EulerGamma,
Sum[(n^(-1 + \[Epsilon]) g[-1 + \[Epsilon], n \[Epsilon]] +
n^(-1 - \[Epsilon]) g[-1 - \[Epsilon], n \[Epsilon]])/
2, {n, 1, 1000}]}, {\[Epsilon], 0, 1}, Evaluated -> True]


• The similar limits you're computing need more care because it looks like you're just trying to swap the limit with the sum, which is in fact wrong because you need some kind of uniform bound the terms whose sum also converges. Commented Jun 16, 2015 at 19:52
• Indeed $\lim_{x\to 0^-}\lim_{n\to\infty}\sum_{k=1}^nf(k,x)\neq\lim_{n\to\infty}\lim_{x\to 0^-}\sum_{k=1}^nf(k,x)$ in general, but I'm trying to find the regularized value of divergent series using the first approach. Commented Jun 17, 2015 at 2:30
• For what it's worth, write $\cos(k^3x)=e^{ik^3x}+e^{-ik^3x}{2}$. Multiply out the sum to get two sums with summands like $k^3\exp(k^3x(1\pm i))$. Now notice that if you instead look at the sum of $\exp(k^3a x)$ and take it's derivative in $x$ you'll get what you're interested in. I think this kind of a sum is related to Jacobi theta functions. Commented Jun 17, 2015 at 2:39
• If you haven't seen this post by Terry Tao, it's worth a look. His approach isn't identical to yours — he uses a compactly supported cutoff function in the place of your exponentials — but I think it rhymes. In particular, he also recovers negative zeta values as the constant terms in asymptotic expansions of his sums, which are dominated by some larger non-constant term — very much like your "further information" identity. Commented Jun 18, 2015 at 21:31
• No limit approach is needed. The sum is in fact infinite. The fact that the $\zeta$ function extends the map $s\mapsto\sum n^{-s}$ for $s>0$ or $s>1$, doesn't mean that we can say the same on the extension. This is an "abuse of continuation", just like how you often write $f(\infty)$ to mean the value of $\lim_{x\to\infty}f(x)$. If that makes you cringe, then the abuse of $\zeta$ is equally cringe-worthy. If neither makes you cringe, then perhaps it's time to reflect on a career in physics. :-) Commented Jun 19, 2015 at 8:29

In Terence Tao's blog post about relating divergent series to negative zeta values, he uses Euler-Maclaurin summation to show that if $\eta$ is smooth and compactly-supported, and $\eta(0)=1$, then $$\sum_{n=0}^\infty n^s \eta(n/N)=C_{\eta,s} N^{s+1} + \zeta(-s) + O(1/N)$$ for any positive integer $s$, where $C_{\eta,s}=\int_0^\infty x^s \eta(x) \, dx$. In fact, it's possible to extend his proof to the case where $\eta$ is a Schwartz function — see below for a proof.

In particular, consider $\eta(x)=e^{-x}\cos x$. After a change of variables we can write your sum in the above form, with this $\eta$. Tao discusses the connection between this asymptotic form and the analytic continuation of the zeta function; I haven't gone through that part of his blog post in detail, but assuming that it, too, could be made to work when $\eta$ was Schwartz-class, that would presumably answer your first question about the connection between your sum and the zeta-regularized form of the original divergent sum.

To answer your second question, we can look at the leading coefficient $C_{\eta,s}=\int_0^\infty x^s e^{-x} \cos x \, dx$. After a bunch of annoying but straightforward integration by parts, it's possible to compute that $\int_0^\infty x^s e^{-x} \cos x=0$ when $s \equiv 1 \pmod{4}$, and $\int_0^\infty x^s e^{-x} \sin x=0$ when $s \equiv 3 \pmod{4}$. That is, whenever $s \equiv 1 \pmod{4}$, the $\eta(x)$ you've chosen coincidentally causes the leading term of the asymptotic expansion to drop out, exposing the zeta value which is the constant term of the expansion.

You can expose $\zeta(-4k)$ in a similar fashion by taking $\eta(x)=e^{-x} \frac{\sin x}{x}$, though since it's zero, that's not particularly interesting. For example, here's $\zeta(-4)$. Similarly, you can expose $\zeta(-4k+2)$ by taking $\eta(x)=e^{-x}(\cos x + \sin x)$. The interesting question — to which I don't know the answer — would be whether there's a nice choice of $\eta$ that exposes $\zeta(-4k+1)$, since that's the other case where it's actually nonzero. But this provides at least a partial answer to your third question.

Also note that if we take $\eta(x)=e^{-x}$, Tao's formula gives your "further information" expansion, as $\int_0^\infty x^n e^{-x} \, dx = n!$.

To prove Tao's formula whenever $\eta$ is smooth and rapidly decreasing, we'll use the infinite-series version of the Euler-Maclaurin formula (as found, e.g., in these notes): if $f(t)$ and all its derivatives tend to zero as $t \to \infty$, then $$\sum_{n=0}^\infty f(n) = \int_{0}^\infty f(t) \, dt + \frac{1}{2} f(0) - \sum_{i=2}^k \frac{b_i}{i!} f^{(i-1)}(0)-\int_0^\infty \frac{B_k(\{1-t\})}{k!} f^{(k)}(t) \, dt$$ for any fixed integer $k$, where the $b_k$ are the Bernoulli numbers, the $B_k$ are the Bernoulli polynomials, and $\{x\}$ denotes the fractional part of $x$.

In our case, we set $f(t)=t^s \eta(t/N)$ and use the $(s+2)$nd-order expansion. Taking each term in the formula one-by-one, note that:

$$\int_0^\infty f(t) \, dt=\int_0^\infty t^s \eta(t/N) \, dt = N^{s+1}\int_0^\infty x^s \eta(x) \, dx=C_{\eta,s}N^{s+1}$$ after making the substitution $x=t/N$; $$\frac{1}{2}f(0)=0$$ is immediate; \begin{align} -\sum_{i=2}^{s+2} \frac{b_i}{i!} f^{(i-1)}(0)&=-\frac{b_{s+1}}{(s+1)!} s! \eta(0)-\frac{b_{s+2}}{(s+2)!}\left(\frac{(s+1)!}{N} \eta'(0)\right)\\&=-\frac{b_{s+1}}{s+1} + O(1/N)\\&=\zeta(-s) + O(1/N) \end{align} as the first $s-1$ derivatives of $f$ vanish, and \begin{align} -\int_0^\infty \frac{B_{s+2}(\{1-t\})}{(s+2)!} f^{(s+2)}(t) \, dt &= -\int_0^\infty \frac{B_{s+2}(\{1-t\})}{(s+2)!} \frac{d^{s+2}}{dt^{s+2}}\left(t^s \eta(t/N)\right) \\ &=-\int_0^\infty \frac{B_{s+2}(\{1-Nx\})}{(s+2)!} \frac{1}{N^{s+2}}\frac{d^{s+2}}{dx^{s+2}}\left(N^sx^s \eta(x)\right) N \, dx \\ &=-\frac{1}{N} \int_0^\infty \frac{B_{s+2}(\{1-Nx\})}{(s+2)!} \frac{d^{s+2}}{dx^{s+2}}\left(x^s \eta(x)\right) \, dx \end{align} where we again make the substitution $x=t/N$.

This last quantity is bounded in magnitude by $\frac{1}{N}\frac{M}{(s+2)!}\int_0^\infty \frac{d^{s+2}}{dx^{s+2}}\left(x^s \eta(x)\right) \, dx$, where $M$ is the maximum value of $B_{n+2}(x)$ on $[0,1]$. Thus it too is $O(1/N)$ so long as that integral converges.

So Tao's formula works, not only for compactly supported $\eta$, but for $\eta$ such that $f(t)=t^s \eta(t)$ as well as all of $f$'s derivatives are integrable on $[0,\infty)$. In particular, it will work for all $s$ whenever $\eta$ is Schwartz on $[0,\infty)$.

• This is excellent! The problem is therefore reduced to finding a Schwartz function $\eta$ such that $\int_0^\infty x^3\eta(x)\,\mathrm{d}x=0$ and $\lim_{x\to0}\eta(x)=1$, correct? Commented Jun 24, 2015 at 0:09
• If you're looking for a family of identities which is as nice as the ones you came up with for the other interesting residue class, you need $\eta$ such that $\int_0^\infty x^{3+4k} \eta(x) \, dx=0$ for all $k$. Or you could just content yourself with knowing that for any $k$, you can choose a different $\eta$ (in which case there's probably always some linear combination of, say, $e^{-x}$ and $e^{-2x}$ that'd do the trick)... Commented Jun 24, 2015 at 2:15
• I added a fractional regulator to my question based on your answer. Commented Jul 1, 2019 at 4:49

Explanation for the Behavior of the Series (Parts $2$ and $3$)

Note that \begin{align} \sum_{k=0}^\infty e^{ikx} &=\frac1{1-e^{ix}}\\ &=\frac12+\frac i2\cot\left(\frac x2\right)\\ &=\frac ix+\frac12-i\,\underbrace{\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\substack{\text{odd function of x}\\\text{analytic on (-\pi,\pi)}}}\tag{1} \end{align}

Odd Exponents

Taking $2n-1$ derivatives of $(1)$, we get $$\sum_{k=0}^\infty k^{2n-1}e^{ikx}=\frac{(-1)^n(2n-1)!}{x^{2n}}+(-1)^n\underbrace{\frac{\mathrm{d}^{2n-1}}{\mathrm{d}x^{2n-1}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\substack{\text{even function of x}\\\text{analytic on (-\pi,\pi)}}}\tag{2}$$

Exponents which are $\boldsymbol{1\bmod{4}}$

If $n$ is odd, so that $2n-1\equiv1\pmod{4}$, substitute $x=(1-i)y$ in $(2)$ to get $$\sum_{k=0}^\infty k^{2n-1}e^{k(1+i)y} =\underbrace{\frac{(2n-1)!}{(2i)^ny^{2n}}}_{\substack{\text{pure imaginary}\\\text{since n is odd}}} -\underbrace{\frac{\mathrm{d}^{2n-1}}{\mathrm{d}x^{2n-1}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\text{plug in x=(1-i)y}}\tag{3}$$ Thus, the real part of $(3)$ is $$\sum_{k=0}^\infty k^{2n-1}e^{ky}\cos(ky)=\left.-\frac{\mathrm{d}^{2n-1}}{\mathrm{d}x^{2n-1}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)\right|_{x=(1-i)y}\tag{4}$$ Taking the limit of $(4)$ as $y\to0$ yields the identity given in the question \begin{align} &\text{When n is odd; that is, }2n-1\equiv1\pmod{4}:\\ &\bbox[5px,border:2px solid #C0A000]{\lim_{y\to0^-}\sum_{k=0}^\infty k^{2n-1}e^{ky}\cos(ky)=\left.-\frac{\mathrm{d}^{2n-1}}{\mathrm{d}x^{2n-1}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)\right|_{x=0}}\tag{5} \end{align}

Exponents which are $\boldsymbol{3\bmod{4}}$

If $n$ is even, so that $2n-1\equiv3\pmod{4}$, the situation is a bit more complicated. We can still substitute $x=(1-i)y$ in $(2)$, getting $$\sum_{k=0}^\infty k^{2n-1}e^{k(1+i)y} =\underbrace{\frac{(2n-1)!}{(2i)^ny^{2n}}}_{\text{real since n is even}} +\underbrace{\frac{\mathrm{d}^{2n-1}}{\mathrm{d}x^{2n-1}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\text{plug in x=(1-i)y}}\tag{6}$$ We can also substitute $x=-iy$ in $(2)$, getting $$\sum_{k=0}^\infty k^{2n-1}e^{ky}=\frac{(2n-1)!}{y^{2n}}+\underbrace{\frac{\mathrm{d}^{2n-1}}{\mathrm{d}x^{2n-1}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\text{plug in x=-iy}}\tag{7}$$ Taking a linear combination of $(6)$ and $(7)$ to cancel the singular part and maintain the constant part, we get \begin{align} &\text{When n is even; that is, }2n-1\equiv3\pmod{4}:\\ &\bbox[5px,border:2px solid #C0A000]{\lim_{y\to0^-}\sum_{k=0}^\infty k^{2n-1}e^{ky}\frac{2^n\cos(ky)-(-1)^{n/2}}{2^n-(-1)^{n/2}} =\left.\frac{\mathrm{d}^{2n-1}}{\mathrm{d}x^{2n-1}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)\right|_{x=0}}\tag{8} \end{align}

Even Exponents

Taking $2n$ derivatives of $(1)$, we get $$\sum_{k=0}^\infty k^{2n}e^{ikx} =(-1)^n\frac{(2n)!}{x^{2n+1}}i -(-1)^ni\underbrace{\frac{\mathrm{d}^{2n}}{\mathrm{d}x^{2n}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\substack{\text{odd function of x}\\\text{analytic on (-\pi,\pi)}}}\tag{9}$$ Substitute $x=(1-i)y$ in $(9)$ to get $$\sum_{k=0}^\infty k^{2n}e^{k(1+i)y} =-\frac{(2n)!}{y^{2n+1}}\frac{1+i}{(2i)^{n+1}} -(-1)^ni\underbrace{\frac{\mathrm{d}^{2n}}{\mathrm{d}x^{2n}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\text{plug in x=(1-i)y}}\tag{10}$$ We can also substitute $x=-iy$ in $(9)$ to get $$\sum_{k=0}^\infty k^{2n}e^{ky} =-\frac{(2n)!}{y^{2n+1}} -(-1)^ni\underbrace{\frac{\mathrm{d}^{2n}}{\mathrm{d}x^{2n}}\left(\frac1x-\frac12\cot\left(\frac x2\right)\right)}_{\text{plug in x=-iy}}\tag{11}$$ Taking a linear combination of $(10)$ and $(11)$ to cancel the singular part and maintain the constant part, we get $$\bbox[5px,border:2px solid #C0A000]{\lim_{y\to0^-}\sum_{k=0}^\infty k^{2n}e^{ky}\frac{2^{n+1}\cos(ky)+(-1)^{\lfloor(n-1)/2\rfloor}}{2^{n+1}+(-1)^{\lfloor(n-1)/2\rfloor}}=0}\tag{12}$$ where $(-1)^{\lfloor(n-1)/2\rfloor}=\mathrm{Re}\left(\frac{i-1}{i^n}\right)$.

Unified Expression

Combining $(5)$, $(8)$, and $(12)$ yields $$\bbox[5px,border:2px solid #C0A000]{\lim_{x\to0^-}\sum_{k=0}^\infty k^ne^{kx}\frac{2^{\lceil n/2\rceil}\cos(kx)+c_n}{2^{\lceil n/2\rceil}+c_n}=\zeta(-n)}\tag{13}$$ where $c_n$ is given by $$c_n=-\tfrac14(-1)^{\lceil n/4\rceil}\left(2+(-1)^{\lceil n/2\rceil}\left(1-(-1)^n\right)\right)$$ or by $$c_n=-(-1)^{\lceil n/4\rceil}\sin^2\left(\tfrac\pi4(n-1)\right)$$ or in the following table $$\small\begin{array}{c|r|r|r|r|r|r|r} n\bmod8&1&2&3&4&5&6&7&8\\ \hline c_n&\hphantom{-}0&\hphantom{-}\frac12&\hphantom{-}1&\hphantom{-}\frac12&\hphantom{-}0&-\frac12&-1&-\frac12 \end{array}\tag{14}$$

• Could this regulator be generalized to fractional $n$? Commented Jun 30, 2019 at 0:43
• There might be, but I am not sure how to approach that.
– robjohn
Commented Jun 30, 2019 at 4:05
• I added a regulator that works for fractional exponents to my question. Commented Jul 1, 2019 at 4:38

The infamous $-\frac{1}{12}$ pops up as the coefficient of of a range of expansions, such as

$\dfrac{z}{{\mathrm e}^{-z}-1}=-1-\dfrac{1}{2}z-\dfrac{1}{12}z^2+{\mathcal O}(z^3)$

like in the Todd class or similarly in the Baker–Campbell–Hausdorff formula and so on. In your post, you are led to it by ${\mathrm e}^{kx}\cos(kx)$, which after the reparametrization $x\mapsto \log(z)$ I'd write as as $z^k\cos(\log(z^k))$. But in any case, the higher order expansion terms of the cosine make this a not so nice choice and you get a lot of cases where the limit procedure doesn't work. I've worked out some broader perspectives..

The analytical continuation, e.g. involved in defining the Zeta function*

$\zeta(s) = \dfrac{1}{\Gamma(s)} \int_{0}^{\infty} \dfrac{x ^ {s-1}}{{\mathrm e} ^ x - 1}\mathrm{d}x,$

is a process which smoothly probes local values of the functional expressions $f$ in $\sum_{k=0}^\infty f(k)$. You can formulate the regularization of the sum in a way that reflects that, e.g. by using a local mean like

$\langle f(k)\rangle:=\int_{k}^{k+1}f(k')\,{\mathrm d}k'$.

So the limit $\lim_{z\to 1}$ of the sum

$0+1\,z^1+2\,z^2+3\,z^3+\dots$

diverges, because of the pole in

$\sum_{k=0}^\infty k\,z^k=z\dfrac{{\mathrm d}}{{\mathrm d}z}\sum_{k=0}^\infty z^k=z\dfrac{{\mathrm d}}{{\mathrm d}z}\dfrac{1}{1-z}=\dfrac{z}{(z-1)^2}, \hspace{1cm} z\in(0,1)$

So let's consider the sum of smooth deviations. With

$\langle k\,z^k\rangle=z\dfrac{{\mathrm d}}{{\mathrm d}z}\langle z^k\rangle=z\dfrac{{\mathrm d}}{{\mathrm d}z}\langle {\mathrm e}^{k \log(z)}\rangle=\left.z\dfrac{{\mathrm d}}{{\mathrm d}z}\dfrac{z^{k'}}{\log(z)}\right|_{k}^{k+1}$.

we find the sum $\sum_{k=0}^n\langle k\,z^k\rangle$ involves canceling upper and lower bounds and we're left with $\dfrac{z^0}{\log(z)^2}$ plus terms suppressed by $z^n$. Finally, using the expansion

$r^2\dfrac{1}{\log(1+r)^2}=\dfrac{1}{1-r+\left(1-\frac{1}{1!\,2!\,3!}\right)r^2+{\mathrm{O}}(r^3)}=1+r+\dfrac{1}{1!\,2!\,3!}r^2+{\mathrm{O}}(r^3),$

we find

$\sum_{k=0}^\infty \left(k\,z^k-\langle k\,z^k\rangle\right)=\dfrac{z}{(z-1)^2}-\dfrac{1}{\log(z)^2}=-\dfrac{1}{12}+{\mathcal O}\left((z-1)^1\right).$

The picture shows the two functions $\dfrac{z}{(z-1)^2}$ and $\dfrac{1}{\log(z)^2}$, as well as their difference (blue, red, yellow). While the functions themselves clearly have a pole at $z=1$, their difference converges against

$$-\frac{1}{1!\,2!\,3!}=-\frac{1}{12}=-0.08{\dot 3}.$$

And here's the Mathematica calculation for higher $n$, as well as the corresponding values of the Zeta function (which match)

I've tried to generalize those in various directions and stumbled upon some futher relations. For example, when comparing finite differences to their first order approximations, we find

$\dfrac{f'(x)\,h}{f(x+h)-f(x)}=1-\dfrac{f''(x)}{2!}\left(\dfrac{h}{f'(x)}\right)+\left(\dfrac{f''(x)\,f''(x)}{2!\,2!}-\dfrac{f'(x)\,f'''(x)}{1!\,3!}\right)\left(\dfrac{h}{f'(x)}\right)^2+{\mathcal O}(h^3).$

Note that $\dfrac{1}{2!\,2!}-\dfrac{1}{1!\,3!}=\dfrac{1}{2!\,3!}(3-2)=\dfrac{1}{12}$, which is the coeffient $\dfrac{1}{2!}B_2$ in the related MacLaurin formula

And the log-substraction scheme can be performed more generally for higher order poles $\frac{1}{(z-1)^n}$ using

$\dfrac{1}{\log(z)^n}=\dfrac{1}{(z-1)^n}\left(1+\frac{n}{2}(z-1)+\frac{n}{2}\frac{3n-5}{12}(z-1)^2+\frac{n}{2}\frac{(n-2)(n-3)}{24}(z-1)^3+\dots\right)$

Here, e.g. plug in $n=2$ to get $\dfrac{n}{2}\dfrac{3n-5}{12}(z-1)^2=\dfrac{1}{12}$.

*As I see from your profile you're interested a physics (I'm a physicist myself) I like to point out how the body of the integrand in that Zeta function representation is the expression of Planck's law :). The exponential function ${\mathrm e}^x\cdot{\mathrm e}^y={\mathrm e}^{x+y}$, say in statistical physics where $-\frac{1}{12}$ come to play a role, enters here because the probabilities of divided systems behave multiplicative, while the energy, by definition, behaves additive.