What is the probability two random maps on n symbols commute? It is well known that two randomly chosen permutations of $n$ symbols commute with probability $p_n/n!$ where $p_n$ is the number of partitions of $n$. This is a special case of the fact that in a group, the probability that two elements chosen uniformly at random (with repetition allowed) commute is the number of conjugacy classes divided by the size of the group.


Question. What is the probability that two mappings of $n$ symbols chosen uniformly at random commute?


I suspect an exact answer would be difficult and would be happy to learn of reasonably tight asymptotic results.
Added. This probability should go to zero quickly because Misha Berlinkov recently showed that with probability going to 1 as $n$ goes to infinity, two random elements generate a subsemigroup containing a constant map and so if they commute they generate a unique constant map. This should happen almost never (and most likely has been proven).
Added based on Brendan McKay's answer. Computing the probability that an element of a monoid $M$ commutes with an element of its groups of units $G$ is no harder than the commuting probability in a group.  Namely, $G$ acts on $M$ by conjugation; let's call the orbits conjugacy classes.  Then the probability that an element of $G$ commutes with an element of $M$ is the number of conjugacy classes of $M$ divided by the number of elements of $M$. The proof is the same as for groups.  If $Fix(g)=\{m\in M\mid gmg^{-1}=m\}$, then
$$\frac{|\{(g,m)\in G\times M\mid gm=mg\}|}{|G||M|} = \frac{1}{|M|}\frac{1}{|G|}\sum_{g\in G}|Fix(g)| = \frac{\text{number of conjugacy classes}}{|M|}$$ by the Cauchy-Burnside-Frobenius orbit counting formula.
For $M=T_n$ the monoid of all mappings on $n$ symbols and $G=S_n$ the symmetric group, conjugacy classes correspond to isomorphism classes of functional digraphs on $n$ vertices.  A functional digraph is a digraph (loops allowed) in which each vertex has outdegree $1$.  Each mapping $f$ gives a functional digraph by drawing an edge from $i$ to $f(i)$. It is obvious that $f,g$ are conjugate iff their corresponding digraphs are isomorphic (it is the same proof that permutations are determined up to conjugacy by cycle type).
According to the book of Flajolet and Sedgewick, the number of unlabelled functional digraphs grows likes $O(\rho^{−n}n^{−1/2})$ where $\rho\approx .29224$. So the probability of a random mapping commuting with a random permutation is pretty small. Brendan raises the nice question of how different the probability of a random permutation commuting with a random mapping is from the probability of a random mapping commuting with a random mapping.  My guess is the latter goes to $0$ qualitatively faster.
 A: The probability is $e^{-2n (\log \log n + 1 - o(1))}$. $\DeclareMathOperator{\kernel}{kernel} \DeclareMathOperator{\range}{range}$
Lower Bound
The above bound is achieved by random $f$ and $g$ with $\range(f) ⊆ \operatorname{kernel}(g)$ and $\range(g) ⊆ \kernel(f)$ (i.e. $∀x \, f(g(x))=g(f(x))=0$); typically for such $f$ and $g$, $|\kernel(f)| ≈ |\kernel(g)| ≈ n/\log n$ (up to $1±o(1)$ factor).
Commutativity is extremely unlikely absent special patterns, and (as we show below) the most likely is (approximately) the above annihilation pattern.
Besides choosing the value of "zero", the one change I have to improve the probability is to have $≈\log^2 n$ exceptions for which commutativity happens accidently:  For a typical input $x$, there are $≈n^2/\log^2 n$ ordinary choices of $(f(x),g(x))$, and $≈n$ choices for accidental commutativity (i.e. $f(g(x))≠0$).  I suspect that this model captures $1-o(1)$ fraction of commuting $f$ and $g$.
With this model, we can compute the asymptotics quite precisely, though it is easy to make a mistake about a constant or otherwise small factor.  I came up with the following probability, where $h(p) = -p \log p - (1-p) \log(1-p)$ (entropy for probability $p$):
$(1±o(1)) \max\limits_{n^{-1/3}≤p<1} p^3 n e^{ -2n(p \log n - (1-p) \log p - h(p)) \, + \, (1-p)^4 p^{-2} \, - \, 2(1/p-1)}$
Intuition for $p$: $p ≈ |\kernel(f)|/n ≈ |\kernel(g)|/n ≈ \log^{-1} n$.
Explanation of all factors:
$n$ - number of choices for "zero"
$p = (p^{1/2})^2$ - stems from the error tolerance (stdev) for $p$ (computed using the second derivatives) being $≈p^{-1/2}$ times smaller than $\sqrt{p/n}$; the square is because the correction is for both $f$ and $g$ (and the covariance of the errors is small enough).
$p^2$ - both kernels include 0
$e^{-2(1/p-1)} = n^{-2±o(1)}$ - probability that $f$ and $g$ are never 0 outside of what we designated as their kernels
$e^{-2n p \log n}$ - kernel
$e^{-2n(1-p) \log p}$ - range restriction
$e^{2h(p)n}$ - effective number of choices for the kernels (before the $p^3$ correction)
$e^{(1-p)^4 p^{-2}} ≈ (1 + \frac{(n(1-p))^2/n}{(pn)^2})^{n(1-p)^2} = n^{(1±o(1))\log n}$ is the correction for the exceptions with accidental commutativity; $≈(1-p)^2$ is the probability of being outside both kernels; $(n(1-p))^2/n$ is accidental commutativity ($n(1-p)$ choices each times $n^{-1}$ collision chance); $n^{-1/3}≤p$ (chosen arbitrarily) more than suffices for the accuracy of the approximation by the exponential.
Numerical results: For $n≤20$, the exact values (times $n^{2n}$) are published as an OEIS sequence A181162 and are in an order-of-magnitude agreement with the predicted values (the above formula underestimates in about 4-5 times there, but (for $n>10$) improving with $n$).  This testing is limited because $p$ goes to zero very slowly.
Upper Bound
For a mapping $f$, the probability of commutativity with a random mapping $g$ is at most $\min_{S⊆\operatorname{domain}(f)} n^{-|f[S]  \setminus S|}$ since for commuting $f$ and $g$,   $g\restriction (f[S] \setminus S)$ is uniquely determined by $f$ and $g \restriction S$.
Thus, for random commuting mappings, with high probability, $|\range(f)|$ and $|\range(g)|$ are both $O(n / \log n ⋅ \log \log n)$ (we will use this below).  This range restriction has probability $e^{-(2±o(1)) n \log \log n}$.  If $|\range(f)|$ were too large here, then unless $f$ is identity on too many points to be likely, we get a large $|f[S] \setminus S|$ and thus commute probabilities well below the lower bound.
We now get a more precise upper bound.  For $x∈\range(g)$ define the equivalence class $[x]_f = \{y∈\range(g): f(x)=f(y)\}$, and analogously define $[x]_g$ (for $x∈\range(f)$).  If $f$ and $g$ commute, we get a corresponding bijection $P$ between the set of all $[x]_f$ and the set of all $[x]_g$: $P([g(x)]_f) = [f(x)]_g$.  Let $S = \{x∈\range(g): |[x]_f|+|P([x]_f)|≥ne^{-\sqrt{\log n}}\}$, and set $T = P[S] = ∪\{P([x]_f):x∈S\}$ (here, $e^{-\sqrt{\log n}}$ can be arbitrary within bounds).  Let $m = |\range(f)∪\range(g)|$.
The probability that random $f$ and $g$ commute is at most $\sup_{|S|≤m, \, |T|≤m, \, m=O(n / \log n ⋅ \log \log n)} e^{o(n)} n^{-|S|-|T|} (\frac{|S||T|}{n^2} + e^{-\sqrt{\log n}})^{n-m} = $ $e^{-2n (\log \log n + 1 - o(1))}$ as desired.  Here, $n^{-|S|-|T|}$ is the probability that $f$ and $g$ both have few enough values on $S$ and $T$ respectively (ignoring $e^{o(n)}$ factors).  Next, after setting $f$ and $g$ on inputs in $\range(f)∪\range(g)$, we test $f(g(x))=g(f(x))$ for $x$ outside of $\range(f)∪\range(g)$ (hence the $n-m$ exponent) with $f$ and $g$ random on $x$.  We either hit a large equivalence class (the $|S||T|$ factor) or deal with the low likelihood that $f$ and $g$ hit matching classes (hence the $e^{-\sqrt{\log n}}$ above).
Furthermore, the $|S||T|$ above is an overestimate unless there is a single dominant equivalence class:  The $|S||T|$ can be refined to $\sum_i |S_i| |P(S_i)|$ where $S_i⊆S$ is an equivalence class (and with a further reduction if $|S_i||P(S_i)|$ does not reflect the probability of $f(x)∈S_i$).  To match the $e^{-2n (\log \log n + 1 - o(1))}$ bound, $|S|$, $|T|$, and for $n \! - \! o(n)$ $x$, $|[g(x)]_f|$ and $|[f(x)]_g|$ must all be $(1±o(1))n/\log n$.
Thus, with high probability, for random commuting $f$ and $g$, the composition $fg$ is constant, except on $o(n)$ inputs.
Related Problems
As you note, the probability that a random permutation $f$ commutes with a random mapping $g$ is $Θ(n^{-n-1/2} ρ^{-n})$ for $ρ≈0.29224$.  This probability should be dominated by $f$ for which the number of $k$-cycles decreases roughly exponentially with $k$.  Here is a derivation of the weaker $e^{Θ(n)}/n!$ bound.
  Briefly, for commutativity, a $k$-cycle of $f$ must be mapped to a cycle whose length divides $k$, and the total number of commuting mappings $g$ is $\prod_k (\sum_{q|k} qc_q)^{c_k}$ where $c_k$ is the number of $k$-cycles of $f$ (so $\sum_k kc_k=n$).  We get the $e^{Θ(n)}/n!$ lower bound by considering cycles of length 1 and 2.  For the $e^{Θ(n)}/n!$ upper bound, if $f$ has $m$ cycles, the number of choices of $f$ is (at most) $e^{O(n)} n^{n-m}$, while the commute probability with a random $g$ is $≤n^{-(n-m)}$.
One can also ask about the likelihood of commutativity in the absence of special patterns:
- For almost all mappings $f$, I think the number of commuting mappings $g$ is $e^{Θ(n)}$:  Once we set $g(x)$ for a given $x$, we can go back through all values of $f^{-1}(x)$, $f^{-2}(x)$, ..., and typically have $O(1)$ choices for $g$ per point.
- For almost all mappings $f$, the number of commuting permutations $g$ is $e^{Θ(n)}$. For the lower bound, for every $x∉f[\range(f)]$, one can arbitrarily permute $f^{-1}(x)$.  The upper bound follows from the bound on the probability that a random permutation and a random mapping commute.
- For almost all permutations $f$, the number of commuting mappings $g$ should be $n^{(\frac{1}{2}±o(1)) \log n}$:  A random permutation $f$ has an expected $≈\log n$ cycles, with an average of $1/k$ cycles of length $k$; and a large cycle of $f$ should be typically mapped to itself by $g$.  Furthermore, for every constant $ε>0$, for $1-ε$ fraction of permutations $f$, a $Θ(1)$ fraction of commuting $g$ should be permutations.
A: The expected value in Brendan McKay's answer on the probability of $f(g(1))=g(f(1))$ is correct. Just count the quintuples $(a,b,c,f,g)$ where $g(1)=a$, $f(1)=b$, $f(a)=c$, $g(b)=c$. For instance, there are $n(n-1)^2$ triples $(a,b,c)$ with $a\ne1$, $b\ne1$, and for each such triple there are $n^{n-2}$ possibilities for $f$ and $g$ each, contributing $n(n-1)^2n^{2n-4}$ to the possibilities. In the cases $a\ne1$, $b=1$ we must have $c=a$, so the contribution is $(n-1)n^{2n-3}$. The same for $a=1$, $b\ne1$. Finally, if $a=b=1$, then $c=1$, and that case contributes $n^{2n-2}$.
A: Let $f$ and $g$ be random mappings. If they commute, then $f(g(1))=g(f(1))$, and this happens with probability $n^{-1}$. Now $f(g(2))=g(f(2))$ also holds with probability $n^{-1}$, but these events need not be independent. To make them independent, suppose that the first equation holds. Fix $f(1), g(1), f(g(1)), g(f(1))$, without loss assume that 2 is not among them. With probability $\leq\frac{16}{n}$ we have that one of these elements equals one of $f(2), g(2), f(g(2)), g(f(2))$, in which case things get messy. If not, then $f(g(2))=g(f(2))$ holds with probability $\frac{1}{n-|\{f(1), g(1), f(g(1)), g(f(1))\}|}<\frac{2}{n}$. Hence the second equation holds true with probability $\leq\frac{18}{n}$.
Now continue. In the $k$-th step the conditional probability that the element you pick refutes commutativity subject to the condition that the mappings satisfy all equations already fixed is $\frac{1}{n-\#\{f(1), \ldots, g(x_k)\}}\leq \frac{1}{n-4k}$, if things are nice. The probability that things are messy is $\leq\frac{16 k}{n}$. Hence the probability that $f$ and $g$ commute is bounded above by
$$
\prod_{16k<n}\left(\frac{1}{n-4k}+\frac{16 k}{n}\right)\leq e^{-n/18}.
$$
A: The number of ordered pairs of commuting functions is A181162.  I agree with those counts up to n=7.  There is little in OEIS that helps to answer the asymptotics question.
Incidentally, the probability that $f(g(1))=g(f(1))$ is not $1/n$.  I think it is $1/n + (n-1)/n^3~$ though I might have miscalculated.  That formula works up to $n=7$.
ANOTHER relevant fact:  If $f$ is a permutation, then any function $g$ commuting with $f$ is determined by the image of one element of each cycle of $f$.  So the number of such $g$ is at most $n^{C(f)}$ where $C(f)$ is the number of cycles of $f$. Random permutations have on average only $\ln n+O(1)$ cycles, so the probability of a random function commuting with a random permutation might be at most something like $n^{-n+\ln n+O(1)}$ (which is an abuse of expectations but might be something akin to the truth).  Is a random function more or less likely to commute with another random function or with a random permutation?  [NOTE: I added "at most" since some assignments don't work: the image of a point in a cycle of length $k$ must lie in a cycle whose length is a divisor of $k$.]
