Let (for concreteness) $a = 2$, $b = \sqrt{5}$ and $\varphi = (\sqrt{5}+1)/2$. I am interested in solutions $(w,x,y,z) \in \mathbb{Z}[\varphi]^4$ of the system

$$ w^2 - ax^2 -by^2 + abz^2 = 1 $$ $$ \lvert w^2 + ax^2 +by^2 + abz^2 \rvert \ll \infty $$ $$ \lvert \bar{w}^2 + a\bar{x}^2 -b\bar{y}^2 - ab\bar{z}^2 \rvert \le C $$ for some constant $C$. Here $\overline{\alpha + \beta\sqrt{5}} = \alpha - \beta\sqrt{5}$ and ``$\ll \infty$'' means that ideally I would like to enumerate solutions in increasing order of this value. (Restriction of scalars turns this problem into a system of two quadratic equations and two inequalities in eight variables in $\mathbb{Z}$; if someone wants to see it, I can write it out including potential mistakes).

- What is the best (or even any practical) way to produce these?

I am aware that there is a lot of classical mathematics associated to this question but I don't quite manage to put it together. Perhaps a subquestion is:

- Can one enumerate the squares $s$ in $\mathbb{Z}[\varphi]$ with $\lvert \bar{s} \rvert \le C$ in increasing order of $\lvert s \rvert$?

Context: Let $k = \mathbb{Q}(\varphi)$ and let $A$ be the quaternion algebra $(\frac{a,b}{k})$ with norm $\nu$. With the above values the algebra $A$ is a skew field but tensoring with $\mathbb{R}$ in the two possible ways (taking $\sqrt{5}$ to $\pm\sqrt{5}$) gives an isomorphism with $M_2(\mathbb{R})$ which we equip with the map
$$
\left\lVert\left(\begin{array}{cc}\alpha&\beta\\\gamma&\delta\end{array}\right)\right\rVert = \frac{1}{2}\left(\alpha^2 + \beta^2 + \gamma^2 + \delta^2\right).
$$
The above system then asks for solutions $\lambda$ in ~~the maximal~~ an order of $A$ for $\nu(\lambda) = 1$, $\lVert\lambda\rVert_{\sqrt{5} \mapsto \sqrt{5}} \ll \infty$ and $\lVert \lambda \rVert_{\sqrt{5} \mapsto -\sqrt{5}} \le C$.