Sums of two squares are strongly biased towards quadratic residues

Chebyshev famously observed empirically that more often than not, there are more primes of the form 3 mod 4 up to x than of the form 1 mod 4. This was confirmed theoretically much later by Rubinstein and Sarnak in a logarithmic density sense. Our understanding of this is conditional on the generalized Riemann hypothesis as well as on the linear independence of the zeros of L-functions.We investigate similar questions for sums of two squares in arithmetic progressions. We find a significantly stronger bias than in primes, which happens for almost all integers in a natural density sense. Because the bias is more pronounced, we do not need to assume linear independence of zeros, only a Chowla-type conjecture on nonvanishing of L-functions at 21. To illustrate, we have under GRH that the number of sums of two squares up to x that are 1 mod 3 is greater than those that are 2 mod 3 100% of the time in natural density sense.