Method of proof of $\sum\limits_{n=1}^{\infty}\tfrac{\coth n\pi}{n^7}=\tfrac{19}{56700}\pi^7$

Since $(7)$ from this answer is valid for any $z\in\mathbb{C}$, we have $$ \begin{align} \pi\coth(\pi n) &=\sum_{k\in\mathbb{Z}}\frac1{n+ik}\\ &=\frac1n+2n\sum_{k=1}^\infty\frac1{n^2+k^2}\tag{1} \end{align} $$ Therefore, $$ \begin{align} \sum_{n=1}^\infty\frac{\pi\coth(\pi n)}{n^7} &=\sum_{n=1}^\infty\frac1{n^8}+2\sum_{n=1}^\infty\sum_{k=1}^\infty\frac1{n^6(n^2+k^2)}\\ &=\zeta(8)+2\sum_{n=1}^\infty\sum_{k=1}^\infty\frac1{k^2}\left(\frac1{n^6}-\frac1{n^4(n^2+k^2)}\right)\\ &=\zeta(8)+2\zeta(2)\zeta(6)-2\sum_{n=1}^\infty\sum_{k=1}^\infty\frac1{k^2}\frac1{n^4(n^2+k^2)}\tag{2}\\ &=\zeta(8)+2\zeta(2)\zeta(6)-2\sum_{n=1}^\infty\sum_{k=1}^\infty\frac1{k^4}\left(\frac1{n^4}-\frac1{n^2(n^2+k^2)}\right)\\ &=\zeta(8)+2\zeta(2)\zeta(6)-2\zeta(4)\zeta(4)+2\sum_{n=1}^\infty\sum_{k=1}^\infty\frac1{k^4n^2(n^2+k^2)}\tag{3}\\[6pt] &=\zeta(8)+2\zeta(2)\zeta(6)-\zeta(4)\zeta(4)\tag{4}\\[12pt] &=\frac{19\pi^8}{56700}\tag{5} \end{align} $$ where $(4)$ is the average of $(2)$ and $(3)$. Also. we've used the values of $\zeta(2k)$ computed in this answer. Thus, $$ \sum_{n=1}^\infty\frac{\coth(\pi n)}{n^7}=\frac{19\pi^7}{56700}\tag{6} $$

Suppose we seek to show that $$\sum_{n\ge 1} \frac{\coth(n \pi)}{n^7} = \frac{19\pi^7}{56700}.$$

Using $$\coth(x) = \frac{e^x+e^{-x}}{e^x-e^{-x}} = 1 + 2\frac{e^{-x}}{e^x-e^{-x}}$$ this is the same as $$2\sum_{n\ge 1} \frac{1}{n^7} \frac{e^{-n\pi}}{e^{n\pi}-e^{-n\pi}} = -\zeta(7) + \frac{19\pi^7}{56700}.$$

The sum term may be evaluated using harmonic summation techniques. Since this method has not been presented I will detail this calculation here.

Put $$S(x) = \sum_{n\ge 0} \frac{1}{n^7} \frac{e^{-nx}}{e^{nx}-e^{-nx}}.$$

We will evaluate $S(\pi)$ using a functional equation for $S(x)$ that is obtained 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 = \frac{1}{k^7}, \quad \mu_k = k \quad \text{and} \quad g(x) = \frac{e^{-x}}{e^x-e^{-x}}.$$

We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$\int_0^\infty \frac{e^{-x}}{e^x-e^{-x}} x^{s-1} dx = \int_0^\infty \frac{e^{-2x}}{1-e^{-2x}} x^{s-1} dx \\ = \int_0^\infty \sum_{q\ge 0} e^{-2x} e^{- 2 q x} x^{s-1} dx = \sum_{q\ge 0} \int_0^\infty e^{-2(q+1)x} x^{s-1} dx \\= \Gamma(s) \sum_{q\ge 0} \frac{1}{2^s (q+1)^s} = \frac{1}{2^s} \Gamma(s) \zeta(s)$$ with fundamental strip $\langle 1, \infty\rangle.$

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

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

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

Fortunately the trivial zeros of the two zeta function terms cancel the poles of the gamma function term. Shifting to $\Re(s) = -7 -1/2$ we get $$S(x) = \frac{\pi^8}{18900}\frac{1}{x} - \frac{1}{2} \zeta(7) + \frac{\pi^6 x}{5670} - \frac{\pi^4 x^3}{8100} + \frac{\pi^2 x^5}{5670} + \frac{4}{45} \zeta'(-6) x^6 + \frac{1}{18900} x^7 \\+ \frac{1}{2\pi i} \int_{-15/2-i\infty}^{-15/2+i\infty} Q(s)/x^s ds.$$

We will turn this into the promised functional equation.

Substitute $s = -6 - t$ in the remainder integral to get $$- \frac{1}{2\pi i} \int_{3/2+i\infty}^{3/2-i\infty} \frac{1}{2^{-6-t}} \Gamma(-6-t) \zeta(-6-t) \zeta(1-t) x^{t+6} dt$$ which is $$\frac{x^6}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} 2^{6+t} \Gamma(-6-t) \zeta(-6-t) \zeta(1-t) x^t dt$$

In view of the desired functional equation we now use the functional equation of the Riemann zeta function on $Q(s)$ to prove that the integrand of the last integral is in fact $-Q(t)/\pi^{6+2t}.$

Start with the functional equation $$\zeta(1-s) = \frac{2}{2^s\pi^s} \cos\left(\frac{\pi s}{2}\right) \Gamma(s) \zeta(s)$$ and substitute this into $Q(s)$ to obtain $$Q(s) = 2^{-s} \frac{\zeta(1-s) 2^s \pi^s}{2\cos\left(\frac{\pi s}{2}\right)} \zeta(s+7) = \frac{1}{2} \pi^s \frac{\zeta(s+7)}{\cos\left(\frac{\pi s}{2}\right)} \zeta(1-s).$$ Apply the functional equation again (this time to $\zeta(s+7)$) to get $$Q(s) = \frac{1}{2} \frac{\pi^s}{\cos\left(\frac{\pi s}{2}\right)} \frac{2}{2^{-6-s} \pi^{-6-s}} \cos\left(\frac{\pi (-6-s)}{2}\right) \Gamma(-6-s) \zeta(-6-s) \zeta(1-s)$$ Observe that $$\frac{\cos\left(-3\pi-\frac{\pi s}{2}\right)} {\cos\left(\frac{\pi s}{2}\right)} = - \frac{\cos\left(-\frac{\pi s}{2}\right)} {\cos\left(\frac{\pi s}{2}\right)} = -1$$ so we finally get $$Q(s) = - 2^{6+s} \pi^{6+2s} \Gamma(-6-s) \zeta(-6-s) \zeta(1-s),$$ thus proving the claim.

Return to the remainder integral and re-write it as follows: $$\frac{(x/\pi)^6}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} 2^{6+t} \pi^{6+2t} \Gamma(-6-t) \zeta(-6-t) \zeta(1-t) (x/\pi^2)^t dt.$$ so that the fact of it being a multiple of the defining integral of $S(x)$ becomes readily apparent.

Using the fact that $4/45 \times \zeta'(-6) = -1/2\times\zeta(7)/\pi^6$ we have established the functional equation $$S(x) = \frac{\pi^8}{18900}\frac{1}{x} - \frac{1}{2} \zeta(7) + \frac{\pi^6 x}{5670} - \frac{\pi^4 x^3}{8100} + \frac{\pi^2 x^5}{5670} - \zeta(7) \frac{1}{2\pi^6} x^6 + \frac{x^7}{18900} \\ - \frac{x^6}{\pi^6} S(\pi^2/x).$$

Now the value $x=\pi$ is obviously special here and we get $$S(\pi) = \pi^7 \left(\frac{1}{18900} + \frac{1}{5670} - \frac{1}{8100} + \frac{1}{5670} + \frac{1}{18900}\right) -\zeta(7)- S(\pi)$$ which gives $$2 S(\pi) = \pi^7 \frac{19}{56700} -\zeta(7)$$ as was to be shown.

The inspiration for this calculation is from the paper "Mellin Transform and its Applications" by Szpankowski.

$\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

There is a straightforward way to evaluate the $\texttt{initial (k,n)-sum}$ which is the starting point of @robjohn fine answer:

\begin{align} \sum_{n = 1}^{\infty}{\pi\coth\pars{\pi n} \over n^{7}} & = \sum_{n = 1}^{\infty}{1 \over n^{8}} + 2\sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty}{1 \over n^{6}\pars{n^{2} + k^{2}}} \\[5mm] & = \zeta\pars{8} + \sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty}\bracks{% {1 \over n^{6}\pars{n^{2} + k^{2}}} + {1 \over k^{6}\pars{k^{2} + n^{2}}}} \\[5mm] & = \zeta\pars{8} + \sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty} {k^{6} + n^{6} \over k^{6}n^{6}\pars{n^{2} + k^{2}}} = \zeta\pars{8} + \sum_{n = 1}^{\infty}\sum_{k = 1}^{\infty} {k^{4} - k^{2}n^{2} + n^{4} \over k^{6}n^{6}} \\[5mm] & = \zeta\pars{8} + \sum_{n = 1}^{\infty}{1 \over n^{6}}\sum_{k = 1}^{\infty}{1 \over k^{2}} - \sum_{n = 1}^{\infty}{1 \over n^{4}}\sum_{k = 1}^{\infty}{1 \over k^{4}} + \sum_{n = 1}^{\infty}{1 \over n^{2}}\sum_{k = 1}^{\infty}{1 \over k^{6}} \\[5mm] & = \zeta\pars{8} + 2\zeta\pars{6}\zeta\pars{2} - \zeta^{2}\pars{4} \end{align}