The question is in the title. Here is a short motivation. The general quadratic Diophantine equation is
$$
x^TAx+bx+c=0,
$$
where $x$ is a vector of $n$ variables, $A$ is $n \times n$ matrix with integer entries, $b$ vector, $c$ integer. We can try simplify it by linear substitution $x=Hy$, where $y$ are new variables and $H$ is $n\times n$ matrix with integer entries and determinant $\pm 1$. The equation becomes
$$
y^T(H^TAH)y + bHy + c = 0.
$$
The first entry in matrix $H^TAH$ is $h^TAh$, where $h$ is the first column of $H$. If there is a non-zero vector $h$ such that $h^TAh=0$, then in the new equation the $y_1^2$ term vanishes, and equation becomes linear in $y_1$. In many cases, this can be used to solve the equation.
The equation $x_1^2+x_2^2=3x_3^2+1$ is the simplest quadratic equation for which this method does not work. Here $h^TAh=0$ reduces to $h_1^2+h_2^2-3h_3^2=0$, which does not have non-trivial integer solutions. Hence, no linear substitution makes the equation linear in any of the variables. How to solve such equations?
There is a general method outlined by by Grunewald and Segal that should work for all quadratic equations: the idea is that (i) the integral orthogonal group of A is finitely generated, and there is an algorithm for listing the generators, and (ii) there is a finite set of solutions to the equation such that all other solutions can be constructed from this finite set by actions of this integral orthogonal group. However, how to implement this idea for this specific equation? Also, maybe there is an easier method?
I remark that I am looking for more or less explicit description of all integer solutions. The answer "Try $z=0,\pm 1,\pm 2,\dots$ and for each $z$ list all representations of $3z^2+1$ as the sum of squares" does not count.
Best Answer
Ok, I now was able to solve the equation myself.
If $(x,y,z)$ is any solution, then $(y,x,z)$, $(x,-y,z)$, $(x,y,-z)$, and $(x,3z-2y,2z-y)$ are also solutions. To check the last one, observe that $$ x^2 + (3z-2y)^2 - 3(2z-y)^2 = x^2 + y^2 - 3z^2 = 1. $$ All these transformations are invertible: if we apply any of them twice, we go back to the solution we have started with. To check this for the last one, observe that $$ 3(2z-y)-2(3z-2y) = y, \quad 2(2z-y)-(3z-2y)=z. $$ Let us prove that all integer solutions to our equation can be produced from $(x,y,z)=(1,0,0)$ by a sequence of this transformations. This gives a complete description of the solution set.
Let us say that two solutions are in the same orbit if they can be obtained from each other by a sequence of these transformations. We will prove that all solutions are in the same orbit as $(1,0,0)$. In any orbit, choose a solution $(x,y,z)$ with $|z|$ minimal. By swapping $x$ and $y$ and changing signs of $y$ and $z$ if needed, we may assume that $0\leq |x|\leq y$ and $0\leq z$. Because this solution has minimal $|z|$ in the orbit, we must have $|2z-y|\geq z$. If $2z-y\geq 0$, this implies that $2z-y \geq z$, or $z \geq y$. But then $2z^2\geq 2y^2 \geq x^2+y^2 = 3z^2+1$, which is a contradiction. If $2z-y<0$, then $|2z-y|\geq z$ reduces to $y-2z\geq z$, or $y\geq 3z$. But then $9z^2 \leq y^2 \leq x^2+y^2 = 3z^2+1$, which is possible only if $z=0$. Hence, $(x,y,z)=(0,1,0)$, and, by applying the operation $(x,y,z)\to (y,x,z)$ for the last time, we arrive at $(x,y,z)=(1,0,0)$. Hence, all solutions are at the same orbit.