Evaluating $\int_0^y \frac{z^2e^{az}}{e^z-1}\,\mathrm dz$.

Question

Given $a,y\in\mathbb{C}$, I want to find a closed form for $$ \int_0^y \frac{z^2e^{az}}{e^z-1}\,\mathrm dz. $$

Background:

This integral was obtained from an infinite sum involving Bernoulli number and I want to simplify it. I noticed using wxMaxima software that it has a closed form for $a=1$ and $y\in\mathbb{R}$ in terms of the polylogarithm, but I have not been able to find such an expression for another $a$.

Such integral comes from the following: For each $n\in\mathbb{N}$ let $B_n$ denote the $n$-th Bernoulli numbers, which can be defined as $$\frac{x}{e^x-1} = \sum_{n\in\mathbb{N}} \frac{B_n}{n!}(x).$$ For some $x,w\in\mathbb{C}$ I want to take the limit $$\lim_{\Re(w)\to\infty}e^w\sum_{n\in\mathbb{N}} \frac{B_n}{n!}(-x)\sum_{k=0}^n \frac{w^k}{k!}.$$ Substituting a known formula I found on Uniform Convergence Behavior of the Bernoulli Polynomials, John Mangual page 2 $$\sum_{k=0}^n \frac{w^k}{k!} = \int_0^w z^{n+1}e^z dz,$$ related to a Gábor Szegő's paper in 1928, then gives $$ \sum_{n\in\mathbb{N}} \frac{B_n}{n!}(-x)\sum_{k=0}^n \frac{w^k}{k!} = \sum_{n\in\mathbb{N}} \frac{B_n}{n!}(-x) \int_0^w z^{n+1}e^z dz = \int_0^w \sum_{n\in\mathbb{N}} z \frac{B_n}{n!}(-xz)^{n}e^z dz \\ = \int_0^w z\frac{-xz}{e^{-xz}-1} e^z dz $$ Finally the main integral comes after the change of variable $z' = \frac{z}{-x}.$

Answer 1

We have $$ I(a,y)=\int_0^y \frac{z^2e^{az}}{e^z-1}\,\mathrm dz. $$ We will first perform integration by parts with $u=z^2$ and $\mathrm dv=e^{az}/(e^z-1)\,\mathrm dz$. Clearly, $\mathrm du=2 z\,\mathrm dz$ so we write $$ \begin{align} v &=\int \frac{e^{az}}{e^z-1}\,\mathrm dz\\ &\overset{x=e^z}{=}-\int\frac{x^{a-1}}{1-x}\,\mathrm dx\\ &=-\int_0^x\frac{w^{a-1}}{1-w}\,\mathrm dw\\ &\overset{w=tx}{=}-x^a\int_0^1\frac{t^{a-1}}{1-xt}\,\mathrm dt. \end{align} $$ Comparing with DLMF 15.6.1 and reintroducing $x=e^z$ yields $$ v=-\frac{1}{a}e^{az}F(1,a;1+a;e^z), $$ such that $$ I(a,y)=-\frac{1}{a}z^2e^{az}F(1,a;1+a;e^z)\bigg|_{z=0}^y+\frac{2}{a}\underbrace{\int_0^yze^{az}F(1,a;1+a;e^z)\,\mathrm dz}_{J(a,y)}, $$ where $F(\alpha,\beta;\gamma;z)$ is the Gauss hypergeometric function. The limit as $z\to 0$ is zero which is due to $F(1,a;1+a;e^z)\sim-a\log(1-e^z)$, c.f. DLMF 15.4.21, and $\lim_{z\to 0} z\log(1-e^z)=0$. Hence, $$ I(a,y)=-\frac{1}{a}y^2e^{ay}F(1,a;1+a;e^y)+\frac{2}{a}J(a,y). $$

Our next task is to determine $J(a,y)$. As before we will use integration by parts with $u=z$ and $\mathrm dv=e^{az}F(1,a;1+a;e^z)\,\mathrm dz$ so that $$ \begin{align} v &=\int e^{az}F(1,a;1+a;e^z)\,\mathrm dz\\ &\overset{x=e^z}{=}\int x^{a-1}F(1,a;1+a;x)\,\mathrm dx\\ &=\int_0^x w^{a-1}F(1,a;1+a;w)\,\mathrm dw\\ &\overset{w=xt}{=}x^a\int_0^1 t^{a-1}F(1,a;1+a;xt)\,\mathrm dt. \end{align} $$ Comparing with DLMF 16.5.2 and again reintroducing $x=e^z$ gives $$ v=\frac{1}{a}e^{az}{_3F_2}\left({1,a,a\atop 1+a,1+a};e^z\right), $$ and $$ J(a,y)=\frac{1}{a}ze^{az}{_3F_2}\left({1,a,a\atop 1+a,1+a};e^z\right)\bigg|_{z=0}^y-\frac{1}{a}\int_0^ye^{az}{_3F_2}\left({1,a,a\atop 1+a,1+a};e^z\right)\,\mathrm dz, $$ where ${_pF_q}(\mathbf a;\mathbf b;z)$ is the generalized hypergeometric function. The limit at zero presents no problems this time and so we evaluate the limit term by direct substitution. With the help of this identity we then write the ${_3F_2}(\cdot)$ function in terms of the Lerch transcendent, $\Phi(z,s,a)$, and subsequently obtain $$ J(a,y)=aye^{ay}\Phi(e^y,2,a)-a\int_0^ye^{az}\Phi(e^z,2,a). $$ To tidy things up a bit we substitute the current form of $J(a,y)$ into $I(a,y)$ and write $$ I(a,y)=-\frac{1}{a}y^2e^{ay}F(1,a;1+a;e^y)+2ye^{ay}\Phi(e^y,2,a)-2\underbrace{\int_0^ye^{az}\Phi(e^z,2,a)}_{J^\ast(a,y)}. $$ Our last task is to find $J^\ast(a,y)$. Substituting $t=e^z$ gives $$ J^\ast(a,y)=\int_1^{e^y}t^{a-1}\Phi(t,2,a)\,\mathrm dt. $$ The antiderivative of the integral is given here yielding $$ \begin{align} J^\ast(a,y) &=t^a\sum_{k=0}^\infty\frac{t^k}{(k+a)^3}\bigg|_{t=1}^{e^y}\\ &=t^a\Phi(t,3,a)\bigg|_{t=1}^{e^y}\\ &=e^{ay}\Phi(e^y,3,a)-\Phi(1,3,a)\\ &=e^{ay}\Phi(e^y,3,a)-\zeta(3,a), \end{align} $$ where $\zeta(s,a)$ is the Hurwitz zeta function. Substituting $J^\ast$ into $I$ then gives the final result $$ \bbox[5px,border:2px solid #C0A000]{% I(a,y)=2\zeta(3,a)+2e^{ay}(y\Phi(e^y,2,a)-\Phi(e^y,3,a))-\frac{1}{a}y^2e^{ay}F(1,a;1+a;e^y).% } $$

Here is the Mathematica code for ploting the result:

Int[a_, y_] := 
  2 HurwitzZeta[3, a] + 
   2 E^(a y) (y HurwitzLerchPhi[E^y, 2, a] - 
      HurwitzLerchPhi[E^y, 3, a]) - 
   1/a y^2 E^(a y) Hypergeometric2F1[1, a, 1 + a, E^y];
Plot[{Int[1, y], Int[2, y], Int[3, y]}, {y, 0, 2}, PlotRange -> All, 
 PlotLegends -> "Expressions"]

...and the plot it produces: