What does the Hecke bound $|a_n| < C n$ say about $\theta(z)^4$? I am learning modular forms from notes of Zagier and I am trying to understand the predictions of the Hecke bound.

Let $f(z)$ be a cusp form of weight $k$ on $\Gamma$ with Fourier expansion $\sum a_n q^n$.  Then $|a_n| < Cn^{k/2}$ for all $n$ for some constant $C$ depending only on $f$.

Let's test this out on a modular function.  Consider the theta function: $\theta(z) = \sum q^{n^2}$ taking the 4th power:
$$ \theta(z)^4 = \sum r_4(n) q^n $$
where $r_4(n)$ is the number of representation of $n$ as the sum of 4 squares. 
$\theta(z)^4$ is a modular form of weight 2 on $\Gamma_0(4)$ and the Hecke bound should predict $\boxed{r_4(n) \leq Cn}$ Is this true?
$$ r_4(n) = 8 \sigma(n) - 32 \sigma(n/4) $$
This is taken from Wikipedia and $\sigma(n)$ is the sum of divisors function (OEIS) unfortunatly Wikipedia also has that $\frac{1}{n}r_4(n)$ can be arbitrarily large.
$$ \{ n : \frac{r_4(n)}{n} > N \} \neq \varnothing $$
so it looks like the Hecke bound is wrong in this case... or have I misused it?
 A: In this case
$\theta^4(z)$ is a combination of Eisenstein series:
$$\theta^4(z) =  8 E_2(z) - 32 E_2(4z),$$
where
$$E_2(z) := \dfrac{-1}{24} + \sum_{n=1}^{\infty} \sigma(n) q^n$$
denotes the "weight $2$ Eisenstein series of level one".
You can think of this as an identity of generating functions that
is equivalent to the formula for $r_4(n)$ in terms of divisor sums.
Actually, $E_2(z)$ is not a weight $2$ modular form of level one (which is why I used quotes above); it satisfies a slightly more complicated transformation rule under $z \mapsto -1/z$.  But $E_2^*(z) := E_2(z) - 2 E_2(2z)$ is
a weight $2$ modular form on $\Gamma_0(2)$, and so
$E_2^*(z)$ and $E_2^*(2z)$ are weight $2$ modular forms on $\Gamma_0(4)$.
The above identity may then be rewritten as
$$\theta^4(z) = 8 E_2^*(z) + 16 E_2^*(2z),$$
which is an equation bewteen weight $2$ modular forms on $\Gamma_0(4)$.
Now you can check that there are no cusp forms of weight $2$ on $\Gamma_0(4)$
(e.g. because the modular curve $X_0(4)$ has genus zero), and so the
space of weight $2$ modular forms on $\Gamma_0(4)$ is spanned by the Eisenstein series $E_2^*(z)$ and $E_2^*(2z)$.  Thus $\theta^4(z)$ must be a linear combination of these two Eisenstein series, and to determine the coefficients, it is enough to compare constant and linear terms in the $q$-expansions.
Thus, just from knowing that $r_4(0) = 1$ and $r_4(1) = 8$, we can derive the above formula for $\theta^4(z)$ in terms of Eisenstein series, and thus
prove the general formula for $r_4(n)$ in terms of divisor sums.
(My understanding is that this is more-or-less how Jacobi originally proved
this formula.)
However, if we consider $r_{2k}(n)$ for $k > 2$, then we obtain weight $k$ modular forms (on $\Gamma_0(4)$ if $k$ is even, and on $\Gamma_1(4)$ is $k$ is odd).
If $k \leq 4$, then there are no cuspforms of weight $k$, and so we may again write $\theta^{2k}$ as a linear combination of appropriate Eisenstein series,
leading to the standard formulas for $r_{2k}(n)$ in terms of $\sigma_{k-1}(n)$
when $k \leq 4$.  (The case of $8$ squares, i.e. of $k = 4$, is discussed
here for example.)
Once $k > 4$, you can't expect $\theta^{2k}$ to be a linear combination of Eisenstein series: rather, it will be a linear combination of Eisentein series, plus a cuspform.   This means that $r_{2k}(n)$ will be expressed in terms of $\sigma_{k-1}(n)$, together with the coefficient $a_n$ of a cuspform.
So now Hecke's bound is useful: since $\sigma_{k-1}(n)$ grows at
least as rapidly as $n^{k-1}$, we see from Hecke's bound
that $a_n$ can be thought of
as an error term in the formula for $r_{2k}(n)$, since $\sigma_{k-1}(n)$ grows at least as rapidly as $n^{k-1}$, while the Hecke bound shows that $a_n$
grows at a slower rate than this.
This is a standard way of estimating $r_{2k}(n)$ (and analogous representation number problems for other quadratic forms) when $k$ is large.
