Evaluating $\sum_{n=1}^\infty \frac{n^3}{e^{2\pi n}-1}$ using inverse Mellin transform

inspiration on the post Evaluating $\sum_{n=1}^{\infty} \frac{n}{e^{2 \pi n}-1}$ using the inverse Mellin transform it is possible to calculate in close form $$\sum _{k=1}^{\infty } -\frac{k^3}{e^{2 \pi k}-1}=\frac{3840 \pi ^4 \psi _{e^{2 \pi }}^{(1)}(1)+480 \pi ^2 \psi _{e^{2 \pi }}^{(3)}(1)-704 \pi ^6-5760 \pi ^5+3 \Gamma \left(\frac{1}{4}\right)^8}{23040 \pi ^6}$$ using euler sum I appreciatte some comment. I like to give another series it will interesting using elliptic theta function theory to prove it $$\sum _{k=1}^{\infty } \frac{\left(k \left(-\log \left(\frac{\pi }{2}\right)\right)\right)^3}{e^{2 \pi \left(k \log \left(\frac{\pi }{2}\right)\right)}+1}$$ $$\frac{\log ^4\left(\frac{\pi }{2}\right) \psi _{e^{-\frac{\pi }{\log \left(\frac{\pi }{2}\right)}}}^{(3)}(1)-\log ^4\left(\frac{\pi }{2}\right) \psi _{e^{-\frac{\pi }{\log \left(\frac{\pi }{2}\right)}}}^{(3)}\left(-\frac{\left(i \pi -\frac{\pi }{\log \left(\frac{\pi }{2}\right)}\right) \log \left(\frac{\pi }{2}\right)}{\pi }\right)}{16 \pi ^4 \log \left(\frac{\pi }{2}\right)}-\frac{1}{240} \log ^3\left(\frac{\pi }{2}\right)-\frac{7}{1920 \log \left(\frac{\pi }{2}\right)}$$ sorry for the latex type i do not to improve

This is an alternative approach which is too long for comment. If we put $q = e^{-\pi}$ then the desired sum is $$-\sum_{n = 1}^{\infty}\frac{n^{3}}{q^{-2n} - 1} = -\sum_{n = 1}^{\infty}\frac{n^{3}q^{2n}}{1 - q^{2n}} = \frac{1 - Q(q^{2})}{240}$$ and we know that $$Q(q^{2}) = \left(\frac{2K}{\pi}\right)^{4}(1 - k^{2} + k^{4})$$ For $q = e^{-\pi}$ we have $k = 1/\sqrt{2}$ so that $(1 - k^{2} + k^{4}) = 1 - 1/2 + 1/4 = 3/4$ and $$K(k) = \frac{1}{4\sqrt{\pi}}\Gamma\left(\frac{1}{4}\right)^{2}$$ so that $$Q(q^{2}) = \frac{3}{4}\left(\frac{2K}{\pi}\right)^{4} = \frac{3}{64\pi^{6}}\Gamma\left(\frac{1}{4}\right)^{8}$$ and therefore the desired sum is equal to $$\frac{64\pi^{6} - 3\Gamma(1/4)^{8}}{15360\pi^{6}}$$ The expression for $Q(q^{2})$ in terms of $K, k$ is derived here. Using the same technique and expression for $R(q^{2})$ we can get the surprisingly simple and beautiful result $$\sum_{n=1}^{\infty}\frac{n^{5}}{e^{2\pi n} - 1} = \frac{1}{504}$$ For the second sum mentioned in the question we let $$q = \exp(-2\pi\log(\pi/2)) = \left(\frac{2}{\pi}\right)^{2\pi}$$ and then the desired sum is equal to $$S = -\left(\log\left(\frac{\pi}{2}\right)\right)^{3}\sum_{n = 1}^{\infty}\frac{n^{3}q^{n}}{1 + q^{n}} = -\left(\log\left(\frac{\pi}{2}\right)\right)^{3}\cdot A$$ where the sum \begin{align} A &= \sum_{n = 1}^{\infty}\frac{n^{3}q^{n}}{1 + q^{n}}\notag\\ &= \sum_{n = 1}^{\infty}n^{3}q^{n}\left(\frac{1}{1 + q^{n}} - \frac{1}{1 - q^{n}} + \frac{1}{1 - q^{n}}\right)\notag\\ &= \sum_{n = 1}^{\infty}n^{3}q^{n}\left(\frac{1}{1 - q^{n}} - \frac{2q^{n}}{1 - q^{2n}}\right)\notag\\ &= \sum_{n = 1}^{\infty}\frac{n^{3}q^{n}}{1 - q^{n}} - 2\sum_{n = 1}^{\infty}\frac{n^{3}q^{2n}}{1 - q^{2n}}\notag\\ &= \frac{Q(q) - 1}{240} - \frac{Q(q^{2}) - 1}{120}\notag\\ &= \frac{1 + Q(q) - 2Q(q^{2})}{240}\notag\\ &= \frac{1}{240} + \frac{1}{240}\left(\frac{2K}{\pi}\right)^{4}(1 + 14k^{2} + k^{4} - 2 + 2k^{2} - 2k^{4})\notag\\ &= \frac{1}{240} - \frac{1}{240}\left(\frac{2K}{\pi}\right)^{4}(1 - 16k^{2} + k^{4})\notag\\ \end{align} and hence $$S = \frac{1}{240}\left(\log\left(\frac{\pi}{2}\right)\right)^{3}\left(\frac{2K}{\pi}\right)^{4}(1 - 16k^{2} + k^{4}) - \frac{1}{240}\left(\log\left(\frac{\pi}{2}\right)\right)^{3}$$ I doubt if it can be put into a closed form which is as simple as that for the previous sum.

Suppose we seek a functional equation for the sum term
$$S(x) = \sum_{k\ge 1} \frac{k^5}{\exp(kx)-1}$$

which is harmonic and may be evaluated by inverting its Mellin transform. We are interested in possible fixed points of the functional equation especially $S(2\pi)$.

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) = \mathfrak{M}\left(g(x);s\right)$ is the Mellin transform of $g(x).$

In the present case we have $$\lambda_k = k^5, \quad \mu_k = k \quad \text{and} \quad g(x) = \frac{1}{\exp(x)-1}.$$

We need the Mellin transform $g^*(s)$ of $g(x)$ which is computed as follows:

$$g^*(s) = \int_0^\infty \frac{1}{\exp(x)-1} x^{s-1} dx = \int_0^\infty \frac{\exp(-x)}{1-\exp(-x)} x^{s-1} dx \\ = \int_0^\infty \sum_{q\ge 1} \exp(-qx) x^{s-1} dx = \sum_{q\ge 1} \frac{1}{q^s} \Gamma(s) = \Gamma(s) \zeta(s).$$

Hence the Mellin transform $Q(s)$ of $S(x)$ is given by $$ Q(s) = \Gamma(s) \zeta(s) \zeta(s-5) \quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1} \frac{k^5}{k^s} = \zeta(s-5)$$ where $\Re(s) > 6$.

Intersecting the fundamental strip and the half-plane from the zeta function term we find that the Mellin inversion integral for an expansion about zero is $$\frac{1}{2\pi i} \int_{13/2-i\infty}^{13/2+i\infty} Q(s)/x^s ds$$ which we evaluate in the left half-plane $\Re(s)<13/2.$

The two zeta function terms cancel the poles of the gamma function term and we are left with just

$$\begin{align} \mathrm{Res}(Q(s)/x^s; s=6) & = \frac{8\pi^6}{63x^6} \quad\text{and}\\ \mathrm{Res}(Q(s)/x^s; s=0) & = \frac{1}{504}. \end{align}$$

This shows that $$S(x) = \frac{8\pi^6}{15x^6} + \frac{1}{504} + \frac{1}{2\pi i} \int_{-1/2-i\infty}^{-1/2+i\infty} Q(s)/x^s ds.$$

To treat the integral recall the duplication formula of the gamma function: $$\Gamma(s) = \frac{1}{\sqrt\pi} 2^{s-1} \Gamma\left(\frac{s}{2}\right) \Gamma\left(\frac{s+1}{2}\right).$$

which yields for $Q(s)$

$$\frac{1}{\sqrt\pi} 2^{s-1} \Gamma\left(\frac{s}{2}\right) \Gamma\left(\frac{s+1}{2}\right) \zeta(s) \zeta(s-5)$$

Furthermore observe the following variant of the functional equation of the Riemann zeta function: $$\Gamma\left(\frac{s}{2}\right)\zeta(s) = \pi^{s-1/2} \Gamma\left(\frac{1-s}{2}\right) \zeta(1-s)$$

which gives for $Q(s)$ $$\frac{1}{\sqrt\pi} 2^{s-1} \pi^{s-1/2} \Gamma\left(\frac{s+1}{2}\right) \Gamma\left(\frac{1-s}{2}\right) \zeta(1-s)\zeta(s-5) \\ = \frac{1}{\sqrt\pi} 2^{s-1} \pi^{s-1/2} \frac{\pi}{\sin(\pi(s+1)/2)} \zeta(1-s)\zeta(s-5) \\ = 2^{s-1} \frac{\pi^s}{\sin(\pi(s+1)/2)} \zeta(1-s)\zeta(s-5).$$

Now put $s=6-u$ in the remainder integral to get

$$\frac{1}{x^6} \frac{1}{2\pi i} \int_{13/2-i\infty}^{13/2+i\infty} 2^{5-u} \frac{\pi^{6-u}}{\sin(\pi(7-u)/2)} \zeta(u-5)\zeta(1-u) x^u du \\ = \frac{64 \pi^6}{x^6} \frac{1}{2\pi i} \int_{9/2-i\infty}^{9/2+i\infty} 2^{u-1} \frac{\pi^{u}}{\sin(\pi(7-u)/2)} \zeta(u-5)\zeta(1-u) (x/\pi^2/2^2)^u du.$$

Now $$\sin(\pi(7-u)/2) = \sin(\pi(-u-1)/2+4\pi) \\ = \sin(\pi(-u-1)/2) = - \sin(\pi(u+1)/2).$$

We have shown that $$S(x) = \frac{8\pi^6}{63x^6} +\frac{1}{504} - \frac{64\pi^6}{x^6} S(4\pi^2/x).$$

In particular we get $$S(2\pi) = \frac{1}{63\times 8} + \frac{1}{504} - S(2\pi)$$

or $$S(2\pi) = \frac{1}{504}.$$

Remark. Unfortunately this method does not work for $$S(x) = \sum_{k\ge 1} \frac{k^3}{\exp(kx)-1}$$

We get the functional equation $$S(x) = \frac{\pi^4}{15x^4} -\frac{1}{240} + \frac{16 \pi^4}{x^4} S(4\pi^2/x).$$

which yields $$S(2\pi) = \frac{1}{15\times 16} -\frac{1}{240} + S(2\pi)$$

which holds without providing any data about the value itself.