Closed-form of $\int_0^\infty \frac{t^s}{(e^t-1)^z}dt$ I am looking for a closed form for the integral $$\int_0^\infty \frac{t^s}{(e^t-1)^z}dt$$ valid for $s,z$ being both complex numbers, hopefully using complex analysis. I have already evaluated this integral when $s$ is complex and $z$ is a positive integer. In that case, the result is $$\frac{\Gamma(s+1)}{\Gamma(z)}\sum_{k=1}^z s(z,k)\zeta(s-k+2),$$ where the coefficients $s(z,k)$ are the Stirling numbers of the first kind.
Edit: as I said in the comments, I also found the closed form using
the generalized hypergeometric function $$\frac{\Gamma(s+1)}{z^{s+1}}\,_{s+2}F_{s+1}(-z,z,z,\dotsc;z+1,z+1,\dotsc;-1).$$ Unfortunately, this is an extension for $z$ but not for $s$.
 A: Actualisation to solution for complex z and s.
By using theorem of my big ego (Integral representation for series of any order)
$\displaystyle   \sum\limits^z_{k=x}     \;\!\!\;\!  \!\!\!\!\!\!\!\!\! \lower -0.2pt {\infty}  \quad \!\!\!     f (k)  =\overbrace {
\sum_{k_{z-1}=x} ^{\infty} ... \sum_{k_1=k_2}^{\infty} \sum_{k_0=k_1}^{\infty}}^{z}f (k_0)
  = \sum_{k=0}^{\infty}  \binom{k+z-1}{z-1}   f (k+x)= \frac {i   }{2}\int_{-\frac {1}{2}-i\infty}^{-\frac {1}{2}+i\infty}  \binom{t+z-1}{z-1}  \cot(\pi t) 
  f (t+x)   dt.$
This is simplified version of theorem of my big ego, becouse yours integral gives no results for divergent series. Now by using of this theorem you can write
$\displaystyle \int_0^\infty \frac{e^{-zt}t^s}{(1-e^{-t})^z}dt =\int_0^\infty \sum\limits^z_{k=z}     \;\!\!\;\!  \!\!\!\!\!\!\!\!\! \lower -0.2pt {\infty}  \quad \!\!\!    e^{-kt}t^s dt  = \Gamma (s+1)   \sum\limits^z_{k=z}     \;\!\!\;\!  \!\!\!\!\!\!\!\!\! \lower -0.2pt {\infty}  \quad \!\!\!    \frac {1}{k^{s+1}}  =   \frac { \Gamma (s+1) i }{2}\int_{-\frac {1}{2}-i\infty}^{-\frac {1}{2}+i\infty}  
  \frac { \binom{t+z-1}{z-1}  \cot(\pi t)  }{(t+z )^{s+1}}  dt.$
All you have to do is write $\frac{(t+z-1)!}{t!} $ as Laurent series around $-z$. What I should do at first place.
$\displaystyle \frac{(t+z-1)!}{t!} = \frac { (t + z)^{-1}   }{Γ(1 - z)}-\frac { \psi^0( 1 - z) + \gamma }{Γ(1 - z)} + \frac {(t + z) (\frac { [\psi^0(1 - z)]^2 }{2}+  \gamma \psi^0( 1 - z) - \frac{ \psi^1(1 - z) }{2}+ \frac {π^2}{12} + \frac {\gamma ^2}{2})}{Γ(1 - z)} -... $
Notice that $\psi^m (1-z)=(-1)^{m+1}m! \zeta (m+1,1-z) $ and  $\frac { i }{2}\int_{-\frac {1}{2}-i\infty}^{-\frac {1}{2}+i\infty}  
  \frac{  \cot(\pi t)  }{(t+z )^{s}}  dt=\zeta (s,z)$, so representation of $ \int_0^\infty \frac{t^s}{(e^t-1)^z}dt$ you can write as series in terms of Hurtwiz zeta function and Gamma function. For positive integer $z $ the equation is reducable to form you were tolking about (that one with Striling numbers of the first kind). Finaly you get
$\displaystyle  \int_0^\infty \frac{t^s}{(e^t-1)^z}dt = \frac {sin (\pi z)\Gamma (s+1)}{\pi}\left [\zeta (s+2,z) -\zeta (s+1,z)(\psi^0( 1 - z) + \gamma)... \right]$
I used reflection formula for Gamma function just in case. To be honest I can't count integral over closed line (used in Laurent series) so I will not write it as series form but it for sure exist. As you wanted there are only 'zeta function' type functions. Imo, your finte sum is secondary to my solution so I would treat it in the same way. My work here is done, rest is up to you Samurai
A: Here's my take on this integral.
$$
I(s, z) = \int_0^\infty t^s (e^t - 1)^{-z}dt
$$
To start off, we'll use the following substitution.
$$u = e^{-t}\\ t = -\ln u \\ dt = -u^{-1} du$$
Which then we'll have
$$
-\int_1^0 (-\ln u)^s (u^{-1} - 1)^{-z} u^{-1}du \\
= (-1)^s \int_0^1 \ln^s u (u^{-1} - 1)^{-z} u^{-1} du \\
= (-1)^s \int_0^1 \ln^s u [u^{-1}(1 - u)]^{-z} u^{-1} du \\
= (-1)^s \int_0^1 \ln^s u (1 - u)^{-z} u^{z-1} du
$$
Now let's take a look at something similar, namely the Beta function, which is defined as the following integral.
$$
\text{B}(x, y) = \int_0^1 t^{x-1} (1-t)^{y-1} dt
$$
What I noticed is that if we take the function's partial derivative with respect to $x$, we'll get
$$
\partial_x \text{B}(x, y) = \partial_x \int_0^1 t^{x-1} (1-t)^{y-1} dt
$$
Using the Leibniz integral rule, we'll get
$$
\partial_x \text{B}(x, y) = \int_0^1 \partial_x t^{x-1} (1-t)^{y-1} dt = \int_0^1 t^{x-1} \ln t (1-t)^{y-1} dt
$$
Similarly, if we take the second order derivative, we'll get
$$
\partial_x^2 \text{B}(x, y) = \int_0^1 t^{x-1} \ln^2 t (1-t)^{y-1} dt
$$
And if we continue this to the $s$-th order derivative, we'll get
$$
\partial_x^s \text{B}(x, y) = \int_0^1 t^{x-1} \ln^s t (1-t)^{y-1} dt
$$
Plug in $x=z$, $y=1-z$ and multiply the whole thing with $(-1)^s$, we'll get back our original integral.
$$
(-1)^s \partial_x^s \text{B}(x, 1-z) \Bigg|_{x = z} = (-1)^s \int_0^1 t^{z-1} \ln^s t (1-t)^{-z} dt = I(s, z)
$$
And for simplicity in the future, I'll use the result from this to turn the Beta function into a series for the sake of ease in the future.
$$
(-1)^s \partial_x^s \sum_{k=0}^\infty \binom{k+z-1}{k} (x+k)^{-1} \Bigg|_{x = z}
$$
The condition for this function (which is also our integral) to converge is for
$$
\begin{cases}
\Re(x) > 0 \\ 
\Re(y) > 0
\end{cases}
$$
Which means we'll have
$$
\begin{cases}
\Re(z) > 0 \\ 
\Re(1-z) > 0
\end{cases}
\Rightarrow
\begin{cases}
\Re(z) > 0 \\ 
\Re(z) < 1
\end{cases}
\Rightarrow 0 < \Re(z) < 1
$$
For the $s$-th order derivative, I noticed that
$$
\partial_x (x+k)^{-1} = -1!(x+k)^{-2} \\
\partial_x^2 (x+k)^{-1} = 2!(x+k)^{-3} \\
\partial_x^3 (x+k)^{-1} = -3!(x+k)^{-4} \\
\partial_x^4 (x+k)^{-1} = 4!(x+k)^{-5}
$$
And if we continue this to the $s$-th order, we'd get
$$
\partial_x^s (x+k)^{-1} = (-1)^s s! (x+k)^{-(s + 1)}
$$
But this only works for natural $s$. However, there's a way around that. By using the Gamma function, we can extend the formula to the real numbers.
$$
\partial_x^s (x+k)^{-1} = (-1)^s \Gamma(s+1) (x+k)^{-(s + 1)} \\
$$
Plugging this back into our original integral, we'll get
$$
(-1)^s \sum_{k=0}^\infty \binom{k+z-1}{k} (-1)^s \Gamma(s+1) (x+k)^{-(s+1)} \Bigg|_{x = z} \\
= (-1)^{2s} \Gamma(s+1) \sum_{k=0}^\infty \binom{k+z-1}{k} (z+k)^{-(s+1)}
$$
There's one problem though, this whole thing diverges for all negative integer $s$. I tried some other ways to evaluate for that case but nothing worked, so we'll just ignore that case in the mean time. I may or may not update this answer.
Now, if you just stopped there, you'd have a good answer for real $s$. But since you asked for complex $s$, my idea is simply to just use the same answer, but instead for real $s$ only, you also use it for complex $s$
Which means our final result for real $s$ is
$$
I(s, z) = (-1)^{2s} \Gamma(s+1) \sum_{k=0}^\infty \binom{k+z-1}{k} (z+k)^{-(s + 1)} \\
\text{for $s \in \mathbb{R}\backslash\mathbb{Z}^-$, $z \in \mathbb{C}$ and $0 < \Re(z) < 1$}
$$
And our final result for complex $s$ is
$$
I(s, z) = (-1)^{2s} \Gamma(s+1) \sum_{k=0}^\infty \binom{k+z-1}{k} (z+k)^{-(s + 1)} \\
\text{for $s \in \mathbb{C}$, $\Re(s) \not\in \mathbb{Z}^-$, $z \in \mathbb{C}$ and $0 < \Re(z) < 1$}
$$
