Is there/can there be a model-theoretic proof of this theorem of arithmetic ?

I don't think this can be done. The problem is that if $a$ is not a square, then there are still infinitely many primes $p$ for which $a$ is a square modulo $p$.

Let us make this argument more formal. Suppose we have such a formula $P(x)$ that expresses that $x$ is an integer. As you proposed, we can now look at the theory given by

• the theory of fields,
• the statement $\exists x: (x^2 = 2 \wedge P(x))$ and
• for each prime $p$, the statement $p \neq 0$.

The point is now that this theory does not have a model, since there is no field of characteristic zero in which 2 is the square of an integer. However, every finite subtheory does have a model. Namely, $2$ is the square of an integer in every field $\mathbb{F}_p$ for $p \equiv 1,7 \pmod 8$, and it is well-known that there are infinitely many such primes, so we can just take $p$ big enough. This contradicts the compactness theorem, hence such a formula $P(x)$ does not exists.

Maybe another approach would lead to a model-theoretic proof, but I don't see how one could avoid the above problem.

I don't believe this can work. It is hard to prove that something can't be done, but here is what persuaded me: For every prime $p$, we can solve $x^2+y^2 \equiv -1 \bmod p$, but there is no solution to $x^2+y^2 = -1$ in the integers.

To give an example that uses only one variable, for all primes $p$, we can solve $(z^2-2)(z^2-3)(z^2-6) \equiv 0 \bmod p$, but we can't solve $(z^2-2)(z^2-3)(z^2-6) =0$ in the integers.

So there is no logical principle that let us go from $\forall_p \exists_{x,y} \ (\mbox{polynomial relation}) \bmod p$ to $\exists_{x,y} \ (\mbox{same polynomial relation})$ in the integers. You have to use some special properties of the quadratic polynomial in question, and that strikes me as too specific for model theory.

Proof of the statement that $x^2+y^2 \equiv -1 \bmod p$ is solvable: Clear for $p=2$, so we focus on $p$ odd. The function $x^2+1$ takes $(p+1)/2$ distinct values, and $-y^2$ does as well, so by the pigeonhole principle, we can solve $x^2+1 \equiv -y^2 \bmod p$.

Proof of the statement that $(z^2-2)(z^2-3)(z^2-6) \equiv 0 \bmod p$ is solvable: Clear for $p=2$ or $3$. For $p$ larger than that, note that $$\left( \frac{2}{p} \right) \left( \frac{3}{p} \right) \left( \frac{6}{p} \right) = \left( \frac{36}{p} \right) =1$$ so the residue symbols on the left hand side cannot all be $-1$.