What is the probability of randomly selecting $ n $ natural numbers, all pairwise coprime? 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 $.)
 A: 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.)
