What is the probability that a natural number is a sum of two squares?

Without going into detail what a "random natural number" might be, we could consider the density of such numbers.

The possibility to express $n$ depends on the prime factorization of $n$: It may be divisible by $2$ and by primes $\equiv 1\pmod 4$ as much as it likes, but each prime $\equiv 3\pmod 4$ must occur in an even power.

The prime $3$ "spoils" all numbers divisible by $3$ (that's $\frac13$), except those divisible by $9$ (that's $\frac19$), though it does spoil those divisible by $27$, execept ... All in all the prime $3$ spoils $\frac13-\frac19+\frac1{27}-\frac1{81}\pm\ldots = \frac14$ (geometric series) of all numbers in a large range. In general, a prime $p\equiv 3\pmod 4$ spoils $\frac1{p+1}$ and the effects of distinct primes are independent. Hence the limit density of natural numbers expressible as sum of squares is $$ \prod_{p\equiv 3\pmod 4}\left(1-\frac1{p+1}\right)$$ This product however diverges to $0$ because it is well-known that $\sum_p\frac1p$ diverges (where it doesn't matter that we use only half of all primes).

Thus: For sufficiently large $N$, the probability that a number picked uniformly from $\{1,\ldots,N\}$ is the sum of two squares becomes arbitrarily small. It is for any $\let\epsilon\varepsilon\epsilon>0$ we can find $N_0$ such that for all $N>N_0$ said probability is $<\epsilon$.


Though this is not an answer to the original question, I'll write it in answer to VividD's question under Hagen von Eitzen post (and I believe it is not totally unrelated to the question).

Let $A_n$ be the number of different pairs $(x,y)$ of non-negative integers solving the equation $x^2 + y^2 = n$. Then it is natural to define the expected number of compositions of a randomly chosen natural number into sum of two squares as
$$ A = \lim_{n\to\infty} \frac1n \sum_{k=1}^n A_k. $$ But the latter quantity is $\frac1n$ times the number of points of the form $(x/\sqrt{n},y/\sqrt{n})$, where $x,y$ are integers, inside the non-negative quarter of the unit circle centered at the origin (minus one, but this does not matter). So here we have a Jordan approximation of the area of the quarter. Therefore, $A = \frac\pi4$.

So, quite paradoxically, the expected value of the number of compositions of a randomly chosen non-negative integer into two squares is positive despite this number is zero "with probability 1".


Consider the function $r: \mathbb{N} \to \mathbb{C}$: $$ r(n) = \begin{cases} 1 &n \text{ is the sum of two squares (0 is a square)} \\ 0 &\text{otherwise}. \end{cases} $$ One way to define the "probability that a random integer is the sum of two squares" would be to consider the distribution on the integers where $n$ selected with probability proportional to $n^x$, for $x > 1$, and then to take the limit as $x \to 1$. That is, we can try to compute: $$ \lim_{x \to 1} \frac{\sum_{i=1}^\infty \frac{r(n)}{n^x}}{\sum_{i=1}^\infty \frac{1}{n^x}} \tag{1}. $$

As mentioned here, $r(n)$ is $1$ iff every prime of the form $4k-1$ occurs an even number of times in $n$. It follows that $r(n)$ is multiplicative, and \begin{align*} \frac{\sum_{i=1}^\infty \frac{r(n)}{n^x}}{\sum_{i=1}^\infty \frac{1}{n^x}} &= \frac{\prod_{p \text{ prime}} \sum_{i=0}^\infty \frac{r(p^i)}{p^{ix}}}{\prod_{p \text{ prime}} \sum_{i=0}^\infty \frac{1}{p^{ix}}} \\ &= \prod_{p \text{ prime}} \frac{1 + \frac{r(p)}{p^x} + \frac{r(p^2)}{p^{2x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1 + 0 + \frac{1}{p^{2x}} + 0 + \frac{1}{p^{4x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \prod_{p \equiv 1,2 \pmod{4}} \frac{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1 + \frac{1}{p^{2x}} + \frac{1}{p^{4x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1 \big/ \left(1 - \tfrac{1}{p^{2x}}\right)}{1 \big/ \left(1 - \tfrac{1}{p^{x}}\right)} \\ &= \prod_{p \equiv 3 \pmod{4}} \frac{1}{1 + p^{-x}}. \\ \end{align*} Now if you plug in $x = 1$, you get the product $$ \prod_{p \equiv 3 \pmod{4}} \frac{p}{p + 1} = \prod_{p \equiv 3 \pmod{4}} \left( 1 - \frac{1}{p + 1} \right) = 0, $$ by the reasoning of Hagen von Eitzen. So by this definition of probability, the probability that a random integer is a sum of two squares is zero.