(this is not a proof)
Apparently, for squarefree $k$, the limit is $l(k) = \frac 1 2 \prod_{\text {prime }p} f_k(p)$, where :
$$ f_k(p) = \left\{ \begin{array}{ll} \lim_{m \to \infty} E[ \#\{c \in \mathbb Z/p^m\mathbb Z, a^2+b^2 = c^2+k \}]_{(a,b)\in \mathbb (Z/p^m\mathbb Z)^2} & \text {for } p=2 \\
\lim_{m \to \infty} E[1+v_p(c^2+k)]_{c \in \mathbb Z/p^m\mathbb Z} & \text {for } p\equiv 1 \mod 4\\
\lim_{m \to \infty} E[\frac 1 2 (1 + (-1)^{v_p(c^2+k)})]_{c \in \mathbb Z/p^m\mathbb Z} & \text {for } p\equiv 3 \mod 4 \end{array} \right.$$
The argument for this formula is the following : the number of ways to write an integer $n$ as a sum of two squares is the number of elements of norm $n$ in $\mathbb Z[i]$. Modulo the units of $\mathbb Z[i]$, it is the product of the $(1+v_p(n))$ for $n \equiv 1 \mod 4$ and of the $\frac 1 2 (1 + (-1)^{v_p(n)})$ for $n \equiv 3 \mod 4$ (it is a direct consequence of the ideal group structure of $\mathbb Z[i]$). So the factors with $p>2$ express the expected value of the number of elements whose norm is a "randomly chosen" integer. The factor for $p=2$ is there to correct the conspiracy of $a^2+b^2$ not hitting things uniformly in $\mathbb Z_{(2)}$ (for example it never hits things of the form $4q+3$). Finally, the $\frac 1 2$ in front of the product is there because we are counting elements of norm $c^2+k$ not only up to units, but also up to conjugation.
Let $\chi(p) = 1$ if $p \equiv 1 \mod 4$ and $\chi(p) = -1$ if $p \equiv 3 \mod 4$, and I will suppose that $k$ is squarefree.
If $p$ is an odd prime dividing $k$, then $v_p(c^2+k) = \min (v_p(c^2,1)) = 1$ if $p$ divides $c$, $= 0$ otherwise. In this case, $f_k(p) = \frac{p-1}p + \frac 2 p = \frac {p+1} p$ when $\chi(p)=1$ ; and $f(p) = \frac{p-1}p$ when $\chi(p)=-1$ : $f_k(p) = \frac{p+\chi(p)}p$
If $p$ is an odd prime not dividing $k$ and if $-k$ is a square mod $p$, then $-k$ has two distinct square roots in $\mathbb Z_{(p)}$ and we have $v_p(c^2+k) = v_p(c+\sqrt{-k}) + v_p(c- \sqrt{-k})$, so that the "probability" of $v_p(c^2+k) \ge l$ is $2/p^l$. Thus, we get $f_k(p) = 1 + 2\chi(p)/p + 2/p^2 + 2\chi(p)/p^3 + \ldots = 2/(1- \chi(p)\frac 1p)-1 = \frac{p+\chi(p)}{p-\chi(p)}$
If $p$ is an odd prime not dividing $k$ and if $-k$ is not a square mod $p$, then $c^2+k$ can never be a multiple of $p$, so $v_p(c^2+k) = 0$, and $f_k(p) = 1$.
Let $\chi_d(p) = 1$ if $d$ is a square mod $p$, $-1$ if it is not, and $0$ if $p$ divides $d$. We see that we always have $$f_k(p) = \frac{p+\chi(p)}{p-\chi_{-k}(p)\chi(p)} = \frac{p+\chi(p)}{p-\chi_k(p)} = \frac{p^2-1}{(p-\chi_k(p))(p-\chi(p))} $$
Then we recognize an Euler product, which converges (except when $k = 1$, where $\chi_k$ is the trivial character) and which we can compute thanks to the class number formula (be careful because so far we implicitly want to have $2$ in the conductor of $\chi_k$) $$\prod f_k(p) = \frac{f_k(2) L(\chi_k,1)L(\chi,1)}{\frac 34 \zeta(2)} = \frac{f_k(2) h(k) R(k) \phi_2(k) }{(w_k/2) \sqrt{|k|} \phi_\infty(k)}$$
where $h(k)$ is the class number of $\mathbb Q(\sqrt{k})$, $R(k)$ is its regulator,$w_k$ its number of roots of unity, $\phi_2(k) = 3$ if $k \equiv 5 \mod 8$, $1$ otherwise, $\phi_\infty(k) = \frac \pi 2$ if $k>0$, $1$ otherwise.
Now, $f_k(2)$ is not difficult to compute, because the expectation of the number of solutions stabilizes very quickly. The formula agrees with the results presented so far
and also works for positive $k$ :
$\begin{array}{ccccccccccccc} k & -13 & -11 & -10 & -7 & -6 & -5 & -3 & -2 & -1 & +2 & +3 & +5 \\
f_k(2) & 1/2 & 1 & 1/2 & 2 & 1/2 & 1/2 & 1 & 1/2 & 1/2 & 1/2 & 1/2 & 1 \\
l(k) & \frac 1 {2 \sqrt{13}} & \frac 3 {2 \sqrt {11}}
& \frac 1 {2 \sqrt{10}} & \frac 1 {\sqrt 7} & \frac 1 {2 \sqrt{6}}
& \frac 1 {2 \sqrt{5}} & \frac 1 {2 \sqrt{3}} & \frac 1 {4 \sqrt{2}} & \frac 18
& \frac {\log(1+\sqrt 2)}{2 \pi \sqrt 2} & \frac {\log(2+\sqrt 3)}{2 \pi \sqrt 3}
& \frac {3 \log(\frac {1+\sqrt 5}2)}{\pi \sqrt 5}\end{array} $
