Which quaternary quadratic form represents $n$ the greatest number of times? Let $Q$ be a four-variable positive-definite quadratic form with integer coefficients and let $r_{Q}(n)$ be the number of representations of $n$ by $Q$. The theory of modular forms explains how $r_{Q}(n) = a_{E}(n) + a_{C}(n)$ is broken into the Eisenstein series piece, and the cusp form piece. Knowing what numbers are represented by $Q$ boils down to getting a formula for $a_{E}(n)$ which is of size ''about'' $n$ (as long as local conditions are satified - the formula of Siegel makes this precise), and bounding $|a_{C}(n)|$ (which is $\leq C_{Q} d(n) \sqrt{n}$, where $d(n)$ is the number of divisors of $n$). The hardest part of this is determining $C_{Q}$, and all the methods for doing this rely on bounds for $r_{Q}(n)$.
My question is

If $n$ is a fixed positive integer, what is ${\rm max}_{Q}~r_{Q}(n)$? Here the max runs over all 4-variable positive-definite integral quadratic forms.

My suspicion is that this max always occurs for a form $Q$ with low discriminant, and in particular I would guess the max is always $\leq 24 \sigma(n)$. Equality is achieved for odd $n$ with $Q = x^{2} + y^{2} + z^{2} + w^{2} + xw + yw + zw$, the smallest discriminant quaternary form. (Edit: There are some $n$ that are represented more ways with a form of discrimiant $5$. It seems the right conjecture is $r_{Q}(n) \leq 30 \sigma(n)$.)
I would even be satisfied with a bound of the type $r_{Q}(n) \leq C n \log(n)$ with a value of $C$ that doesn't depend on $Q$. It's possible that one could prove something like this for quadratic polynomials by induction on the number of variables, but I don't see how to make that work. (It seems that $2$ variables is the naturally starting point.)
Another possible approach is something like the circle method, which can recover similar bounds to those which the theory of modular forms gives. In the papers on this subject though, the dependence on the form in the error term seems to make a result of the type I seek difficult.
 A: Too long for comment.
$r_Q(n) \le 24 \sigma(n)$ is false according to my computations.
Take $Q=x^2+y^2+z^2+w^2+x y+x z+x w+y z+y w+z w$.
The sigma inequality fails at least for $n \in \{13, 17, 23, 37, 43\}$
and the ratio $r_Q(n)/\sigma(n)$ is increasing as $n$ increases.
Checked for the unknowns in the range $[-n,n]$.
A: Theorem. Let $Q(x_1,\dots,x_k)$ be a positive definite integral quadratic form in $k\geq 2$ variables. Then the number of integral representations $Q(x_1,\dots,x_k)=n$ satisfies 
$$r_Q(n)\ll_{k,\epsilon}n^{k/2-1+\epsilon}.$$ 
The implied constant depends only on $k$ and $\epsilon$, so it is independent of the actual coefficients of $Q$.
Remark. The "Added 1" section, posted with the permission of Valentin Blomer, contains a more precise result for $k=4$.
Proof. We induct on $k$, and for simplicity we do not indicate the dependence of implied constants on $k$. The case $k=2$ is classical and goes back to Gauss (see the "Added 2" section for more details). So let $k\geq 3$, and assume that the statement holds with $k-1$ in place of $k$. We can assume that
$$ Q(x_1,\dots,x_k)=\sum_{1\leq i,j\leq k} a_{ij}x_ix_j$$
is Minkowski reduced. In particular, $a_{ij}=a_{ji}$ and
$$ 0<a_{11}\leq a_{22}\leq\dots\leq a_{kk}. $$
Then we have a decomposition
$$ Q(x_1,\dots,x_k)=\sum_{i=1}^k h_i\left(\sum_{i\leq j\leq k}c_{ij}x_j\right)^2,$$
where $h_i\asymp a_{ii}$, $c_{ii}=1$ and $c_{ij}\ll 1$ (see Theorem 3.1 and Lemma 1.3 in Chapter 12 of Cassels: Rational quadratic forms). In particular, the coefficients of $Q$ satisfy
$$ a_{11}\dots a_{kk}\asymp h_1\dots h_k=\det(Q),$$
hence also $a_{ij}\ll a_{kk}\ll\det(Q)$ and $h_k\asymp a_{kk}\gg\det(Q)^{1/k}$.
We fix the positive integer $n$, and we consider the integral representations $Q(x_1,\dots,x_k)=n$. The number of representations with $x_k=0$ is 
$\ll_{\epsilon}n^{k/2-3/2+\epsilon}$ by the induction hypothesis, so we can focus on the representations with $x_k\neq 0$. From the above, we see immediately that $x_k\ll\sqrt{n}\det(Q)^{-1/(2k)}$, and then also that $x_{k-1}\ll\sqrt{n}$, then $x_{k-2}\ll\sqrt{n}$, and so on, finally $x_3\ll\sqrt{n}$. It follows that there are $\ll n^{(k-2)/2}\det(Q)^{-1/(2k)}$ choices for the $(k-2)$-tuple $(x_3,\dots,x_k)$ such that $x_k\neq 0$. Fixing such a tuple, 
we are left with an inhomogeneous binary representation problem
$$ a_{11}x_1^2 + 2a_{12}x_1x_2 + a_{22}x_2^2 + d_1 x_1 + d_2 x_2 + e = 0 $$
with fixed integral coefficients $d_1,d_2\ll\sqrt{n}\det(Q)$ and $e\ll n\det(Q)$. Using Lemma 8 in this paper of Blomer and Pohl, it follows that the number of choices for $(x_1,x_2)$ is $\ll_\epsilon n^\epsilon\det(Q)^\epsilon$. Summing up, we get
$$ r_Q(n)\ll_{\epsilon} n^{k/2-3/2+\epsilon} + n^{(k-2)/2+\epsilon}\det(Q)^{-1/(2k)+\epsilon} \ll n^{k/2-1+\epsilon},$$
and we are done.
Added 1. I have been in touch with Valentin Blomer about the original question, and my answer above incorporated a key input from him. More importantly, he realized that the above argument together with some automorphic input allows one to prove, for the case of $k=4$ variables, the striking uniform upper bound (with an absolute implied constant)
$$r_Q(n) \ll \sigma(n).$$ 
Here is the argument of Valentin Blomer, posted with his permission.
For $n\leq\det(Q)^{10}$, the last line of the inductive proof above gives
$$ r_Q(n)\ll_{\epsilon} n^{1/2+\epsilon} + n^{1+\epsilon}\det(Q)^{-1/8+\epsilon} \ll n^{79/80+2\epsilon},$$
so we can (and we shall) assume that $n>\det(Q)^{10}$. We decompose the $\theta$-series of $Q$ uniquely as
$$\theta_Q(z) = E(z) +S(z) = \sum_{n=1}^\infty a(n) e(nz) + \sum_{n=1}^\infty b(n) e(nz)$$ into an Eisenstein series and a cusp form of weight $2$ and level $N$, which is the level of $Q$. Accordingly, $r_Q(n)=a(n)+b(n)$, so it suffices to show that
$a(n)\ll\sigma(n)$ and $b(n)\ll\sigma(n)$. The first bound was proved by Gogishvili (Georgian Math. J. 13 (2006), 687-691.), as follows from (2) and (13)-(14) in his paper. Therefore, it suffices to prove the second bound. We write
$$S = \sum_{f \in B} c(f) f$$
in terms of an orthonormal Hecke eigenbasis $B$ for $S_2(N, \chi)$, where $\chi$ is a quadratic character and the inner product is given by
$$(f, g) = \int_{\Gamma_0(N)\backslash \mathcal{H}} f(z)\bar{g}(z) \frac {dx\, dy}{y^2}.$$
We write $f(z) = \sum_n \lambda_f(n) e(nz)$, so that $b(n) = \sum_f c(f) \lambda_f(n)$. We avoid any use of Eichler/Deligne, among other things because it would require us to deal with oldforms carefully. Instead, we use the Petersson formula and Weil's bound for Kloosterman sums (together with Cauchy-Schwarz and Parseval):
$$\begin{split}
|b(n) |^2 \| S \|_2^{-2} n^{-1} & \leq  n^{-1} \sum_f |\lambda_f(n)|^2  \ll  1 + \sum_{c} \frac{1}{c} S_{\chi}(n, n, c) J_1\left(\frac{4\pi n}{c}\right)\\
& \ll 1 + \sum_{c} \frac{(n, c)^{1/2}\tau(c)}{c^{1/2}} \min\left(\frac{n}{c}, \frac{c^{1/2}}{n^{1/2}}\right) \ll_\epsilon n^{1/2 + \epsilon},
\end{split}$$
so that
$$b(n) \ll_\epsilon \| S \|_2  n^{3/4 + \epsilon}.$$
We have, by Lemma 4.2 of Blomer (Acta Arith. 114 (2004), 1-21.),
$$\| S \|_2 \ll_\epsilon \det(Q)^{2+\epsilon},$$
whence in the end
$$b(n) \ll_\epsilon \det(Q)^{2+\epsilon} n^{3/4 + \epsilon} \leq n^{19/20+2\epsilon}.$$
This concludes the proof. We note that for the twisted Kloosterman sum, the Weil-Estermann bound  is not always true for higher prime powers, see Section 9 of Knightly-Li (Mem. Amer. Math. Soc. 224 (2013), no. 1055), but it is true for the case of quadratic characters that we are using here.
Added 2. Let me provide the details for the case $k=2$. Without loss of generality, 
$$Q(x,y)=ax^2+bxy+cy^2$$
is a reduced form. That is,
$$|b|\leq a\leq c,\qquad\text{whence also}\qquad a\ll\det(Q)^{1/2}.$$
The equation $Q(x,y)=n$ can be rewritten as
$$(2ax+by)^2+4\det(Q)y^2=4an.$$
We can assume that there are (integral) solutions with nonzero $y$, for otherwise there are at most two solutions. In this case,
$$n\geq\det(Q)/a\gg a.$$
The equation factors in the ring of integers of an imaginary quadratic number field, hence a standard divisor bound argument combined with the previous display yields that the number of solutions is 
$$\ll_\epsilon(an)^\epsilon\ll_\epsilon n^{2\epsilon}.$$
A: GH from MO gave in his answer a bound due to himself and Valentin Blomer of $O(\sigma(n))$. I thought it would be interesting to compute the $O$. Here I'm looking for an effective bound of the form $\leq C \Delta^{-\delta} \sigma(n) + o(\sigma(n))$ where the $C$ is explicit but the little $o$ need not be, so we only need be concerned with the Eisenstein term.
The formula in the cited paper is $$\frac{\pi^2 n}{\sqrt{\Delta/16}} \prod_p \chi(p)$$
Define $f_Q(p)$ to be the max over $n$ of $\frac{\chi(p)}{p^{v_p(\Delta)/2}} (\sum_{t=0}^{v_p(n)} p^{-t}) $
Then the main term
$$\frac{4\pi^2 n}{\sqrt{\Delta}} \prod_p \chi(p) \leq  4\pi^2 \sigma(n) \prod_p f_Q(p)$$
So our goal is to upper bound $f_Q(p)$. According to the paper, for primes not dividing the discriminant $\Delta$, $\chi(p)= \left(1- \left(\frac{\Delta}{p} \right)p^{-2}\right) \sum_{t=0}^{v_p(n)}\left(\frac{\Delta}{p} \right)^t p^t$ so $f_Q(p) =  1- \left(\frac{\Delta}{p} \right)p^{-2}$. So the primes not dividing the discriminant contribute $\prod_p\left(1- \left(\frac{\Delta}{p} p^{-2} \right) \right)= L(\chi_{\sqrt{\Delta}},2)^{-1}$. However for a crude upper bound, we instead use $\chi(p) \leq 1+1/p^2$.
For odd primes dividing the discriminant, the bounds given for $\chi(p)$ depend on whether $p$ divides $n$ or not. If $p$ does not divide $n$, they depend further on $n_1$, which is the rank of the quadratic form mod $p$. In this case $\chi(p) \leq 2$ if the rank is $1$ but is at most $1+1/p$ otherwise. The rank can only be $1$ if $v_p(d) \geq 3$, because each term in the formula for the determinant $d$ of the symmetric matrix will be divisible by $p^3$. If $p$ does not divide $n$, the bound for $f_Q(p)$ is $p^{v_p(d)/6} (1+p^{-2}) (1+p^{-1})$. So we have
$$ f_Q(p) \leq \max \left( \frac{1 + 1/p}{p^{v_p(d)/2}}, \frac{1 + 1_{v_p(d) \geq 3}}{p^{v_p(d)/2}}, \frac{ 1 + p^{-2}}{p^{v_p(d)/3}}\right)$$
The third contribution is always the greatest as long as $(1+1/p) \leq p^{1/6} (1+ 1/p^2)$, which happens for all $p>2$, and $2 \leq \sqrt{p} (1+1/p^2)$, which happens for $p>3$ and only fails by a factor of $.962$ for $p=3$.  
The remaining contribution is the local contribution at the prime $2$. The paper gives the bound $4 \cdot 2^{ (v_2(\Delta)-4)/6}$ for $\chi(2)$ and hence $f_Q(2) \leq 2^{4/3} 2^{-v_2(\Delta)/3}$.
Hence $$\prod_p f_Q(p) \leq  \frac{2^{4/3}}{(1+1/4)}\frac{2}{\sqrt{3} (1+ 1/9)}\frac{1}{\Delta^{1/3}} \prod_p (1+1/p^2)$$
Using
$$\prod_p (1+ p^{-2}) = \frac{\zeta(2)}{\zeta(4)}=\frac{15}{\pi^2}$$
we get an upper bound for the main term of
$$ \frac{2^{4/3}}{(1+1/4)}\frac{2}{\sqrt{3} (1+ 1/9)}\frac{60 \sigma(n)}{\Delta^{1/3}}= 125.697\dots \frac{\sigma(n)}{\Delta^{1/3}}$$
To get the ratio under $30$, then, we need $\Delta>73$.
Combined with Valentin and GH's arguments, I believe this implies that there are only finitely many counterexamples to "the number of representatives is at most $30 \sigma(n)$" with $\Delta>73$.
It might be possible to prove a sharp bound by:
1) bounding the error term explicitly to eliminate large $\Delta$ counterexamples.
2) Explicitly calculating the main term for medium $\Delta$, instead of using this crude bound, to eliminate medium $\Delta$ counterexamples.
3) Explicitly calculating the main term and showing the error term vanishes for small $\Delta$, explicitly calculating the highest examples.
