On proving that $\sum\limits_{n=1}^\infty \frac{n^{13}}{e^{2\pi n}-1}=\frac 1{24}$

Suppose we seek to evaluate

$$S = \sum_{n\ge 1} \frac{n^{13}}{e^{2\pi n}-1}.$$

This sum may be evaluated using harmonic summation techniques.

Introduce the sum $$S(x; p) = \sum_{n\ge 1} \frac{n^{4p+1}}{e^{nx}-1}$$ with $p$ a positive integer and $x\gt 0.$

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 = k^{4p+1}, \quad \mu_k = k \quad \text{and} \quad g(x) = \frac{1}{e^x-1}.$$

We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$\int_0^\infty \frac{1}{e^{x}-1} x^{s-1} dx = \int_0^\infty \frac{e^{-x}}{1-e^{-x}} x^{s-1} dx \\ = \int_0^\infty \sum_{q\ge 1} e^{-q x} x^{s-1} dx = \sum_{q\ge 1} \int_0^\infty e^{-q x} x^{s-1} dx \\= \Gamma(s) \sum_{q\ge 1} \frac{1}{q^s} = \Gamma(s) \zeta(s).$$

It follows that the Mellin transform $Q(s)$ of the harmonic sum $S(x,p)$ is given by

$$Q(s) = \Gamma(s) \zeta(s) \zeta(s-(4p+1)) \\ \text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1} k^{4p+1} \frac{1}{k^s} = \zeta(s-(4p+1))$$ for $\Re(s) > 4p+2.$

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

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=4p+2) & = \Gamma(4p+2) \zeta(4p+2) / x^{4p+2} \quad\text{and}\\ \mathrm{Res}(Q(s)/x^s; s=0) & = \zeta(0) \zeta(-(4p+1)). \end{align}$$

Computing these residues we get

$$(4p+1)! \frac{B_{4p+2} (2\pi)^{4p+2}}{2(4p+2)! \times x^{4p+2}} = \frac{B_{4p+2} (2\pi)^{4p+2}}{2\times (4p+2) \times x^{4p+2}}$$ and $$- \frac{1}{2} \times -\frac{B_{4p+2}}{4p+2}.$$

This shows that $$S(x;p) = \frac{B_{4p+2} (2\pi)^{4p+2}}{(8p+4)\times x^{4p+2}} + \frac{B_{4p+2}}{8p+4} + \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-(4p+1))$$

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-(4p+1)) \\ = \frac{1}{\sqrt\pi} 2^{s-1} \pi^{s-1/2} \frac{\pi}{\sin(\pi(s+1)/2)} \zeta(1-s)\zeta(s-(4p+1)) \\ = 2^{s-1} \frac{\pi^s}{\sin(\pi(s+1)/2)} \zeta(1-s)\zeta(s-(4p+1)).$$

Now put $s=4p+2-u$ in the remainder integral to get

$$- \frac{1}{x^{4p+2}} \frac{1}{2\pi i} \int_{4p+5/2+i\infty}^{4p+5/2-i\infty} 2^{4p+1-u} \\ \times \frac{\pi^{4p+2-u}}{\sin(\pi(4p+3-u)/2)} \zeta(u-(4p+1))\zeta(1-u) x^u du \\ = \frac{2^{4p+2} \pi^{4p+2}}{x^{4p+2}} \frac{1}{2\pi i} \int_{4p+5/2-i\infty}^{4p+5/2+i\infty} 2^{u-1} \\ \times \frac{\pi^{u}}{\sin(\pi(4p+3-u)/2)} \zeta(u-(4p+1))\zeta(1-u) (x/\pi^2/2^2)^u du.$$

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

We have shown that $$\bbox[5px,border:2px solid #00A000] {S(x;p) = \frac{B_{4p+2} (2\pi)^{4p+2}}{(8p+4)\times x^{4p+2}} + \frac{B_{4p+2}}{8p+4} - \frac{(2\pi)^{4p+2}}{x^{4p+2}} S(4\pi^2/x;p)}.$$

In particular we get

$$S(2\pi; p) = \frac{B_{4p+2}}{8p+4}.$$

The sequence in $p$ starting from $p=1$ is

$${\frac{1}{504}},{\frac{1}{264}},1/24, {\frac{43867}{28728}},{\frac{77683}{552}}, {\frac{657931}{24}},{\frac{1723168255201}{171864}}, \ldots$$

We thus have for $p=3$ as per request from OP

$$\bbox[5px,border:2px solid #00A000]{ \sum_{n\ge 1} \frac{n^{13}}{e^{2\pi n}-1} = \frac{1}{24}.}$$

References, as per request, are: Flajolet and Sedgewick, Mellin transform asymptotics, INRIA RR 2956 and Szpankowski, Mellin Transform and its applications, from Average Case Analysis of Algorithms on Sequences.

It is the weight $14$ Eisenstein series $$G_{14}(z)=\sum_{(n,m)\ne (0,0)} \frac1{(zn+m)^{14}}= 2\zeta(14)+\sum_{n\ne 0} \frac{1}{13!} \frac{d^{13}}{dz^{13}}\frac{2i\pi}{e^{2i\pi n z}-1}$$ $$=2\zeta(14)+\sum_{n\ge 1} \frac{4i\pi}{13!} \sum_{m\ge 1} (2i\pi m)^{13}e^{2i\pi mz}=2\zeta(14)+(2i\pi)^{14}\frac{2}{13!}\sum_{k\ne 1}\frac{k^{13}}{e^{-2i\pi kz}-1} $$

$$G_{14}(z)= z^{-14}G_{14}(-1/z)\implies \qquad G_{14}(i)=0$$

$$\boxed{(2i\pi)^{14}\frac{2}{13!}\sum_{k\ne 1}\frac{k^{13}}{e^{2\pi kz}-1}+2\zeta(14)=0 }$$ $2\zeta(14)=-\frac{B_{14}(2\pi)^{14}}{(14)!} $

For your curiosity !

I do not know if these results are known but, beside this one, $$ \sum_{n=1}^\infty \frac{n^{5}}{e^{2\pi n}-1}=\frac 1{504}=\frac 1{21 \times 24}\qquad\text{and} \qquad \sum_{n=1}^\infty \frac{n^{9}}{e^{2\pi n}-1}=\frac 1{264}=\frac 1{11 \times 24}$$

If they are known, please tell me where I could find them.