I wonder if closed expression can be found for sums of harmonic numbers with a squared argument.
Examples are
$$s_{1}=\sum_{k=1}^\infty \frac{ H(k^2)}{k^2} \simeq 3.28709\tag{1}$$
$$s_{2}=\sum_{k=1}^\infty (-1)^{k+1} \frac{ H(k^2)}{k}\simeq 0.456221\tag{2}$$
$$s_{3}(x)=\sum_{k=1}^\infty x^k H(k^2)\tag{3}$$
I tried the representation
$$H(z) = \int_0^1 \frac{1-x^z}{1-z}\,dx\tag{4}$$
which leads to the sums of the type
$$p_1=\sum _{k=1}^{\infty } \frac{y^{k^2}}{k}\tag{5a}$$ $$p_2=\sum _{k=1}^{\infty } \frac{y^{k^2}}{k^2}\tag{5b}$$
related to Jacobi elliptic functions which I could not evaluate. Can you do better?