Repeated draws from multinomial distribution (This is a cross-post from Math StackExchange https://math.stackexchange.com/questions/609641/multinomial-distribution-sum-of-squared-probabilities)
Let $\vec X = (X_1, \dots, X_k)$ be a draw from a fair multinomial distribution with $n$ trials, i.e. $P(X_1 = x_1, \dots, X_k = x_k) = \binom{n}{x_1, \dots, x_k} k^{-n}$
Let $\vec Y$ be an independent draw from the same distribution, i.e. $\vec X$ and $\vec Y$ follow the same law.
My questions is this: what is the probability that $\vec X = \vec Y$?
It appears that this probability, for fixed $k$, is roughly $(c_k + o(1)) n^{-(k-1)/2}$. This makes sense; the first $k-1$ coordinates vary on the order of $\sqrt{n}$ and are "independent-ish". (The last coordinate is of course fixed by the first $k-1$). Can one compute $c_k$?
 A: This is essentially a product of local CLTs (since one can uncover the multinomial
distribution variable by variable: decide first how many of the $n$ variables equal "1" with probability of success $1/k$ for each one and thus variance $(n/k)(1-1/k)$; Then, of the roughly $n(k-1)/k$ remaining variables, decide how many are "2" with probability of success 
$1/(k-1)$ each, and variance $(n/k)(1-1/(k-1))$, and so on). 
** Edit: As pointed out by Will Nelson,
there is a  factor of $1/\sqrt{2}$ in each step that comes from the fact that 
we are using the local CLT not at the mean but rather at the location of the
first multinomial $X$.
Thus, the probability is asymptotic to
$$\left(\frac{k}{4\pi n}\right)^{(k-1)/2}\cdot \frac{1}{\sqrt{\prod_{j=2}^k (1-j^{-1})}}$$ 
For $k=2$, this is just 
coming from 
$$\left(\frac{1}{\sqrt{2\pi n\cdot (1/4)}}\right)^2\int e^{-2\cdot (x^2/(2\cdot (1/4))}dx=1/\sqrt{\pi n}$$ 
A: I roughly outline a computation for $k=2$ below. Assuming I have not made a mistake for $k=2$, generalization to $k>2$ will likely be straightforward. If there seems to be demand, I plan on posting a proof of at least the $k=2$ case when I get a chance.
Direct computer simulation suggests $c_2 \approx 0.564$, which is roughly $\frac{1}{\sqrt{\pi}}$.
Here is the computation (again, this isn't a proof, but it's in the right direction):
\begin{eqnarray}
\sqrt{n}\sum_{i=0}^n \frac{1}{2^{2n}} \left(\begin{array}{c} n \\ i\end{array}\right)^2 
 &=& \frac{1}{2\pi\sqrt{n}}\left( \sum_{i=0\ ; \ \ x\equiv i/n}^n \frac{1}{x (1-x)} e^{-2 n (x\log x + (1-x) \log(1-x) + \log 2)} + O(n^{-1})\right) \\
 &=& \frac{1}{2\pi\sqrt{n}}\left( \sum_{i=0\ ; \ \ x\equiv i/n}^n 4 e^{-4 n (x-\frac{1}{2})^2} + O(n^{-1}) \right) \\
 &\to& \frac{4}{2\pi} \int_{-\infty}^{\infty} e^{-4 x^2} \ dx \\
 &=& \frac{1}{\sqrt{\pi}}.
\end{eqnarray}
The first line uses Stirling's approximation. The second uses a Taylor expansion at $x=\frac{1}{2}$. The third utilizes a straightforward limit of an approximation of the Riemann integral given on the third line.

For $k>2$, the same direct calculation approach yields:
$$
c_k = \frac{k^{\frac{k}{2}}}{(4\pi)^{\frac{k-1}{2}}}.
$$
Not surprisingly, this is the same result obtained by @oferzeitouni (though simplified). The local CLT approach may ultimately be simpler, but the straight analytic proof from Stirling's formula and a Taylor expansion is not difficult.
A: Write the probability as $k^{-2n}Q(k,n)$, where $Q(k,n)$ is the sum of the squares of the multinomial coefficients.  With the help of OEIS we find that these have been studied before:


*

*k=2: http://oeis.org/A000984 

*k=3: http://oeis.org/A002893

*k=4: http://oeis.org/A002895

*k=5: http://oeis.org/A169714

*k=6: http://oeis.org/A169715
The first few links give single summations and asymptotic formulas.  The last link says that $Q(k,n)$ is something called $W_k(2n)$ for which several references are given, and a generating function which I am guess-generalising here:
$$ \sum_{n=0}^\infty \frac{Q(k,n)}{n!^2}x^n
   = \left( \sum_{j=0}^\infty \frac{x^j}{j!^2} \right)^k. $$
I didn't prove this but I bet it is proved in one of the references given in the last link.  (Actually I think it is elementary: equate coefficients on each side with the original formulation of the problem in mind.) It looks like the ants' pants for deriving rigorous asymptotics.  The thing inside the large parens is of course a Bessel function.
This also gives a recurrence:
$$ Q(k,n) = n!^2 \sum_{j=0}^n \frac{Q(k-1,n-j)}{j!^2 (n-j)!^2}, $$
starting with $Q(0,n)=\delta_{0,n}$.
