Divisor sums over values of binary forms of primes Let $\tau$ be the divisor function, that is 
$$
\tau(n)=\sharp\{d \in \mathbb{N}, d|n\}.
$$
I was wondering if anyone has ever proved an asymptotic estimate
for the sum
$$S(x):=\sum_{p,q\leq x}\tau(p^2+q^2),$$
where the summation is taken over pairs of primes.
One obviously expects $$S(x)\sim c\frac{x^2}{\log x}$$ as $x \to \infty,$
where $c$ is a positive constant which is an infinite product of Euler factors.
This is based on the heuristic that each of the $\pi(x)^2$
terms present in $S(x)$ is approximated by a constant multiple 
of
$\log x$ on average.
Brun-Titchmarsch and Bombieri-Vinogradov 
can be used to prove the upper and the lower bound
$$ 
c\frac{x^2}{\log x} (\frac{1}{2}+o(1))
\leq S(x) \leq 
c\frac{x^2}{\log x} (2+o(1)),
$$ as $x\to \infty$ respectively.
But the question remains, $\textit{can we prove an asymptotic?}$
 A: An answer regarding the use of large sieve suggested by Lucia (too long for a comment). 
I guess her/his thought was along the following lines: 
The divisors $d$ of $p^2+q^2$ are of order $x^2$ and the hyperbola trick reduces to estimating sums of the form $$\sum_{d\leq x}\sum_{\lambda^2=-1 (d)} \sum_{\substack{p,q \leq x \\ p= \lambda q (d)}} 1. $$ Inserting multiplicative characters in the sum over primes we will be left with an error term coming from the primitive characters, which looks like 
$$
\sum_{d\leq x}\frac{1}{\phi(d)}
\sum^*_{\chi(d)}
\sum_{\lambda^2=-1(d)}\overline{\chi(\lambda)}
\ \Big|
\sum_{p\leq x}\chi(p)\Big|^2 ,$$ 
which behaves as
$$\frac{1}{x}
\sum_{d\leq x}\tau(d)\frac{d}{\phi(d)}
\sum^*_{\chi(d)}
\ \Big|
\sum_{p\leq x}\chi(p)\Big|^2.$$
Ideally we would like to show that this is $o(x^2/\log x)$.
However the large sieve in the form of
[Th.4,p.g. 160, multiplicative Davenport]
gives a bound for this quantity which is of order $O(x^2)$
and is therefore inadequate. Essentially this is the level of distribution problem.
What we described can show that the contribution of $d\leq x/(\log x)^A$ to the sums with primitive characters is indeed $o(x^2/\log x)$ for an appropriate value of $A>0$ but I am not sure whether the micro logarithmic savings coming from Hooley's-Delta function can be used to say something about the remaining range.
A: When computing e.g. an asymptotic for $\sum_{p\leq x}d(p-1)$ you would like to estimate the number of primes $p$ such that $n$ divides $p-1$ by the prime number theorem as $\sim\frac{x}{\varphi(n)\log x}$. We do not know GRH, so we can't use this estimate for all $n<x^{1/2-\epsilon}$, but we do have Bombieri-Vinogradov, so we can use this estimate for almost all $n$ in the relevant range.
This approach does not work here, since for $\sum d(p^2+q^2)$ we have to take $n$ as large as $x^{1-\epsilon}$, however, the distribution of $p^2+q^2$ is a lot nicer than the distribution of $p$. So instead of a non-trivial bound for $\max_a\left|\pi(x, a, n) - \frac{x}{\varphi(n)\log x}\right|$ it suffices to bound 
$$
\sum_{a^2+b^2\equiv 0\pmod{n}} \pi(x,a,n)\pi(x, b, n) - \frac{x^2}{\varphi(n)^2\log^2 x}\sum_{a^2+b^2\equiv 0\pmod{n}} 1
$$
for almost all $n$. Since we are now averaging over $a, b$, we should be able to avoid Bombieri-Vinogradov but use the Barban-Davenport-Halberstam theorem, which is applicable to all $n$ up to $x/\log^A x$, which suffices.
Filling in the detail will probably be quite some work, but I guess that two or three pages of Cauchy-Schwarz and character calculations should do the trick.
