This is not homework.
My question is:
Prove or disprove: $$\lim_{n \to +\infty} \frac{1}{n^2} \sum_{a,b=1}^n \frac{1}{ \mathrm{gcd} (a,b)} = \frac{\zeta(3)}{\zeta(2)}$$
This would represent the probability as $n \to +\infty$ that, after having picked randomly $a,b,c \in \{ 1, \dots ,n \}$, the line $$ax+by+c=0$$ has some integer coordinate point $(x,y \in \Bbb Z^2)$. Indeed, fixed $a,b \in \{ 1, \dots ,n \}$, then $c=-ax-by$ for some $x,y \in \Bbb Z$ is equivalent to $$c \mathrm{\ is \ a\ multiple\ of\ gcd}(a,b)$$ because of Bezout identity; this event happens with probability $1/\gcd (a,b)$.
What I tried first: clearly $$0 \le \frac{1}{n^2} \sum_{a,b=1}^n \frac{1}{ \mathrm{gcd} (a,b)} \le \frac{1}{n^2} \sum_{a,b=1}^n 1 = 1$$ so, our limit belongs to $[0,1]$ (if it exists).
Then I made some euristic argument:
I know that the probability, as $n \to +\infty$, that two random numbers $a,b\in \{ 1, \dots ,n \}$ are coprime is $$\zeta(2)^{-1} = \frac{6}{\pi^2},$$ thus the probability that $\gcd(a,b)=d$ is equal to $$\frac{\zeta(2)^{-1}}{d^2} = \frac{6}{d^2\pi^2}$$ Thus our probability is (and here is where I have doubts in my argument) $$\sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2} \to \frac{1}{\zeta(2)} \cdot \zeta(2) =1$$
EDIT: or maybe my last step should be $$\sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2} \frac{1}{d} \to \frac{\zeta(3)}{\zeta(2)}$$
There is no need to involve probabilities, just a bit of analytic number theory suffices. I don't know how familiar you are with $L$-functions, so I'll leave some details open for now. If any part needs clarification, let me know.
For every positive integer $n$ define $$f(n):=\sum_{a,b=1}^n\frac{1}{\gcd(a,b)} \qquad\text{ and }\qquad g(n):=\sum_{d\mid n}\frac{\varphi(\tfrac{n}{d})}{d}$$ so that for all $n>1$ we have \begin{eqnarray*} f(n)-f(n-1) &=&-\frac1n+2\sum_{c=1}^n\frac{1}{\gcd(c,n)}\\ &=&-\frac1n+2\sum_{d\mid n}\frac{\varphi(\tfrac{n}{d})}{d}\\ &=&-\frac1n+2g(n). \end{eqnarray*} Because $\lim_{n\to\infty}\frac{1}{n^2}\sum_{k=1}^n-\frac1k=0$, this already shows that $$\lim_{n\to\infty}\frac{1}{n^2}f(n)=2\lim_{n\to\infty}\frac{1}{n^2}\sum_{k=1}^ng(n).$$ Note that $g$ is a multiplicative function, and in fact $g=\iota\ast\varphi$ where $\iota(n):=\tfrac1n$ and $\ast$ denotes the Dirichlet convolution of multiplicative functions. This immediately shows that the $L$-function associated to $g$ converges for all $s\in\Bbb{C}$ with $\operatorname{Re}(s)>2$ and that for all such $s\in\Bbb{C}$ we have $$L_g(s)=L_{\iota}(s)L_{\varphi}(s)=\zeta(s+1)L_{\varphi}(s).$$ It is easy to see that $L_{\iota}(s)=\zeta(s+1)$ for these $s$, and it is a basic exercise to show that $$L_{\varphi}(s)=\frac{\zeta(s-1)}{\zeta(s)},$$ whenever $\operatorname{Re}(s)>2$. Then by a Tauberian theorem, for example Wiener-Ikehara, we have \begin{eqnarray*} \lim_{n\to\infty}\frac{1}{n^2}f(n) &=&2\operatorname{Res}(L_g,2)\\ &=&2\cdot\operatorname{Res}(\zeta(s+1),2)\cdot\operatorname{Res}(L_{\varphi},2)\\ &=&2\cdot\zeta(3)\cdot\frac{1}{2\zeta(2)}\\ &=&\frac{\zeta(3)}{\zeta(2)}. \end{eqnarray*}