For any set of numbers $1$ through $n$, what is the average probability that two randomly selected derangements are derangements of each other?

First, let's acknowledge that if $n=1$ then $P(A_1)=0$, since there are no derangements of a set of order one; and if $n=2$, $P(A_2)=0$, since there is only one derangement of a set of order two. For $n=3$, the derangements of the set {$1,2,3$} are {$2,3,1$} and {$3,1,2$}. Since these are the only examples, and the 1st, 2nd, and 3rd digits are all different, the probability so far is $P(A_3)=1$.

Now take into account $n=4$. The derangements now are

1. {$2,1,4,3$}
2. {$2,3,4,1$}
3. {$2,4,1,3$}
4. {$3,1,4,2$}
5. {$3,4,1,2$}
6. {$3,4,2,1$}
7. {$4,1,2,3$}
8. {$4,3,1,2$}
9. {$4,3,2,1$}

In this example, I think that the chance of picking two derangements that are derangements of each other is $P(A_4)=.5$, since, for example, if 1 is picked, then there is a $\frac{1}{2}$ chance of picking either 5, 6, 8, or 9 from the eight remaining derangements, and I believe the same is found if any other derangement is picked first.

We can define $P(A, n)=\frac{1}{n} \Sigma_{i=1}^n P(A_i)$. What can we expect $\lim_{n\rightarrow\infty}P(A,n)$ to be? If we cut off $n$ at $n=4$, then we get $P(A,4) = \frac{1}{4}(0+0+1+.5) = 0.375$. Can we expect $\lim_{n\rightarrow\infty}P(A,n)$ to converge to zero, or can we expect $\lim_{n\rightarrow\infty}P(A,n)>0$? If a non-zero value of the limit exists, what is it, and how can it be found?

Let $D_n \subset S_n$ be the set of derangements. It is well-known that $$|D_n| / |S_n| \to e^{-1}$$ as $n \to \infty$. In fact $|D_n|$ is exactly equal to $n!/e$ rounded to the nearest integer.

You ask about the probability that two elements $x, y \in D_n$ are derangements of each other, or equivalently $xy^{-1} \in D_n$. Now because $xy^{-1}$ has no apparent special properties one expects it to be roughly the same thing as a uniformly random element of $S_n$, and therefore one expects $P(xy^{-1} \in D_n) \to e^{-1}$.

This guess can be verified using some nontrivial tools from representation theory. In terms of convolution, we are trying to estimate $\langle 1_{D_n} * 1_{D_n}, 1_{D_n}\rangle$, and this can be expressed in terms of characters: $$\langle 1_{D_n} * 1_{D_n}, 1_{D_n} \rangle = \sum_\chi \chi(1)^{-1} \langle \chi, 1_{D_n} \rangle^3.$$ The main term, arising from the trivial character, is $\langle 1, 1_{D_n} \rangle^3 = (|D_n| / |S_n|)^3$. The term arising from the sign character is $\pm (n-1)^3 / n!^3$, which is tiny, and the sum of everything else is small because the degree of every other character is at least $n-1$ and $$ \sum_{\chi \neq 1, \mathrm{sgn}} \chi(1)^{-1} |\langle \chi, 1_{D_n}\rangle|^3 \leq (n-1)^{-1} \big(\sum_\chi |\langle \chi, 1_{D_n} \rangle|^2 \big)^{3/2} = (n-1)^{-1} (|D_n| / |S_n|)^{3/2}.$$ This is also related to quasirandomness.