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}.$$
Best Answer
This is a problem I have thought alot about. I have not seen any of the modern techniques in your list applied to the problem. Part of the issue is that if you represent $\sigma(n)=2n$ as a Diophantine equations in $k$ variables (corresponding to the prime factors--but allowing the powers to vary) then there are lots of solutions (just not where all the variables are simultaneously prime). So the usual methods of trying to show non-existence of solutions just don't cut it. Historically, this multiplicative approach is the one many people have taken, because at least some progress can be made on the problem. My personal feeling is that maybe someday these bounding computations will be tweaked to the point that they lead to the discovery of some principle that will solve the problem. For example, in one of my recent papers, I was led to consider the gcd of $a^m-1$ and $b^n-1$ (where $a$ and $b$ are distinct primes). I would conjecture that this gcd has small prime factors unless $m$ or $n$ is huge. If that happens, many of the computations related to bounding OPNs become much easier.
I have occasionally thought about whether modular forms might say something about this topic (which is why I'm currently sitting in on my colleague's course). Instead of $\sigma(n)$, the `right' function to consider is $\sigma_{-1}(n)=\sigma(n)/n$ and I don't know off the top of my head if it appears in connection with (weakly holomorphic) modular forms. But I know there are some nice techniques about multiplicative functions that decrease over the primes, etc...