What is the probability of randomly selecting $ n $ natural numbers, all pairwise coprime?

Question

It's known that the probability of selecting $ n $ natural numbers randomly and ending up with a greatest common divisor equal to one is $ \prod (1-p^{-n}) = 1/\zeta(n) $. However, a total GCD of 1 does not rule out any of the pairs among the set of $ n $ numbers sharing a common factor. What's the probability none of them share a common factor? (Since there's a possibility of "random selection" being ambiguous, let's take it to mean chosen with uniform probability from {$ 1,2,3,\dots, N$} as $N \to \infty $.)

Answer 1

This answer is on a heuristical level; I don't purport to rigorously prove the result.

The probability $\prod (1-p^{-n})$ can be derived from the fact that for sufficiently large $N$ the probability of each number being divisible by a given prime $p$ is $p^{-1}$, so the probability of all $n$ numbers being thus divisible (which would lead to the gcd not being $1$) is $p^{-n}$. Applying this reasoning to your case, we have that the probability for not more than one of the $n$ numbers being divisible by $p$ is the sum of the probabilities of zero or one of them being divisible by $p$, which is

$$\binom n0\left(\frac{p-1}{p}\right)^n+\binom n1\left(\frac{p-1}{p}\right)^{n-1}\frac1p=\left(1-\frac1p\right)^{n-1}\left(1+\frac{n-1}p\right)\;.$$

So we could expect the probability you're looking for to be given by

$$p_n=\prod_p\left(1-\frac1p\right)^{n-1}\left(1+\frac{n-1}p\right)\;.$$

I don't know whether it's possible to evaluate this Euler product explicitly. Numerically (both evaluating the Euler product and generating large pseudo-random numbers), I get

$$ \begin{array}{|c|c|} n&p_n\\ \hline 1&1\\ 2&0.6079\\ 3&0.2867\\ 4&0.1149\\ 5&0.0409\\ \end{array} $$

Of course $p_2=1/\zeta(2)$, since in that case the two conditions coincide.

P.S.: The Euler product can be evaluated more efficiently numerically if we rewrite it as

$$p_{n+1}=\prod_p\left(1-\frac1p\right)^n\left(1+\frac np\right)=\prod_p\left(1-\frac1{p^2}\right)^n\frac{1+\frac np}{\left(1+\frac 1p\right)^n}=\zeta(2)^{-n} \prod_p\frac{1+\frac np}{\left(1+\frac 1p\right)^n}\;,$$

since the remaining product converges more quickly. (Also it looks more like something that might have a closed form.)