Ramanujan found the following formula:
$$\large \sum_{n=1}^\infty \frac{n^{13}}{e^{2\pi n}-1}=\frac 1{24}$$
I let $e^{2\pi n}-1=\left(e^{\pi n}+1\right)\left(e^{\pi n}-1\right)$ to try partial fraction decomposition and turn the sum into telescoping, but methinks it doesn't lead anywhere and only makes things hairy.
How does one go about proving this? Thanks.
Best Answer
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.