Only the start of an answer, but too long for a comment:
$H_k$ is proportional to the $z^{4k-2}$ coefficient of the $\wp$ function associated to an elliptic curve of $j$-invariant $1728$. The Weierstrass differential equation for $\wp$ gives various ways to compute these coefficients, and thus the $H_k$, without recourse to high-precision floating-point arithmetic. For example, writing $\wp$ as the inverse function of an elliptic integral yields the following gp code that recovers all the tabulated numbers; increase $N$ to get more terms:
N = 20;
x = 1 / serreverse(intformal(1/sqrt(1-r^4/16+O(r^(4*N+1)))))^2;
H = vector(N, n, (4*n)! * polcoeff(x,4*n-2) / (4*n-1))
This connection might yield at least some of the results that David observed experimentally. For a start, most of the divisibility by high powers of small primes $p \equiv 3 \bmod 4$ is explained by the factor $k!/(k-1)$ together with observation that in the power-series expansion of $\int (1-(r^4/16))^{-1/2} dr = \sum_m a_m r^m$ the valuation of $a_m$ is not as negative as the usual $-v_p(m)$ — this must be a manifestation of the supersingularity of the curve at such primes; and the primes in the denominator of $H_k$ should be amenable to a similar analysis.