[proofplan]
The proof has three ingredients. First, Euler's four-square identity shows that the set of sums of four squares is closed under multiplication, reducing the problem to primes. Second, for an odd prime $p$, a pigeonhole argument on the $(p+1)/2$ values of $x^2$ and of $-1 - y^2$ modulo $p$ produces $x, y$ with $x^2 + y^2 + 1 \equiv 0 \pmod p$, so some multiple $mp$ with $1 \leq m < p$ is a sum of four squares. Third, a descent argument due to Lagrange refines $m$ step by step to $m = 1$: assuming $mp = x_1^2 + x_2^2 + x_3^2 + x_4^2$ with $m > 1$, one replaces each $x_i$ by its representative $y_i$ in $(-m/2, m/2]$, uses Euler's identity to produce a new representation at level $m' < m$, and iterates.
[/proofplan]
[step:Dispose of the base cases $n = 0, 1, 2$]
The integer $0 = 0^2 + 0^2 + 0^2 + 0^2$ is a sum of four squares (all zero), and $1 = 1^2 + 0^2 + 0^2 + 0^2$ is a sum of four squares. It therefore suffices to prove the theorem for $n \geq 2$.
Moreover, observe that $2 = 1^2 + 1^2 + 0^2 + 0^2$. Combined with Euler's identity below, this reduces the theorem to the case where $n$ is an odd prime.
[/step]
[step:Reduce to primes via Euler's four-square identity]
We invoke Euler's four-square identity, which asserts that for any $a_1, a_2, a_3, a_4, b_1, b_2, b_3, b_4 \in \mathbb{Z}$,
\begin{align*}
(a_1^2 + a_2^2 + a_3^2 + a_4^2)(b_1^2 + b_2^2 + b_3^2 + b_4^2) = c_1^2 + c_2^2 + c_3^2 + c_4^2,
\end{align*}
where
\begin{align*}
c_1 &= a_1 b_1 + a_2 b_2 + a_3 b_3 + a_4 b_4, \\
c_2 &= a_1 b_2 - a_2 b_1 + a_3 b_4 - a_4 b_3, \\
c_3 &= a_1 b_3 - a_2 b_4 - a_3 b_1 + a_4 b_2, \\
c_4 &= a_1 b_4 + a_2 b_3 - a_3 b_2 - a_4 b_1.
\end{align*}
This identity is verified by direct expansion (one may interpret it as the multiplicativity of the norm in the Lipschitz quaternions, where $|q|^2 = a_1^2 + a_2^2 + a_3^2 + a_4^2$ for $q = a_1 + a_2 i + a_3 j + a_4 k$; the identity is then $|q q'|^2 = |q|^2 |q'|^2$).
Consequently, if $n = p_1^{e_1} \cdots p_k^{e_k}$ is the prime factorisation and each prime $p_j$ (including $p_1 = 2$ if present) is a sum of four squares, then by iterating Euler's identity so is any product, hence so is $n$. It therefore suffices to show: every odd prime $p$ is a sum of four squares.
[/step]
[step:Find a small multiple of $p$ that is a sum of four squares via pigeonhole]
[claim:Existence of $x, y$ with $x^2 + y^2 + 1 \equiv 0 \pmod{p}$]
For every odd prime $p$, there exist $x, y \in \mathbb{Z}$ with $0 \leq x, y \leq (p-1)/2$ and $x^2 + y^2 + 1 \equiv 0 \pmod p$.
[/claim]
[proof]
Consider the two sets
\begin{align*}
A &= \{ x^2 \bmod p : 0 \leq x \leq (p-1)/2 \}, \\
B &= \{ -1 - y^2 \bmod p : 0 \leq y \leq (p-1)/2 \}.
\end{align*}
Each set has $(p+1)/2$ elements. Indeed, for $x_1, x_2 \in \{0, 1, \ldots, (p-1)/2\}$ with $x_1^2 \equiv x_2^2 \pmod p$, we have $(x_1 - x_2)(x_1 + x_2) \equiv 0 \pmod p$; since $p$ is prime, $x_1 \equiv \pm x_2 \pmod p$, and the range $0 \leq x_1, x_2 \leq (p-1)/2$ forces $x_1 = x_2$. Similarly for $B$.
Now $|A| + |B| = (p+1)/2 + (p+1)/2 = p + 1 > p = |\mathbb{Z}/p\mathbb{Z}|$, so $A$ and $B$ cannot be disjoint in $\mathbb{Z}/p\mathbb{Z}$. Pick $x, y$ with $x^2 \equiv -1 - y^2 \pmod p$, giving $x^2 + y^2 + 1 \equiv 0 \pmod p$.
[/proof]
Using the claim, choose such $x, y$. Then $x^2 + y^2 + 1^2 + 0^2 = mp$ for some integer $m \geq 1$. Since $x, y \leq (p-1)/2$,
\begin{align*}
mp = x^2 + y^2 + 1 \leq 2 \cdot \left( \frac{p-1}{2} \right)^2 + 1 < \frac{p^2}{2} + 1 < p^2 \quad \text{for } p \geq 2,
\end{align*}
so $m < p$. Hence there exists $m \in \{1, 2, \ldots, p - 1\}$ such that $mp$ is a sum of four squares.
[/step]
[step:Descend from $mp$ to $p$: reduce $m$ until $m = 1$]
Let $m$ be the **least** positive integer for which $mp$ is a sum of four squares, with $mp = x_1^2 + x_2^2 + x_3^2 + x_4^2$. The previous step shows $1 \leq m < p$. We aim to show $m = 1$.
Suppose for contradiction $m \geq 2$.
**Case A: $m$ is even.** Then $x_1^2 + x_2^2 + x_3^2 + x_4^2 \equiv 0 \pmod 2$. The $x_i^2$ residues mod $2$ are $0$ or $1$, and their sum is even, so an even number of the $x_i$ are odd. Relabel (pairing odd with odd, even with even) so that $x_1 \equiv x_2 \pmod 2$ and $x_3 \equiv x_4 \pmod 2$. Then
\begin{align*}
\frac{x_1 + x_2}{2},\ \frac{x_1 - x_2}{2},\ \frac{x_3 + x_4}{2},\ \frac{x_3 - x_4}{2} \in \mathbb{Z},
\end{align*}
and one verifies
\begin{align*}
\left( \frac{x_1 + x_2}{2} \right)^2 + \left( \frac{x_1 - x_2}{2} \right)^2 + \left( \frac{x_3 + x_4}{2} \right)^2 + \left( \frac{x_3 - x_4}{2} \right)^2 = \frac{x_1^2 + x_2^2 + x_3^2 + x_4^2}{2} = \frac{mp}{2}.
\end{align*}
So $(m/2)p$ is a sum of four squares. Since $m/2 < m$ and $m/2 \geq 1$ (using $m \geq 2$), this contradicts minimality of $m$.
**Case B: $m$ is odd (and $m \geq 3$).** For each $i$, define $y_i$ to be the unique integer with
\begin{align*}
y_i \equiv x_i \pmod m, \qquad -\frac{m}{2} < y_i < \frac{m}{2}.
\end{align*}
Then $|y_i| \leq (m-1)/2$ (since $m$ is odd), and
\begin{align*}
y_1^2 + y_2^2 + y_3^2 + y_4^2 \equiv x_1^2 + x_2^2 + x_3^2 + x_4^2 = mp \equiv 0 \pmod m.
\end{align*}
Write $y_1^2 + y_2^2 + y_3^2 + y_4^2 = m m'$ with $m' \geq 0$. We have
\begin{align*}
m m' = \sum_i y_i^2 \leq 4 \left( \frac{m-1}{2} \right)^2 < m^2,
\end{align*}
so $m' < m$. Moreover, $m' > 0$: otherwise all $y_i = 0$, which means $m \mid x_i$ for all $i$; then $m^2 \mid \sum x_i^2 = mp$, i.e. $m \mid p$, contradicting $1 < m < p$ with $p$ prime.
Now apply Euler's four-square identity to the product
\begin{align*}
(mp)(m m') = (x_1^2 + x_2^2 + x_3^2 + x_4^2)(y_1^2 + y_2^2 + y_3^2 + y_4^2) = z_1^2 + z_2^2 + z_3^2 + z_4^2,
\end{align*}
with
\begin{align*}
z_1 &= x_1 y_1 + x_2 y_2 + x_3 y_3 + x_4 y_4, \\
z_2 &= x_1 y_2 - x_2 y_1 + x_3 y_4 - x_4 y_3, \\
z_3 &= x_1 y_3 - x_2 y_4 - x_3 y_1 + x_4 y_2, \\
z_4 &= x_1 y_4 + x_2 y_3 - x_3 y_2 - x_4 y_1.
\end{align*}
Each $z_j$ is divisible by $m$. Indeed, modulo $m$ we have $y_i \equiv x_i$, so
\begin{align*}
z_1 &\equiv x_1^2 + x_2^2 + x_3^2 + x_4^2 = mp \equiv 0 \pmod m, \\
z_2 &\equiv x_1 x_2 - x_2 x_1 + x_3 x_4 - x_4 x_3 = 0 \pmod m,
\end{align*}
and similarly $z_3, z_4 \equiv 0 \pmod m$. Writing $z_j = m w_j$ with $w_j \in \mathbb{Z}$:
\begin{align*}
m^2 m' p = \sum_j z_j^2 = m^2 \sum_j w_j^2,
\end{align*}
so $m' p = w_1^2 + w_2^2 + w_3^2 + w_4^2$. Since $1 \leq m' < m$, this contradicts minimality of $m$.
Both cases yield a contradiction, so $m = 1$. Hence $p = x_1^2 + x_2^2 + x_3^2 + x_4^2$ is a sum of four squares.
[guided]
The descent argument is the heart of the proof, so we explain each move carefully.
We have fixed $p$ and let $m$ be the least positive integer with $mp$ a sum of four squares. We know $1 \leq m < p$ from the pigeonhole step. Our strategy: if $m > 1$, we produce a strictly smaller $m' \geq 1$ with $m' p$ a sum of four squares, contradicting minimality. This forces $m = 1$, which is what we want.
**Why split on $m$'s parity?** The even case admits a direct halving trick; the odd case requires the full Euler-identity machinery.
**Case A (even $m$) in detail.** We have $\sum x_i^2 = mp$ even. A sum of four squares is even iff an even number of summands are odd (since $x^2 \equiv x \pmod 2$ in $\{0, 1\}$). So either $0, 2, 4$ of the $x_i$ are odd. Reorder so that $x_1 \equiv x_2 \pmod 2$ and $x_3 \equiv x_4 \pmod 2$. Then $(x_1 \pm x_2)/2$ and $(x_3 \pm x_4)/2$ are all integers, and the identity
\begin{align*}
\left( \frac{a+b}{2} \right)^2 + \left( \frac{a-b}{2} \right)^2 = \frac{a^2 + b^2}{2}
\end{align*}
applied twice gives $(m/2)p$ as a sum of four squares. Since $m/2 < m$, minimality is violated.
**Case B (odd $m \geq 3$) in detail.** The key is to construct $y_i$ close to zero with $y_i \equiv x_i \pmod m$. Setting $|y_i| \leq (m-1)/2$ gives $\sum y_i^2 \leq m(m-1)^2 / 1 < m^2$, so $mm' = \sum y_i^2 < m^2$ forces $m' < m$. We must rule out $m' = 0$: if $m' = 0$ then all $y_i = 0$, i.e. $m \mid x_i$ for all $i$. But then $m^2 \mid \sum x_i^2 = mp$, so $m \mid p$. Since $p$ is prime and $1 < m < p$, this is impossible. Hence $m' \geq 1$.
Now we combine $mp$ and $mm'$ using Euler's identity to produce $(mm')(mp) = m^2 m' p$ as a sum of four squares. The subtle point is that each component $z_j$ is divisible by $m$: computing modulo $m$, the $y_i$ are indistinguishable from $x_i$, so $z_1$ reduces to $\sum x_i^2 = mp \equiv 0$, and $z_2, z_3, z_4$ reduce to alternating sums of products $x_i x_j$ that cancel pairwise. Dividing out, $w_j = z_j / m$ gives $\sum w_j^2 = m' p$ with $1 \leq m' < m$. Contradiction.
**Why does $m = 2$ not appear as a tricky subcase?** $m = 2$ is even, handled by Case A. The condition $m \geq 3$ for Case B is needed because the representative range $-m/2 < y_i < m/2$ has precisely one integer of each residue class only when $m$ is odd; for even $m$ the boundary case $|y_i| = m/2$ would need special handling, but Case A handles the entire even case more directly.
[/guided]
[/step]
[step:Assemble the conclusion]
We have established: every odd prime $p$ is a sum of four squares. Together with $0 = 0^2 + 0^2 + 0^2 + 0^2$, $1 = 1^2$, and $2 = 1^2 + 1^2$, Euler's four-square identity (which expresses products of sums of four squares as sums of four squares) implies that every positive integer, being a product of primes, is a sum of four squares. This proves the theorem.
[/step]