Let's consider the double sum over all $\,(n,k)\in\mathbb{Z}^2\,$ except the origin $(0,0)$ :
$$\tag{1}S(s):=\sum_{(n,k)\neq (0,0)}\frac{1}{\left ( n^2+k^2 \right )^s}=4\left(\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^s}+\zeta(2\,s)\right)$$
This is a "Lorenz–Hardy sum" (ref. Lorenz L. $1871$ "Bidrag tiltalienes theori" and Hardy G. H. $1919$ "On some definite integral considered by Mellin" in his Collected papers vol $7$).
Analytical method
The Mellin transform of a function $f$ is defined by :
$$\tag{2}\{\mathcal{M}f\}(s):=\int_0^\infty t^{s-1}f(t)\;dt$$
applied to $\;f:t\mapsto e^{\large{-mt}}\;$ and from the definition of the Gamma function gives :
$$\tag{3}\frac{\Gamma(s)}{\left(m \right)^s}=\int_0^\infty t^{s-1}\,e^{-m\,t}\;dt$$
Supposing that $\Re(s)>1$ we may rewrite our double-sum as :
\begin{align}
\Gamma(s)\,S(s)&=\sum_{(n,k)\neq (0,0)}\frac{\Gamma(s)}{\left ( n^2+k^2 \right )^s}\\
&=\sum_{(n,k)\neq (0,0)}\int_0^\infty t^{s-1}\,e^{-(n^2+k^2)\,t}\;dt\\
&=\int_0^\infty t^{s-1}\sum_{(n,k)\neq (0,0)}\,e^{-(n^2+k^2)\,t}\;dt\\
&=\int_0^\infty t^{s-1}\left(\sum_{n\in\mathbb{Z}}\,e^{-n^2\,t}\sum_{k\in\mathbb{Z}}\,e^{-k^2\,t}-1\right)\;dt\\
\tag{4}&=\int_0^\infty t^{s-1}\left(\theta_3(0,e^{-t})^2-1\right)\;dt\\
\end{align}
i.e. the Mellin transform of $\,f:t\mapsto \theta_3(0,e^{-t})^2-1\,$ where the Jacobi theta function $\theta_3$ is defined by (in this answer we will implicitly suppose that $z=0$ and $q=e^{-t}$) :
$$\tag{5}\theta_3(z,q):=\sum_{n=-\infty}^\infty q^{n^2}e^{2niz}$$
In $1829$ Jacobi published his deep elliptic functions book "Fundamenta nova theoriae functionum ellipticarum" where he obtained numerous identities including the formula $4$) of page $103$ :
$$\tag{6}\theta_3(0,q)^2=\frac {2\,K}{\pi}=1+4\sum_{n=1}^\infty \frac{q^n}{1+q^{2n}}$$
($K$ is the Complete elliptic integral of the first kind but we won't use it here)
The series at the right of $(6)$ is a Lambert series. The expansion of the denominator in power series and substitution of the double-sum in $(4)$ gives :
\begin{align}
\Gamma(s)\,S(s)&=4\int_0^\infty t^{s-1}\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m q^n\,q^{2mn}\;dt\\
&=4\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m \int_0^\infty t^{s-1}e^{-n(2m+1)t}\;dt\\
&=4\,\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m \frac{\Gamma(s)}{(n(2m+1))^s}\\
&=4\,\Gamma(s)\sum_{n=1}^\infty \frac 1{n^s}\sum_{m=0}^\infty \frac{(-1)^m}{(2m+1)^s}\\
&=4\,\Gamma(s)\,\zeta(s)\,\beta(s)\\
\\
&\text{so that }\\
\\
\tag{7}\sum_{(n,k)\neq (0,0)}\frac{1}{\left ( n^2+k^2 \right )^s}&=4\,\zeta(s)\,\beta(s),\quad\Re(s)>1\\
&\text{and}\\
\tag{8}\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^s}&=\zeta(s)\,\beta(s)-\zeta(2s),\quad\Re(s)>1\\
\\\end{align}
$\qquad\qquad\qquad$with $\beta$ the Dirichlet beta function.
This powerful method allows to deduce many closed forms of lattice sums from the known theta identities.
I'll try to add one (or more) alternative derivations when time allows...
References :
- Glasser and Zucker 1980 "Lattice sums" In "Theoretical Chemistry, Advances and Perspectives Vol $5$".
- Borwein and Borwein 1987 "Pi and the AGM".
- Borwein, Glasser, McPhedran, Wan and Zucker 2013 "Lattice Sums Then and Now" (and older papers from these authors)
Best Answer
We have $$\zeta(a)-H_n^{(a)}=\frac{(-1)^{a-1}}{(a-1)!}\int_0^1 \frac{x^n\ln^{a-1}(x)}{1-x}dx$$
Then
$$\sum_{n=1}^\infty\frac{\zeta(a)-H_n^{(a)}}{n}=\frac{(-1)^{a-1}}{(a-1)!}\int_0^1 \frac{\ln^{a-1}(x)}{1-x}\left(\sum_{n=1}^\infty\frac{x^n}{n}\right)dx$$
$$\frac{(-1)^{a-1}}{(a-1)!}\int_0^1 \frac{-\ln^{a-1}(x)\ln(1-x)}{1-x}dx=\frac{(-1)^{a-1}}{(a-1)!}\sum_{n=1}^\infty H_n\int_0^1x^{n}\ln^{a-1}(x)dx$$
$$=\frac{(-1)^{a-1}}{(a-1)!}\sum_{n=1}^\infty H_n\left(\frac{(-1)^{a-1}(a-1)!}{(n+1)^a}\right)=\sum_{n=1}^\infty\frac{H_n}{(n+1)^a}$$
$$=\sum_{n=1}^\infty\frac{H_n}{n^a}-\zeta(a+1)$$
We have the generalized Euler sum
$$\sum_{n=1}^\infty\frac{H_n}{n^a} =\frac{a+2}{2}\zeta(a+1)-\frac12\sum_{j=1}^{a-2}\zeta(a-j)\zeta(j+1)$$
Thus
$$\boxed{\sum_{n=1}^\infty\frac{\zeta(a)-H_n^{(a)}}{n}=\frac{a}{2}\zeta(a+1)-\frac12\sum_{j=1}^{a-2}\zeta(a-j)\zeta(j+1)}.$$
An alternative way using Abel's summation:
$$\sum_{k=1}^n a_k b_k=b_{n}A_n+\sum_{k=1}^{n-1}A_k\left(b_{k}-b_{k+1}\right)$$
where $ A_n=\sum_{i=1}^n a_i$
Let $b_k=\zeta(a)-H_k^{(a)}$ and $a_k=\frac1k$ we have
$$\sum_{k=1}^n \frac{\zeta(a)-H_k^{(a)}}{k}=(\zeta(a)-H_n^{(a)})\sum_{i=1}^n\frac1i+\sum_{k=1}^{n-1}\left(\sum_{i=1}^k\frac1i\right)\left(-H_{k}^{(a)}+H_{k+1}^{(a)}\right)$$
$$=(\zeta(a)-H_n^{(a)})H_n+\sum_{k=1}^{n-1}\left(H_k\right)\left(\frac{1}{(k+1)^a}\right)$$
Let $n\mapsto \infty$ we have
$$\sum_{k=1}^\infty \frac{\zeta(a)-H_k^{(a)}}{k}=0+\sum_{k=1}^\infty\frac{H_k}{(k+1)^a}=\sum_{k=1}^\infty\frac{H_k}{k^a}-\zeta(a+1)$$