I'm interested in a handy way, if it exists, to generate all rational solutions of $x^2 + y^2 = z (z^2 - 1)$.


Clearly there are rational solutions for $z = a/b$ (where $a, b$ are coprime integers with $0 < b < a$) iff $a, b, a - b, a + b$ are each expressible as a sum of two squares (of integers of course), or equivalently that none has a prime divisor of the form $4 Z + 3$ to odd multiplicity.


A perl search script for $0 < b < a < 1000$ showed the proportion of "eligible" coprime pairs $a, b$ to be about 0.01, reducing as the upper bound increases. (It would be interesting to know where this proportion is heading asymptotically: 1/128 perhaps, or something obscure involving $\pi$ ? )


For any eligible z, which gives one rational pair $x, y$, all such pairs can obviously be found by composing the LHS with $p^2 + q^2 = 1$.


I wondered if the latter trick could be exploited to ensure some $x$ and/or $y$ of a special form from which all eligible values of $z$ could be algebraically expressed, in effect a parametric solution. But my efforts to obtain a general two-parameter solution have failed, and I'd be (pleasantly) surprised if there is one, although there are of course special one-parameter solutions.


So, in summary, are any references to this equation? (Google is little help when searching for formulae!) Or perhaps the recently developed methods to apply Mordell-Weil groups to algebraic surfaces may be applicable to it?



Rational solutions of $x^2 + y^2 = z (z^2 - 1)$
I doubt there is a simple complete parametrization. Here are a number of comments:


**There is a two-dimensional family of solutions** Start with two solutions that do not have $x^2+y^2=0$. For concreteness sake, I'll choose $$\left( \frac{3}{8} \right)^2+ \left( \frac{6}{8} \right)^2 = \left( \frac{5}{4} \right)^3 - \frac{5}{4}$$
and
$$\left( \frac{16}{27} \right)^2 + \left(\frac{2}{27} \right)^2=\left( \frac{-4}{9} \right)^3 - \left( \frac{-4}{9} \right).$$


Each such solution generates a curve of solutions: Set 
$$p(t) = \left( \frac{3}{8} \frac{1-t^2}{1+t^2} + \frac{6}{8} \frac{2t}{1+t^2}, - \frac{3}{8} \frac{2t}{1+t^2} + \frac{6}{8} \frac{1-t^2}{1+t^2}, \frac{5}{4} \right)$$
and
$$q(u) = \left( \frac{16}{27} \frac{1-u^2}{1+u^2} + \frac{2}{27} \frac{2u}{1+u^2}, - \frac{16}{27} \frac{2u}{1+u^2} + \frac{2}{27} \frac{1-u^2}{1+u^2}, \frac{-4}{9} \right).$$
Geometric, this is just rotating the $(x,y)$ coordinates by $2 \tan^{-1}(t)$.


For any $t$ and $u$, the line through $p(t)$ and $q(u)$ meets the cubic at three points. Two of them are $p(t)$ and $q(u)$, call the third one $\phi(t,u)$. This is given by some complicated rational function of $t$ and $u$, which I leave it to you to compute.


This is called a "degree $4$ unirational parametrization", because recovering $t$ and $u$ from $\phi(t,u)$ requires solving a degree $4$ equation. There are also degree $2$ unirational parametrizations. (Sketch: Take a smooth rational point on $p$ the cubic. Let $T\_p$ be the tangent plane to $p$ and let $C\_p$ be the intersection of $T\_p$ with the cubic. Then $C\_p$ is a nodal plane cubic, and hence has a rational parametrization $\alpha(t)$. Also, the projective completion $(x^2+y^2)w = z(z^2-w^2)$ contains a line at $(w:x:y:z) = (0:\ast:\ast: 0)$, let $\beta(u)$ parametrize this line. The line through $\alpha(t)$ and $\beta(u)$ intersects the cubic at one additional point.)


**There is no rational parametrization of all solutions** Specifically, there does not exist rational functions $(f(t,u), g(t,u), h(t,u))$ so that $f^2+g^2 = h^3-h$ and, letting $t$ and $u$ range over $\mathbb{Q}^2$, we sweep out all solutions. To see this, note that all solutions have either $-1 \leq z \leq 0$, or $1 \leq z$. (I deliberately chose my example solutions to contain one of each type.) A rational parametrization couldn't jump over this gap, so it will only hit solutions of one of these two types. (This argument is taken from Section 1.4 of [Karen Smith's notes](http://www.math.lsa.umich.edu/~kesmith) on rational varieties.)


Moreover, I highly doubt this is the only restriction; I would guess that each unirational parametrization will hit only a small fraction of the total list of points. I don't have an argument or reference for this, though.


**The density of solutions goes to zero** As you say, the goal is to find $a$ and $b$ coprime such that $a$, $b$, $a+b$ and $a-b$ are all sums of squares. For any prime $p \equiv 3 \mod 4$, the probability that a random $(a,b)$ will not be divisible by $p$ exactly only once is something like $1-4/p+1/p^2$. So, intuitively, the probability that this condition is obeyed for all $p \equiv 3 \mod 4$, should behave like $\prod\_{p \equiv 3 \bmod 4} (1-4/p+1/p^2)$, which diverges to $0$. It should be easy to make this precise.


The big paper on counting points on hypersurfaces seems to be "[The density of rational points on curves and surfaces](http://www.jstor.org/stable/3062125)" by Heath-Brown. His framework is to count integer solutions to 
$$(x^2+y^2)w = z^3 - z w^2$$
with $|w|$, $|x|$, $|y|$, $|z| \leq B$ and $GCD(w,x,y,z)=1$, discarding any that lie on lines contained in the surface. In this case, I think that just means he is discarding $(w,x,y,z) = (0,\ast, \ast,0)$, which you don't care about anyway. His best bound is $O(B^{52/27+\epsilon})$ (Theorem 6); he conjectures that the true bound is $O(B^{1+\epsilon})$.


There are 57 papers in Mathscinet citing this one, so it is obviously an active field; I didn't check to see whether the bounds have been improved. I also didn't think about the question of just counting the number of distinct values for $z/w$.


