Given a positive definite integral binary quadratic form $f$, denote by $r_f(n)=\#\{(x,y)\in \mathbb{Z}^2~|~f(x,y)=n\}$, how can we obtain estimates for $$\sum_{~~~~~n\leq x\\ n\equiv k~\mathrm{mod}~q}r_f^2(n),$$ where $(k,q)=1$.
I'm aware of the results of Tolev here, but he treats the special case $f(x,y)=x^2+y^2$ and studies $\sum r_f(n)$ instead. Is there away to do this in more generality? Maybe the work of Blomer and Granville might help, as I need an explicit dependence on the form $f$ in the main term.
I know $$\sum_{~~~~~n\leq x\\ n\equiv k~\mathrm{mod}~q}r_f^2(n)=\frac{1}{q}\sum_{n\leq x} r_f^2(n)+\frac{1}{q}\sum_{1\leq a\leq q-1}e(-ak/q)\sum_{n\leq x} r_f^2(n)e(an/q).$$ The first term can be estimated by a result due to Muller down here, but how do I show that the error term coming from the other additive characters is small enough? say of order $\ll x^{3/4}$? I believe we have to examine the Dirchlet function $$L(r_f^2,a,s)=\sum \frac{r_f^2(n)e(an/q)}{n^s}$$ where $a\in \mathbb{Z}/q\mathbb{Z}$.
Any help would be appreciated.
Edit: I found this paper of Muller (see theorem 6.1) where he studies $\sum_{n\leq x}r_f^2(n)$ for any quadratic form $f$ ( doesn't have to be necessarily of fundamental discriminant), but his approach uses the mean-squared estimate which eliminates the usage of the characters.