Sum of the inverse squares of the hypotenuse of Pythagorean triangles

What is the sum of the series

$$ S = \frac{1}{5^2} + \frac{1}{13^2} + \frac{1}{17^2} + \frac{1}{25^2} + \frac{1}{29^2} + \frac{1}{37^2} + \cdots $$

where the sum is taken over all hypotenuse of primitive Pythagorean triangles.

By numerical computation, I found the sum to be $0.056840308812554488$ correct to $18$ decimal places. I would like to know if this sum has a closed form.

Using the general formula for primitive Pythagorean triangles, $$ S = \sum_{r>s\ge 1, \\ \gcd(r,s)= 1}\frac{1}{(r^2 + s^2)^2} $$

Trivially, for all primitive and non primitive Pythagorean triangles, the sum will be $\zeta(2) = \pi^2/6$ times the corresponding sum for primitive Pythagorean triangles which turn out to be about $0.09349856033594433852$.

Motivation: We have equated the sum of the square of the sides of a right triangle to the square of the hypotenuse, so I was curious to know what the sum of the reciprocal of the square of the hypotenuse would be. Also since $\zeta(2)$ converges, and the density of hypotenuse is smaller than that of natural numbers, this sum must trivially converge.

Related question: What is the sum of the reciprocal of the hypotenuse of Pythagorean triangles?

Let us use the notations from

http://mathworld.wolfram.com/DoubleSeries.html

We fix a positive definite binary quadratic form $q$ given by $q(m,n)=am^2+bmn+cn^2$, $a,b,c$ integers. We use sumations over the index set $$J=\Bbb Z\times\Bbb Z-\{(0,0)\}\ .$$ We define $$ \begin{aligned} S(q;s) = S(a,b,c;s) &=\sum_{(m,n)\in J} q(m,n)^{-s}=\sum_{(m,n)\in J} (am^2+bmn+cn^2)^{-s}\ ,\\ S_1(q;s) = S_1(a,b,c;s) &=\sum_{(m,n)\in J} \color{blue}{(-1)^m}\; q(m,n)^{-s}\ ,\\ S_2(q;s) = S_2(a,b,c;s) &=\sum_{(m,n)\in J} \color{blue}{(-1)^n}\; q(m,n)^{-s}\ ,\\ S_{12}(q;s) = S_{12}(a,b,c;s) &=\sum_{(m,n)\in J}\color{blue}{(-1)^{m+n}}\; q(m,n)^{-s}\ . \end{aligned} $$ The last three sums are "twisted versions" of the first sum, the "twist" occurs by using a character for the first parameter, for the second one, for both. In our case, $q(m,n)=m^2 +n^2$, and $(a,b,c)=(1,0,1)$, we have a symmetric case (w.r.t. the exchange $a\leftrightarrow c$).

We will drop $q$ below from notations in $S_?(q,s)$, since we use only the above quadratic form $q$. I decided during the edit operation that should bring us quick to numbers we can compute that it is better for the check to introduce the versions $S^+$ for all sums, where the plus index indicates a further restriction to $(m,n)\in J$ with $$(+)\qquad m,n>0\ .$$

From loc. cit. we extract the following relations: $$ \begin{aligned} S(s) &= \sum_{(m,n)\in J}(m^2+n^2)^{-s} \\ &= 4\beta(s)\;\zeta(s)\ ,\\ %S_1(s) =S_2(s) &= \sum_{(m,n)\in J}(-1)^m(m^2+n^2)^{-s} =\dots %= 2^{-s}b_2(2s) = -2^{-s}\cdot 4\beta(2s)\; \eta(2s) %\ , %\\ S_{12}(s) &= \sum_{(m,n)\in J}(-1)^{m+n}(m^2+n^2)^{-s} \\ &= -4 \beta(s) \;\eta(s)=-4\beta(s) \;(1-2^{1-s})\; \zeta(s)\ . \\[2mm] &\qquad\text{ Then the plus versions are:} \\[3mm] S^+(s) &= \beta(s)\;\zeta(s) - \zeta(2s)\ , \\ -S_{12}^+(s) &= \beta(s)\;\eta(s) - \eta(2s) \\ &= \beta(s)\;(1-2^{1-s})\zeta(s) - (1-2^{1-2s})\zeta(2s)\ , \\ &\qquad \text{ which gives} \\ S^+(s)-S_{12}^+(s) &=2\beta(s)\;(1-2^{-s})\zeta(s) - 2(1-2^{-2s})\zeta(2s)\ . \end{aligned} $$

Now let us search for a linear combination of the above sums that correspond to summing $q(m,n)^{-s}$ over the set of $K$ of all $(m,n)$ with positive (components with) different parity. This is $$\frac 12(\ S^+(s)-S^+_{12}(s)\ )\ .$$ So far can we write: $$ \begin{aligned} &\beta(s)\;(1-2^{-s})\zeta(s) - (1-2^{-2s})\zeta(2s) \\ &\qquad=\frac 12(\ S^+(s)-S^+_{12}(s)\ ) \\ &\qquad=\sum_{\substack{(m,n)\in K\\m,n> 0}} q(m,n)^{-s}\\ &\qquad=2\sum_{\substack{(m,n)\in K\\m>n> 0}} q(m,n)^{-s}\\ &\qquad=2\sum_{\substack{(m,n)\in K\\m>n> 0\\ d=(m,n)\text{ odd}}} q(m,n)^{-s}\qquad\text{ and with }M=m/d,\ N=n/d\\ &\qquad=2\sum_{d>0\text{ odd}}d^{-2s} \sum_{\substack{(M,N)\in K\\M>N> 0\\ (M,N)=1}} q(M,N)^{-s}\\ &\qquad= 2(1-2^{-2s})\; \zeta(2s)\; \color{blue}{ \sum_{\substack{(M,N)\in K\\M>N> 0\\ (M,N)=1}} q(M,N)^{-s} } \ . \end{aligned} $$ The isolated sum in the last expresson is the sum we need, let us take it for $s=2$.

The value we obtain is: $$ \color{brown}{ \frac{\beta(2)\;\zeta(2)}{2(1+2^{-2})\zeta(4)} -\frac 12\ = \frac{6C}{\pi^2} - \frac 12.} $$

$$ \color{brown}{ \frac{\beta(2)\zeta(2)}{2(1+2^{-2})\zeta(4)} - \frac{1}{2} = \frac{6C}{\pi^2} - \frac 1 2.} $$

where $C$ is the Catalan constant. Numerically:

sage: E = catalan * zeta(2) / 2 / (1+2^-2) / zeta(4) - 1/2 
sage: E.n()
0.0568403090661582