I see: the condition that $d$ be squarefree was added about four hours after I put the majority of this answer. That is, I am allowing $d$ to have prime factors with exponents allowed to be one or larger than one.
Later still: worth pointing out a selection of facts in the other direction. If we have prime $p \equiv 1 \pmod 4,$
then there is an integer solution to $x^2 - p y^2 = -1.$ Proof in Mordell. There are, however, plenty of $d$ that satisfy the conditions but fail: there is no integer solution to $x^2 - 17 y^2 = -1.$ For an example with odd $d,$ there is no integer solution to $x^2 - 205 y^2 = -1.$
Next Day: it is fairly easy, given that $d > 0$ (but not a square) is not divisible by $4$ or by any prime $q \equiv 3 \pmod 4,$ to show that $d$ is the sum of two squares. The hard part, which people seem to be ignoring, is showing that we can arrange $d = u^2 + v^2$ with $\gcd(u,v) = 1,$ meaning $d$ is the sum of coprime squares.
Well, $d$ cannot be divisible by $4$ either.
I had to look for a while to find a correct proof that your $d$ is not just the sum of two squares, it is the sum of two coprime squares. Most books will not show how to do this; I am paraphrasing from Modern Elementary Theory of Numbers, by Leonard Eugene Dickson, especially Theorem 65 on page 63. I am going to discuss many things that will require further study.
We can handle the even case separately, so let $d$ be odd and the product of (possibly several) $p^{e_p,}$ where the prime $p \equiv 1 \pmod 4$ and the exponent $e_p \geq 1.$ An induction argument shows that there is an integer $n > 0$ such that
$$ n^2 \equiv -1 \pmod d. $$
$$ 4 n^2 \equiv -4 \pmod d, $$
$$ 4 n^2 + 4 = d t$$
with integer $t.$ Actually, as $d$ is odd, we see that $t$ is divisible by $4,$ and we may write
$$ 4 n^2 + 4 = 4 d s$$
with integer $s.$ Or,
$$ 4 n^2 - 4 d s = -4. $$ This means that the discriminant of the binary quadratic form $$ \langle d, 2n, s \rangle $$
is $-4.$
I should probably add that $ \langle a,b,c \rangle $ means the quadratic form
$$ f(x,y) = a x^2 + b x y + c y^2, $$
with discriminant
$$ \Delta = b^2 - 4 a c. $$
This means that the form is $SL_2 \mathbb Z$ equivalent ( this is called Gauss reduction) to $ \langle 1,0,1 \rangle. $ Should this be unfamiliar, it means there is an integer matrix $R$ of determinant $1,$
$$
R =
\left(
\begin{array}{cc}
\alpha & \beta \\
\gamma & \delta
\end{array}
\right)
$$
such that
$$ R^T H R = I, $$
with
$$
H =
\left(
\begin{array}{cc}
d & n \\
n & s
\end{array}
\right)
$$
This means that, taking
$$
Q = R^{-1} =
\left(
\begin{array}{cc}
\delta & -\beta \\
-\gamma & \alpha
\end{array}
\right)
$$
we have
$$ Q^T Q = H $$
with $\gcd(\gamma, \delta) = 1.$ And, you see,
$$ \gamma^2 + \delta^2 = d. $$
The even case is just $2d = (\gamma + \delta)^2 + (\gamma - \delta)^2,$ where both numbers are odd.
