Show that : $\sum\limits_{\sigma \in S_n} (\mbox{number of fixed points of } \sigma)^2= 2 n!$

I came across this result while doing some representation theory of the permutation group $S_n$ $$ \sum\limits_{\sigma \in S_n} (\mbox{number of fixed points of } \sigma)^2 = 2 n!$$

This can be proved very easily using the permutation representation of $S_n$. Does anyone know of a direct, more elementary way of proving this without using any representation theory?

You can use Burnside's lemma. $\text{Fix}(\pi)^2$ is the number of elements fixed by $\pi$ acting on $[n]^2$, where $[n] = \{ 1, 2, ... n \}$, so Burnside's lemma tells you that

$$\frac{1}{n!} \sum_{\pi \in S_n} \text{Fix}(\pi)^2$$

is the number of orbits of $S_n$ acting on $[n]^2$. But there are, by inspection, two such orbits: the ordered pairs $(a, b)$ where $a \neq b$, and the ordered pairs $(a, a)$.

(Exercise: generalize this argument to figure out $\displaystyle \frac{1}{n!} \sum_{\pi \in S_n} \text{Fix}(\pi)^k$.)

We can use basic probability theory. Take a random permutation. Let $X_i=1$ if $i$ is a fixed point of the permutation, and let $X_i=0$ otherwise. Then the number of fixed points is $X_1+\cdots+X_n$. Square. We get $$\sum_i X_i^2+\sum_{(i,j), i\ne j} X_iX_j.$$ Now by the linearity of expectation, $$E((X_1+\cdots +X_n)^2)=\sum_iE(X_i^2)+\sum_{(i,j), i\ne j}E(X_iX_j).\tag{1}$$ We have $\Pr(Xi=1)=\frac{1}{n}$ and $\Pr(X_iX_j=1)=\frac{1}{n(n-1)}$.

It follows that $E((X_1+\cdots+X_n)^2)=2$. For our sum, multiply by the number $n!$ of permutations.

Remark: Another instance of a mean proof.

Here is an answer using Derangements. Purely combinatoric.

For any set of $k$ elements there are $\mathcal{D}(n-k)$ permutations that fix those $k$ elements and derange the others. There are $\binom{n}{k}$ ways to choose the $k$ elements. So the sum of the square of the number of fixed points is $$ \begin{align} \sum_{k=0}^nk^2\binom{n}{k}\mathcal{D}(n-k) &=\sum_{k=0}^nk(k-1)\binom{n}{k}\mathcal{D}(n-k)\\ &+\sum_{k=0}^nk\binom{n}{k}\mathcal{D}(n-k)\\ &=\sum_{k=0}^nn(n-1)\binom{n-2}{k-2}\mathcal{D}((n-2)-(k-2))\\ &+\sum_{k=0}^nn\binom{n-1}{k-1}\mathcal{D}((n-1)-(k-1))\\[6pt] &=n(n-1)(n-2)!\\[12pt] &+n(n-1)!\\[12pt] &=2n!\tag{1} \end{align} $$ for $n\ge2$. We've used formula $(2)$ from the cited answer: $$ n!=\sum_{k=0}^n\binom{n}{k}\mathcal{D}(n-k)\tag{2} $$ For $n=1$, $n(n-1)(n-2)!+n(n-1)!=1$ and for $n=0$, the sum is $0$.

Generalization

Qiaochu asked a question that I had actually thought about, but had missed a fine point: How does this generalize for summing other powers of the number of fixed points? For this, we will use Stirling Numbers of the Second Kind (whose defining equation is below in red) and equation $(2)$ above. $$ \newcommand{\stirtwo}[2]{\left\{{#1}\atop{#2}\right\}} \begin{align} \sum_{k=0}^n\color{#C00000}{k^p}\binom{n}{k}\mathcal{D}(n-k) &=\sum_{k=0}^n\color{#C00000}{\sum_{j=0}^p\stirtwo{p}{j}\binom{k}{j}j!}\binom{n}{k}\mathcal{D}(n-k)\\ &=\color{#00A000}{\sum_{k=0}^n}\sum_{j=0}^p\stirtwo{p}{j}\binom{n}{j}j!\color{#00A000}{\binom{n-j}{k-j}\mathcal{D}((n-j)-(k-j))}\\ &=\sum_{j=0}^p\stirtwo{p}{j}\binom{n}{j}j!\color{#00A000}{(n-j)!}\\ &=\sum_{j=0}^n\stirtwo{p}{j}n!\tag{3} \end{align} $$ For $n\ge p$, $(3)$ can be simplified to $$ \sum_{j=0}^n\stirtwo{p}{j}n!=\mathrm{B}_pn!\tag{4} $$ where $\mathrm{B}_p=\displaystyle\sum_{j=0}^p\stirtwo{p}{j}$ are the Bell Numbers.

The labelled species $\mathcal{Q}$ of permutations with fixed points marked that we are working with here is $$\mathfrak{P}\left(\mathcal{UZ} + \mathfrak{C}_2(\mathcal{Z}) + \mathfrak{C}_3(\mathcal{Z}) + \mathfrak{C}_4(\mathcal{Z})+\cdots\right).$$

Therefore the bivariate generating function corresponding to $\mathcal{Q}$ is given by $$G(z, u) = \exp\left(uz-z+\log\frac{1}{1-z}\right) = \frac{e^{-z}}{1-z} e^{uz}.$$ We want to turn the terms having shape $q\times u^k z^n/n!$ into $q\times k^2 z^n/n!$, so differentiate with respect to $u$ to get $$\frac{d}{du} G(z,u) = \frac{e^{-z}}{1-z} z e^{uz}$$ and multiply by $u$ for $$\left(u\frac{d}{du}\right) G(z,u) = \frac{e^{-z}}{1-z} uz e^{uz}.$$ Differentiate with respect to $u$ one more time to get $$\frac{d}{du}\left(u\frac{d}{du}\right) G(z, u) =\frac{e^{-z}}{1-z} \left(z e^{uz} + uz^2 e^{uz}\right)$$ and finally put $u=1$ for the end result $$\left.\frac{d}{du}\left(u\frac{d}{du}\right) G(z, u) \right|_{u=1} = \frac{e^{-z}}{1-z} \left(z e^{z} + z^2 e^{z}\right) = \frac{z+z^2}{1-z}.$$ Now for $n\ge 2$ we have $$n! [z^n] \frac{z+z^2}{1-z}= 2\times n!,$$ done (the value at $n=1$ is correct also). The collection of proofs on this page and their diversity is remarkable. A closely related species is studied at this recent MSE link.

Addendum. Here is a proof of the generalization by @robjohn. Observe that $$\left.\left(\frac{d}{du}\right)^k G(z, u)\right|_{u=1}$$ is the generating function of the factorial moments of the RV $F$ representing the fixed points, e.g. $$\left.\left(\frac{d}{du}\right)^3 G(z, u)\right|_{u=1} = \sum_{n\ge 0} \mathrm{E}[F_n(F_n-1)(F_n-2)] z^n.$$ Now recall that $$x^q = \sum_{k=0}^q {q\brace k} x^{\underline{k}}.$$ This immediately implies that $$\mathrm{E}[F_n^q] = [z^n] \sum_{k=0}^q {q\brace k} \left.\left(\frac{d}{du}\right)^k G(z, u)\right|_{u=1} = [z^n] \sum_{k=0}^q {q\brace k}\left. \frac{e^{-z}}{1-z} z^k e^{uz}\right|_{u=1} \\ = [z^n] \sum_{k=0}^q {q\brace k} \frac{z^k}{1-z}.$$ This means that for $n\ge q$ we have $$\mathrm{E}[F_n^q] = \sum_{k=0}^q {q\brace k} = B_q$$ as claimed.