Step 1: Generating function. $r_4(n) = \int_0^1 \theta(\alpha)^4 e^{-2\pi in\alpha}\,d\alpha$ where $\theta(\alpha) = \sum_{m \in \mathbb{Z}} e^{2\pi im^2\alpha}$ is the Jacobi theta function.

Step 2: Major arcs. Near $\alpha = a/q + \beta$ with small $q$, the theta function factors: $\theta(a/q + \beta) = \frac{1}{q}G(a, q) \cdot \theta_0(\beta) + O(q^{1/2+\varepsilon})$ where $G(a,q) = \sum_{r=1}^q e^{2\pi iar^2/q}$ is the Gauss sum (with $|G(a,q)| = \sqrt{q}$ for $(a,q)=1$) and $\theta_0(\beta) = \int_\mathbb{R} e^{2\pi i\beta t^2}dt = (2|\beta|)^{-1/2}e^{\pi i \operatorname{sgn}(\beta)/4}$.

Step 3: Singular series. The major arc contribution gives $r_4(n) = \mathfrak{S}_4(n) \cdot \frac{\pi^2 n}{2} + \text{error}$ where $\mathfrak{S}_4(n) = \sum_{q=1}^\infty A_4(q, n)$ with $A_4(q, n) = q^{-4}\sum_{a=1}^q G(a,q)^4 e^{-2\pi ian/q}$.

Step 4: Computing $\mathfrak{S}_4(n)$. By multiplicativity and evaluation of Gauss sums: $\mathfrak{S}_4(n) = \frac{8}{\pi^2 n}\sum_{\substack{d|n \\ 4\nmid d}} d$ for all $n \geq 1$, so $r_4(n) = 8\sum_{\substack{d|n \\ 4\nmid d}} d$. In particular, $r_4(n) \geq 8 > 0$ for all $n \geq 1$ (taking $d = 1$ in the sum).

Step 5: Minor arcs. For the four-square problem, the minor arc bound is actually unnecessary because Jacobi's formula provides an exact identity via modular forms. The theta function $\theta(\tau)^4$ is a modular form of weight 2 for $\Gamma_0(4)$, and the Fourier coefficients are determined exactly by the Eisenstein series decomposition.

Alternatively, one can verify $r_4(n) > 0$ for all $n$ directly from the formula: $r_4(n) = 8(\sigma(n) - 4\sigma(n/4))$ where $\sigma(m) = \sum_{d|m} d$ and $\sigma$ of a non-integer is 0. Since $\sigma(n) \geq n + 1$ and $4\sigma(n/4) \leq 4\sigma(n/4)$ (only relevant when $4|n$), one verifies $r_4(n) \geq 8$ always. $\square$