Sum of $k-$th powers of numbers of fixed points of permutations of $\{1,2,\cdots, n\}$ This originates from a combinatorics competition problem: given a permutation $\sigma$ of $\{1,2,...,2013\}$, let $f(\sigma)$ be the number of fixed points of $\sigma$. Find
\begin{equation*}
\sum_{\sigma\in S} f(\sigma)^4
\end{equation*}
where $S$ is the set of all possible permutations.
The key idea is to express
\begin{equation*}
f(\sigma)^4=\sum_{1\leq x_1,x_2,x_3,x_4\leq 2013} g(\sigma,x_1,x_2,x_3,x_4)
\end{equation*}
in which $g=1$ if $x_1,x_2,x_3,x_4$ are all fixed points and $g=0$ otherwise. Then we can interchange the summations and discuss the number of nonzero $g$ values in $\sum_{\sigma\in S} g(\sigma,x_1,x_2,x_3,x_4)$ for each quadruple $(x_1,x_2,x_3,x_4)$. The final answer is $15(2013!)$.
I have never seen this idea before. Is there a name of this method and any reference? Thanks a lot.
 A: This is a disguised special case of the proof of Burnside's lemma. If $G$ is a finite group acting on a finite set $X$ and $\text{fix}(g)$ denotes the number of fixed points of $g \in G$ acting on $X$, then Burnside's lemma says that
$$\frac{1}{|G|} \sum_{g \in G} \text{fix}(g)$$
is the number of orbits of the action of $G$ on $X$. The proof (one of them anyway) proceeds by writing
$$\text{fix}(g) = \sum_{x \in X} \delta_{gx, x}$$
where $\delta_{gx, x} = 1$ if $g$ fixes $x$ and $0$ otherwise. Then we interchange the sums, giving
$$\frac{1}{|G|} \sum_{g \in G} \sum_{x \in X} \delta_{gx, x} = \frac{1}{|G|} \sum_{x \in X} \text{stab}(x)$$
where $\text{stab}(x)$ here denotes the size of the stabilizer subgroup of $x$. If we organize the sum into orbits of the action of $G$ on $X$ we get
$$\frac{1}{|G|} \sum_{[x] \in X/G} \text{stab}(x) \text{orb}(x)$$
where $\text{orb}(x)$ denotes the size of the orbit containing $x$. By the orbit-stabilizer theorem we have $\text{stab}(x) \text{orb}(x) = |G|$ and the conclusion follows.

For application to your problem we want to take powers of $\text{fix}(g)$ but $\text{fix}(g)^k$ is just the number of fixed points of $G$ acting diagonally on $X^k$, so we get that
$$\frac{1}{|G|} \sum_{g \in G} \text{fix}(g)^k$$
is the number of orbits of the action of $G$ on $X^k$. Applied to the action of the symmetric group $G = S_n$ on $X = \{ 1, 2, \dots n \}$, we get that
$$\frac{1}{n!} \sum_{g \in S_n} \text{fix}(g)^k$$
is the number of orbits of the diagonal action of $S_n$ on $\{ 1, 2, \dots n \}^k$, which it's not hard to see is equal to the number of partitions of $\{ 1, 2, \dots k \}$ into at most $n$ nonempty subsets. This can be written as a certain sum of Stirling numbers, and if $n \ge k$ (which is the case here: $n = 2013, k = 4$) it's equal to the Bell number $B_k$ counting the number of partitions of $\{ 1, 2 \dots k \}$ into subsets. And indeed we have $B_4 = 15$ as desired.
This result previously occurred in this old math.SE answer of mine. It implies that the number of fixed points of a random permutation in $S_n$ has the same first $n$ moments as a Poisson random variable with $\lambda = 1$, and in fact you can prove that as $n \to \infty$ the number of fixed points of a random permutation converges in a suitable sense to such a Poisson random variable!
