Infinite series $n^7/(\exp(2\pi n)-1)$

The above answers do not mention that we can also find closed form expressions by using harmonic sum formulae.

Introduce the sum $$S_m(x) = \sum_{n\ge 1} \frac{(nx)^{2m-1}}{e^{2\pi nx}-1}$$ so that we are looking to find $S_m(1).$

The sum term is harmonic and may be evaluated by inverting its Mellin transform.

Recall the harmonic sum identity $$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) = \left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$ where $g^*(s)$ is the Mellin transform of $g(x).$

In the present case we have $$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad g(x) = \frac{x^{2m-1}}{e^{2\pi x}-1}.$$ We need the Mellin transform of $1/(e^{2\pi x}-1)$ which is $$\int_0^\infty \frac{1}{e^{2\pi x}-1} x^{s-1} dx = \int_0^\infty \frac{e^{-2\pi x}}{1-e^{-2\pi x}} x^{s-1} dx.$$ Expanding the inner sum we obtain $$\int_0^\infty \sum_{q\ge 1} e^{-2\pi qx} x^{s-1} dx = \sum_{q\ge 1} \int_0^\infty e^{-2\pi qx} x^{s-1} dx = \Gamma(s) \sum_{q\ge 1} \frac{1}{(2\pi q)^s} = \frac{1}{(2\pi)^s} \Gamma(s) \zeta(s).$$ It follows that the Mellin transform $g^*(s)$ of $g(x)$ is given by $$\frac{1}{(2\pi)^{s+2m-1}} \Gamma(s+2m-1) \zeta(s+2m-1).$$ Therefore the Mellin transform $Q_m(s)$ of $S_m(x)$ is given by $$\frac{1}{(2\pi)^{s+2m-1}} \Gamma(s+2m-1) \zeta(s)\zeta(s+2m-1).$$

The Mellin inversion integral here is $$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q_m(s)/x^s ds$$ which we evaluate by shifting it to the left for an expansion about zero.

We need to examine which poles of the gamma function term are canceled by the two zeta function terms. The poles of the gamma function term are at all integers less than or equal to $-(2m-1).$ The plain zeta term cancels all the even ones and it also cancels the one at $-(2m-1)+1$ from the compound zeta term when $m>1.$ The compound zeta term cancels all the odd ones starting at $-(2m-1)-2.$ This leaves only two poles at $s=1$ and at $s=-(2m-1).$

Therefore we have the following residues that contribute: $$\mathrm{Res}(Q_m(s)/x^s; s=1) = \frac{1}{(2\pi)^{2m}} \Gamma(2m) \zeta(2m) \frac{1}{x} \\= \frac{1}{(2\pi)^{2m}} (2m-1)! \frac{(-1)^{m+1} B_{2m} (2\pi)^{2m}}{2\times (2m)!} \frac{1}{x} = \frac{(-1)^{m+1} B_{2m}}{2\times 2m} = \frac{(-1)^{m+1} B_{2m}}{4m} \frac{1}{x}$$ and $$\mathrm{Res}(Q_m(s)/x^s; s=-(2m-1)) = \zeta(-(2m-1)) \zeta(0) x^{2m-1} \\= -\frac{1}{2} \times -\frac{B_{2m}}{2m} x^{2m-1}= \frac{B_{2m}}{4m} x^{2m-1}.$$ When $m=1$ there is a pole at $s=0$ which contributes $$ \mathrm{Res}(Q_m(s)/x^s; s=0) = -\frac{1}{4\pi}.$$

We now specialize in the case when $m$ is odd and at least three. Then we observe a fact which is of critical importance, namely that $Q(s)$ is odd on the line $\Re(s) = -(m-1),$ which we now prove.

Recall the formula for $Q(s),$ which was $$\frac{1}{(2\pi)^{s+2m-1}} \Gamma(s+2m-1) \zeta(s)\zeta(s+2m-1).$$ Now put $s = -(m-1)+it$ to get $$\frac{1}{(2\pi)^{m+it}} \Gamma(m+it) \zeta(-m+it+1)\zeta(m+it).$$ Using the functional equation of the Riemann zeta function on the first zeta function term this becomes $$\frac{2}{(2\pi)^{2m}} \Gamma(m+it) \Gamma(m-it) \zeta(m+it) \zeta(m-it) \cos\left(\frac{1}{2}(m-it)\pi\right).$$ Now using the series representation of the zeta function in the half plane $\Re(s)>1$ it is easy to see that the two zeta function terms are conjugates of each other and hence their product is a real number. More importantly the value stays the same when $t$ changes sign. Similarly the limit definition of the gamma function shows that the product of the two gamma function terms is a real number, which also stays the same when $t$ changes sign. The factor in front no longer depends on $t.$ That leaves $$\cos\left(\frac{1}{2}(m-it)\pi\right) = \cos\left(\frac{\pi}{2}m - \frac{\pi}{2}it\right) = (-1)^{\frac{m-1}{2}} \sin\left(\frac{\pi}{2}it\right).$$ This term is indeed odd in $t$ which concludes the argument.

Returning to the original computation we now know to shift the integral from $\Re(s) = 3/2$ to $\Re(s) = -(m-1),$ picking up only the residue at $s=1$ and profiting from the fact that the contribution from the left vertical segment is zero. (Set $x=1$ before the shift.) We have established that for $m$ odd and $m>1,$ $$\sum_{n\ge 1} \frac{n^{2m-1}}{e^{2\pi n}-1} = (-1)^{m+1} \frac{B_{2m}}{4m} = \frac{B_{2m}}{4m}.$$ There is a correction term of $1/\pi/8$ when $m=1$ originating as a half contribution from the pole at zero, which is on the left side of the contour.

Take a look at equation $(6)$ in Ramanujan's Formula for the Logarithmic Derivative of the Gamma Function by David Bradley.

Note that just before the formula it says, "Let $N$ be a positive integer," but the formula is valid for negative $N$ as well.

EDIT: To answer Cody's question in the comments:

When $N$ is a negative even integer, where $N=-2m,$ the sum on the RHS of $(6)$ is taken to be the empty sum, and so is equal to zero. This means that for $N=-4$ we have

$$\sum_{k=1}^\infty \frac{k^7}{ e^{2\pi k} -1 } = \frac{\pi}{8} \sum_{k=1}^\infty \frac{k^8}{ \sinh^2(\pi k) } - \frac{1}{480}$$

and more generally

$$\sum_{k=1}^\infty \frac{k^{4m-1}}{ e^{2\pi k} -1 } = \frac{\pi}{4m} \sum_{k=1}^\infty \frac{k^{4m}}{ \sinh^2(\pi k) } + \frac{B_{4m}}{8m}.$$

The sums you want all appear in $\displaystyle E_{4k}(q)=1-\frac{4k}{B_{2k}}\sum_{n \ge 1}\frac{n^{4k-1} q^n}{1-q^n}\,$, with $q=e^{-2\pi}$ and $B_{2k}$ the ${2k}\,$th Bernoulli number. By the recurrence relation from Wikipedia (which uses $d_k=(4k+6)E_{2k+4}k!\zeta(2k+4)$), all $E_{4k\ge 8}$ can be expressed in terms of $E_4$ and $E_6$ (but $E_6(e^{-2\pi})\,$=$0$, so you only need to consider $E_4$).

Now to find $E_4$, plug in $\tau=i\;$ to the equation for the modular discriminant$$g_2^3(\tau)-27g_3^2(\tau)=(2\pi)^{12}e^{2i\pi\tau}\prod_{k\ge 1}(1-e^{2ik\pi\tau})^{24}$$and use the well known value $\prod_{k\ge 1}(1-e^{-2k\pi})=\frac{e^{\pi/12}\Gamma(1/4)}{2\pi^{3/4}}$ along with $g_2(i)=120\zeta(4)E_4(e^{-2\pi})$ and $g_3(i)=0$ (N.B. $g_3(i)=280\zeta(6)E_6(e^{-2\pi})$)

Finally, the recurrence relations simplify to $E_{4n}=A_n E_4^n$, where $A_n$ is a rather unwieldy rational number (the first few though can be figured easily from the Wikipedia link above, under the subsection "Products of Eisenstein series").

(Edit) I've confirmed the first two sums are $\sum_{k\ge1} \frac{k^3}{e^{2\pi k}-1}=-\frac1{240}+\color{maroon}{\frac{\Gamma^8(1/4)}{5120\pi^6}}$and$\sum_{k\ge1} \frac{k^7}{e^{2\pi k}-1}=-\frac1{480}+\color{maroon}{\frac{3\Gamma^{16}(1/4)}{655360\pi^{12}}}$, but I didn't get consistent results with Maple for the next sum.

(There shouldn't be any weird spaces now.)