**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}.$$