Sum=Integral? Find $m$ so that $\sum _{n=1}^{\infty } \frac{n^m}{e^{2 \pi n}-1}=\int_0^{\infty } \frac{x^m}{e^{2 \pi x}-1} \, dx$

The standard method is to use that $$\frac{\pi^2}{\sin^2(\pi z)}-\sum_m \frac1{(z+m)^2}=0$$ (it is a $1$-periodic entire function which is bounded on $\Re(z)\in [0,1]$ and which vanishes at $i\infty$)

to get that for $k\ge 4$ even and $\Im(z)>0$ $$G_k(z)=\sum_{(n,m)\ne (0,0)} \frac1{(zn+m)^k}= 2\zeta(k)+\frac{4i\pi}{(k-1)!}\sum_{n\ge 1}\sum_{d\ge 1}(2i\pi d)^{k-1}e^{2i\pi dnz}$$

When $k\equiv 2\bmod 4$ we get that $G_k(i)=0$.

On the other hand $$\int_0^\infty \frac{x^{k-1}}{e^{2\pi x}-1}dx=(k-1)!(2\pi)^{-k} \zeta(k),\qquad\sum_{d\ge 1}\frac{d^{k-1}}{e^{2\pi d}-1}=\sum_{d\ge 1}\sum_{n\ge 1}d^{k-1}e^{-2\pi dn}$$

The case $k\equiv 0\bmod 4$ has a solution too. From the theory of modular forms (mainly that $E_4(z)^3-E_6(z)^2$ is a non-vanishing cusp form), with $E_k(z)=\frac{G_k(z)}{2\zeta(k)}$: $$E_k(z)=\sum_{4a+6b=k} c_{k,a,b}E_4(z)^a E_6(z)^b$$ for some rational coefficients $c_{k,a,b}$. As $E_6(i)=0$ we get $$E_k(i) = c_{k,k/4,0}E_4(i)^{k/4}$$ The known closed-form of $E_4(i)$ in term of $\Gamma(1/4)$ implies that $E_k(i)$ is irrational, so it is $\ne 2$ and hence $\int_0^\infty \frac{x^{k-1}}{e^{2\pi x}-1}dx\ne \sum_{d\ge 1}\frac{d^{k-1}}{e^{2\pi d}-1}$.

Acknowledgement: This solution was initiated by my ambition to understand the "stenographic" answer of @reuns. I would not have come to terms without the challenges and nudges of his answer.

I have tried to make the solution as simple as possible, in particular I have made no reference to special functions, and the only non-elementary ingredient is the partial fraction decomposition of the cotangent.

We have to find the values of $k$ for which

$$S(k) := \sum_{n=1}^{\infty}\frac{n^{k-1}}{e^{2 \pi n}-1}=i(k) :=\int_{0}^{\infty}\frac{n^{k-1}}{e^{2 \pi n}-1}\,dn \tag{1}$$

The integral easily calculated with the result

$$i(k) = \frac{(k-1)!}{(2\pi)^{k-1}}\zeta(k)\tag{2} $$

The tough part is the sum. Instead of the original sum with the power on $n$ we consider the expression

$$g(k,t) = \sum_{n=1}^{\infty}\frac{e^{2\pi n t}}{e^{2 \pi n}-1}\tag{3}$$

from which the sum in question can be generated by differentiating thus

$$S(k) = (-\frac{1}{2\pi}\frac{\partial}{\partial t})^{k-1}g(k,t) |_{t\to 0}\tag{4}$$

Now we carry out some transformations on $g(k,t)$

$$\begin{align} g(k,t) &=\sum_{m=1} ^{\infty}\frac{1}{e^{2\pi(m+t)}-1}\tag a\\ &=\frac{1}{2} \sum_{m=1}^{\infty}\left( \coth\left(\pi(m+t)\right)-1 \right) \tag b\\ &=\frac{1}{2} \sum_{m=1}^{\infty}\left(\frac{i}{\pi} \sum_{n=-\infty}^{\infty}\left( \frac{1}{i(m+t)+n} \right)-1 \right)\tag c\\ &=\frac{1}{2} \sum_{m=1}^{\infty}\left( \frac{i}{\pi} \frac{1}{i(m+t)}\\ + \frac{i}{\pi}\sum_{n \ge 1}\left( \frac{1}{i(m+t)+n}+\frac{1}{i(m+t)-n} \right)-1 \right)\tag d\\ \end{align} $$

Explanation:
$(a)$: $\sum_{n\ge 1}\frac{e^{- 2\pi t n}}{e^{2 \pi n}-1}= \sum_{n\ge 1}e^{-2 \pi n t}\sum_{m\ge 1} e^{-2 \pi n m}=\sum_{m\ge 1}\left(\sum_{n\ge 1}e^{-2 \pi n t} e^{-2 \pi n m}\right)$

$(b)$: $\frac{1}{e^{2a}-1}=\frac{e^{-a}}{e^{a}-e^{-a}}+\frac{1}{2}-\frac{1}{2}=\frac{e^{-a}}{e^{a}-e^{-a}}+\frac{1}{2}\frac{e^{a}-e^{-a}}{e^{a}-e^{-a}} -\frac{1}{2} = \frac{e^{-a}+\frac{1}{2}(e^{a}-e^{-a})}{e^{a}-e^{-a}}-\frac{1}{2}=\frac{1}{2}\left( \frac{e^{a}+e^{-a}}{e^{a}-e^{-a}} -1 \right)$

$(c)$: We use the well-known partial fraction decomposition of the cotangent $\cot(\pi z) = \sum_{n=-\infty}^{\infty} \frac{1}{(n+z)}$ and switch to an imaginary argument: $\cot(i \pi x) = i \coth(\pi x)$

$(d)$: We split the two-sided sum into three parts: $\sum_{n=-\infty}^{\infty} c(n) =\sum_{n=0}^{0} c(n)+\sum_{n=-\infty}^{-1} c(n)+\sum_{n=1}^{\infty} c(n)= c(0) + \sum_{n \ge 1}\left( c(n)+ c(-n) \right)$ so that all indices are in the poistive range.

Hence we find from $(4)$ with $(c)$

$$S(k) = \frac{(k-1)!}{(2\pi)^k}\left(\zeta(k) + W(k) \right)\tag{5}$$

where

$$W(k) = \sum_{n,m \ge 1}\left( \frac{1}{(m-i n)^k} + \frac{1}{(m+i n)^k} \right)\tag{6}$$

Comparing $(2)$ and $(5)$ shows that equality $(1)$ holds if

$$W(k) = 0\tag{7}$$

We can easily find a necessary condition by the following symmetry argument: as the summation is symmetric against the interchange of $n$ and $m$ the sum vanishes if the summand

$$a(n,m) = \frac{1}{(m-i n)^k} + \frac{1}{(m+i n)^k} \tag{8}$$

is antisymmetric.

Starting with the summand $(8)$ in reverse order we have

$$\begin{align} a(m,n) &= \frac{1}{(n-i m)^k} + \frac{1}{(n+i m)^k}\\ &=\frac{(i)^{k}}{(m+i n)^k} + \frac{(-i)^k}{(m-i n)^k}\\ &= i^{k}\left((-1)^k \frac{1}{(m-i n)^k}+\frac{1}{(m+i n)^k} \right)\end{align}\tag{9}$$

In order that this is the negative of $(8)$ we need first $k = 2 p$ and then $i^{2p} = -1$ which requires odd $p$ so that we obtain finally

$$k = 2 (2q+1), \text{or } k \equiv 2 \text{ }(\bmod 4 )\tag{10}$$

Notice the argument breaks down for $k=2$ since the sum $W(k)$ is not konvergent. And, indeed, there is no equality of sum and integral in this case but for $k=2$ we have $i(2) = \frac{1}{24}$ and $S(k) = \frac{1}{24}-\frac{1}{8 \pi}$

The first few coinciding integrals (= sums) are in the format (k,S(k))

$$\left( \begin{array}{cc} 6 & \frac{1}{504} \\ 10 & \frac{1}{264} \\ 14 & \frac{1}{24} \\ 18 & \frac{43867}{28728} \\ 22 & \frac{77683}{552} \\ 26 & \frac{657931}{24} \\ \end{array} \right)$$

Alternative closed form of $S(k)$

An alternative closed form expression of the sum $S(k)$ in terms of known functions was given in $(4)$ of the OP. We provide here the missing derivation:

We consider again the generating function (slightly different from the above developments) which can be transformed as follows (notice that $m=k-1$)

$$\begin{align} & \sum _{n=1}^{\infty } \frac{e^{n t}}{e^{2 \pi n}-1}\\ = & \sum _{j=1}^{\infty } \left(\sum _{n=1}^{\infty } e^{n t} \exp (2 \pi (-j) n)\right)\\ =& \sum _{j=1}^{\infty } \frac{1}{e^{2 \pi j-t}-1}=\sum _{j=1}^{\infty } \frac{e^t}{e^{2 \pi j}-e^t}\\ = &\frac{-\psi _{e^{-2 \pi }}^{(0)}\left(1-\frac{t}{2 \pi }\right)-\log \left(1-e^{-2 \pi }\right)}{2 \pi } \end{align}\tag{11}$$

The last equality was found using Mathematica (which refused to do the sum the step before).

Now the $m$-th derivative at $t=0$ gives $(4)$ of the OP.

Notice that $(11)$ also provides the case $m=0$ ($k=1$).

Things are even simpler! As a matter of fact our sum is not only expressible in terms of the q-polygamma function, it is essentially defining it.

Indeed, the q-digamma function is defined as

$$\psi_{q}(z)=\log (q) \sum _{n=0}^{\infty } \frac{q^{n+z}}{1-q^{n+z}}-\log (1-q)\tag{12}$$

Letting

$$q\to e^{-2 \pi }\tag{13}$$

$(12)$ becomes

$$\psi_{q}(z)=-2 \pi \sum _{n=1}^{\infty } \frac{1}{e^{2 \pi (n+z-1)}-1}-\log \left(1-e^{-2 \pi }\right)\tag{14}$$

Notice that we have shifted the summation index $n\to n-1$ which is adjusted by taking $z \to z-1$.

Comparison with $(11)$ shows agreement.