Is every integer $z$ representable in Pell form as $x^2 \pm dy^2 =z$? We know that there are integers that cannot be represented as the sum of two squares (Fermat's Four Square Theorem).
We also know that every natural number can be represented as the sum of four squares (Lagrange's Four Square Theorem).
Is every integer $z$ representable in Pell form as $x^2 \pm dy^2 = z$, with $d$ being a square-free integer with $|d| > 1$? $d$ is not fixed and cannot be equal to $z$ (since $x=0, y=1, d=z$ would be a trivial solution). Similarly, $x^2$ can be taken to be any square and $y = 1$ and $d = \pm(z - x^2)$ would be a trivial solution.
So, the question is are there any non-trivial solutions $(x, y, d)$ for the equation $x^2 \pm dy^2 = z$?
In other words, I am looking for representing $z$ as the sum (or difference) of a square and $d$ repetitions of a square.
Notes:

*

*$d = 1$ is the Two Square Theorem and $d = -1$ is the factorization of $z$

*$d$ is required to be square-free as the equation would reduce to the Two Squares form otherwise

 A: Let me paraphrase your question as follows:

Determine all $z\in\Bbb{Z}$ for which there exist $d,x,y\in\Bbb{Z}$ with $|d|,|y|>1$ and $d$ squarefree such that
$$x^2+dy^2=z.\tag{1}$$

First note that for $z=0$ there are no integral solutions with $d$ squarefree.
If $z\neq0$ then for every integer $x$ we have the trivial solution
$$(d,x,y)=(z-x^2,x^2,1),$$
which of course fails to meet the condition that $|y|>1$. But for sufficiently large values of $x$ we get
$$d=z-x^2<-1,$$
and so $\Bbb{Z}[\sqrt{-d}]$ is a real quadratic ring. By Dirichlet's unit theorem its unit group has rank $1$, so if $u+v\sqrt{-d}\in\Bbb{Z}[\sqrt{-d}]$ is a fundamental unit and $n\in\Bbb{Z}$ is any integer we have
$$N\left((x+y\sqrt{-d})(u+v\sqrt{-d})^n\right)=N(x+y\sqrt{-d})N(u+v\sqrt{-d})^n=z,$$
yielding infinitely many integral solutions to $(1)$: If $a_n,b_n\in\Bbb{Z}$ are such that
$$a_n+b_n\sqrt{-d}=(x+y\sqrt{-d})(u+v\sqrt{-d})^n,$$
then the above shows that indeed
$$a_n^2+db_n^2=z.$$
Moreover, this yields infinitely many integral solutions $(d,x,y)=(d,a_n,b_n)$ with $|y|>1$, because if $b_m=b_n$ then it quickly follows that $m=n$.
All that remains to be shown is that we can choose $x$ sufficiently large such that $d=z-x^2$ is squarefree.
