Constant of integration of exponential-type function

Question

The Laurent series about $0$ of $1/{z(e^z-1)}$ is

$$\frac{1}{z(e^z-1)}=\frac{1}{z^2}-\frac{1}{2z}+\sum_{n=2}^\infty\frac{B_n}{n!}z^{n-2}$$

When we integrate this we have

$$\int\frac{1}{z(e^z-1)}dz=-\frac{1}{z}-\frac{\ln{|z|}}{2}+\sum_{n=2}^\infty\frac{B_n}{n!(n-1)}z^{n-1}+C$$

Something interesting happens when we subtract the expansion from the definite integral form. It equals a constant.

$$C_0=\int_\infty^x\frac{1}{z(e^z-1)}dz-\left(-\frac{1}{x}-\frac{\ln{|x|}}{2}+\sum_{n=2}^\infty\frac{B_n}{n!(n-1)}x^{n-1}\right)$$

Unfortunately the series is only defined on $0\lt x\lt 2\pi$.

What is the value of the constant $C_0$?

I thought about seeing if I can do the expansion about some other point $a$, so here's what I did

$$\begin{align} \frac{1}{z(e^z-1)}&=\frac{1}{a(z/a)(e^{a(z/a)}-1)}\\ &=\sum_{n=0}^\infty\frac{B_n}{n!}a^{n-2}\left(\frac{z}{a}\right)^{n-2}\\ &=\sum_{n=0}^\infty\frac{B_n}{n!}a^{n-k-2}\sum_{k=0}^\infty{n-2\choose{k}}\left(z-a\right)^k\\ &=\sum_{k=0}^\infty c_k\left(z-a\right)^k \end{align}$$

Where

$$c_k=\sum_{n=0}^\infty\frac{B_n}{n!}{n-2\choose{k}}a^{n-k-2}$$

Using the Euler-Maclaurin summation formula with $f(x)=\frac{(-1)^k}{k!}\Gamma(k+1,-x)$ we get

$$c_k=\frac{\left(-1\right)^k}{a^{k+1}k!}\sum_{n=1}^\infty \Gamma\left(k+1,na\right)$$

So by integration we have

$$\int\frac{1}{z(e^z-1)}dz=\sum_{k=0}^\infty \frac{c_k}{k+1}\left(z-a\right)^{k+1}+C$$

Then

$$\int_\infty^x\frac{1}{z(e^z-1)}dz=\lim\limits_{a\to\infty}\sum_{k=0}^\infty \frac{c_k}{k+1}\left(x-a\right)^{k+1}$$

Then we can try subtracting the original series expression that is defined on $(0,2\pi)$ from this expression that is defined on $(0,\infty)$ for some $x\in (0,2\pi)$.

Answer 1

The answer is $\color{blue}{C_0=(\ln2\pi-\gamma)/2}$. We rewrite the definition as \begin{align*} C_0&=\lim_{x\to0^+}\left[\frac1x+\frac{\ln x}2-\int_x^\infty\frac{dz}{z(e^z-1)}\right] \\&=\lim_{x\to0^+}\left[\frac{\ln x}2+\int_x^\infty\left(\frac1z-\frac1{e^z-1}\right)\frac{dz}z\right] \end{align*} and use $\lim\limits_{x\to0^+}\left(\ln x+\int_x^\infty(e^{-z}/z)\,dz\right)=-\gamma$ (see e.g. here) to get $$ C_0=-\frac\gamma2-\color{blue}{\int_0^\infty\left(\frac{e^{-z}}2-\frac1z+\frac1{e^z-1}\right)\frac{dz}z}. $$ The integral looks close to Binet's log-gamma integrals, and can be found the way these are obtained.

Another approach is a link to Riemann's $\zeta$; consider, for $\Re a>0$ and $\Re s>0$, $$ I(s,a)=\int_0^\infty\left(\frac{e^{-z}}2-\frac{e^{-az}}{z}+\frac1{e^z-1}\right)z^{s-1}\,dz. $$ For $\Re s>1$, this is easily integrated termwise, giving $$ I(s,a)=\Gamma(s)\left(\frac12-\frac{a^{1-s}}{s-1}+\zeta(s)\right), $$ which then holds for $\Re s>0$ by analytic continuation.

Now, for $0<\Re s<1$, we can take $a\to0^+$, which gives $$ \int_0^\infty\left(\frac{e^{-z}}2-\frac1z+\frac1{e^z-1}\right)z^{s-1}\,dz=\Gamma(s)\left(\frac12+\zeta(s)\right) $$ (in fact, this holds for $-1<\Re s<1$, again by a.c.).

The integral we're looking for is obtained at $s\to 0$, and equals $\zeta'(0)=-\frac12\ln2\pi$.