Here is the link to the OEIS sequence of the first several terms

This OEIS entry also contains two links to Math StackExchange discussions

I personally love such kinds of problems because $n^n$ denoting the total number of endofunctions (i.e. functions from a set to itself) is also the number of rooted trees with one marked vertices. There are nice combinatorial bijections which can be used for asymptotic analysis.

Let us start with an easier problem.

Clearly, for each of the $n$ points there are $(n-1)$ allowed possibilities of where the function $f$ would be pointing to. The number of functions is therefore $(n-1)^n$. The probability that a randomly chosen function does not have fixed points is asymptotically

In this case, all the fixed points are formed independently, which is not the case in our problem.

Furthermore, if we look at the exact formula from OEIS, we observe

which I guess could be obtained by some kind of inclusion-exclusion. However, asymptotic analysis of such formulae may be sometimes difficult, while can still be useful for numerical simulations. In this case, however, asymptotic analysis is possible because the summands do not give increasingly divergent contributions.

Let us take a closer look at each of the summands as $n \to \infty$.

Note that $\binom{n}{k}$ is asymptotically $n^k / k!$, i.e. a polynomial of degree $k$. The ratio $n! / (n-k)!$ is again a polynomial of degree $k$, so when divided by $n^{2k}$ the remaining part asymptotically does not depend on $n$. The first several terms of this sum are approximately

With some extra work to make this proof rigorous we arrive to the answer $e^{-1}$.

Source: https://math.stackexchange.com/a/983293/208535

Suppose that the items from the set $B$ receive, accordingly, $k_1, k_2, \ldots, k_n$ incoming edges.

We can notice that since $k_1 + \ldots + k_n = n$, we have $\mathbb E k_i = 1$ for each $i$. We could also compute $\mathbb E k_i^2$ and make sure it is finite (asymptotically does not depend on $n$). In fact each $k_i$ is distributed according to a Binomial distribution (although not independently) $\mathrm{Binom}(n,\frac{1}{n})$. Its variance is $\mathbb E k_1^2 - (\mathbb E k_1)^2 = 1 - \frac1n$, and therefore, $\mathbb E k_1^2 \sim 2$ as $n \to \infty$.

First of all, suppose that all $k_i$ are equal to 1. Then, the probability that there are no mutual friendships is the same as in the derangements problem

Now, in the case of general $k_i$, we replace the 1’s in the product by $k_i$’s:

By developing the Taylor series of the log, we obtain

This requires a better justification but heuristics says that on average, the answer does not depend on $k_i$ unless some of the $k_i$ is too large, which does not happen too often. By decomposing the probability up to conditional probabilities, we obtain

Let us solve the derangement problem in a hard way first. From each node there is an arrow to another node. We look at the ensemble of arrows and nodes as a whole and try to find some structure in it.

It turns out that some nodes assemble into cycles, and others attach to nodes assembled as cycles. Everything outside the cycles forms trees ascending to one of the nodes of a cycle.

With some help of the exponential generating functions, we can express the combinatorial property that the endofunctions are sets of these cycle-trees (unicycles).

where $U(z)$ is yet-to-be-determined exponential GF of unicycles. Since the exponential GF of cycles is just simply

the GF of unicycles could be obtained by a functional composition with the cycle function

where the generating function of rooted trees $T(z)$ is defined by a recursive relation “a rooted tree is composed of a root and a set of rooted trees”:

This gives a surprising identity

The above approach allows to easily obtain the celebrated formula for the number of unrooted labelled trees $T_n = n^{n-2}$. From the functional definition of the function $T(z)$, we get, by taking a pointed derivative

which is a linear equation that can be solved:

Therefore,

and it easily follows that the number of rooted trees is $n^{n-1}$, so therefore, the number of non-rooted trees is $n^{n-1}$ divided by $n$, which is $n^{n-2}$.

The form of the generating functions allows one to obtain the Poisson distribution for the number of isolated cycles with almost no extra effort. Since the GF for all cycles is

we can mark the cycles of length 1 by introducing an extra variable $x$:

The generating function of all the endofunctions with isolated cycles marked is now

The heuristics goes as follows:

But wait, we almost forgot about the factor $\dfrac{1}{1 - T(z)}$ in front and it is terrifying! But in order to get the probability generating function for the number of cycles, we can divide the counting GF by the total GF, just like in combinatorial probability, where we divide the number of “good” outcomes by the number of “total” outcomes.

and we know that the multiple $\dfrac{1}{1 - T(z)}$ will introduce some perturbation, but $e^{x}$ is a tiny guy (with asymptotically discrete distribution), so the final answer will not be changed. We, therefore, obtain, by some hand-waving

The form of the generating functions also allows one to get the asymptotic distribution of previously unaccessible patterns. For example, we could ask for the distribution of the number of unicycles in the above graph. The answer comes by marking the unicycles with an extra variable:

We could still get the expected number of the cycles by using this approach.

The asymptotics of this generating function can be obtained by applying transfer theorems of Flajolet-Odlyzko. This is more or less a routine procedure and requires only the knowledge of these generating functions around their respective singularities.

I invite the curious reader to check out a paper of Flajolet and Odlyzko entitled “Random mapping statistics” where they squeeze out everything that is possible to know about random mappings by a repeated application of the transfer theorem and the above symbolic method.

Let us return to the original problem. It seems to be more difficult than just treating one mapping because now we have two mappings. But we will see that there is no need in applying transfer theorems at all, and a simple approach will do.

Now we have vertices of two colours, and we have to make sure that the arrows go between alternating colours.

We will use a variable $z_1$ for the vertices of the first colour, and $z_2$ for the vertices of the second.

The generating function of trees rooted at, respectively, the first or the second colour, would be

They, still, share the same singularity, but now it is trickier because in two dimensions, the singularity is not a single point, but a manifold. Well, if we require $z_1 = z_2$ then the singularity would still be a point $z_1 = z_2 = e^{-1}$. In the end we need to make sure that the number of vertices of each colour is equal, which indirectly points to the symmetry and to the fact that $z_1 = z_2$ in the heuristics.

Since in every cycle we have an even number of vertices, and black and white trees alternate, we therefore, have the generating function of cycles

Excluding the cycles of length 2 would mean cropping the first term in the expansion of the logarithm:

and the composition of two mappings now has a generating function with $x$ marking all the cycles of size $2$, also corresponding to fixed points:

Note that the singular value of $T_1 T_2$ is again equal to 1 at the point $z_1 = z_2 = e^{-1}$.

As in the derangement exercise, we extract the coefficient $[z_1^n z_2^n]$ and this gives us the probability generating function of the number of fixed points: