Asymptotic independence in a multinomial setting. Let $(X_1,\ldots,X_r)$ be a multinomial vector with parameters $n$ and $1/r$, i.e., we throw $n$ balls into $r$ bins, with a uniform probability for each ball to land in each of the bins.  As is well known, $X_i$ and $X_j$ are dependent, with $\mathrm{cov}(X_i, X_j) = -n/r^2$.
Now let $r = r(n) \rightarrow \infty$, so we have a sequence of multinomial vectors (of increasing length), the nth vector being $(X_{1,n},\ldots,X_{r(n),n})$.  For a fixed index $k$, I know how to prove that if $r(n)/n \rightarrow 0$ as $n \rightarrow \infty$, then the normalized binomial RVs $[X_{k,n} - n/r(n)]/\sqrt{n(1/r(n))(1 - 1/r(n))}$ converge in distribution to a $N(0,1)$ RV (using the Lindeberg-Feller Theorem for triangular arrays).
However, I want to show that the resulting limiting RVs, for different $k$'s, are independent of each other.  The right way to think about it, I believe, is to extend each of the multinomial vectors with infinitely many zeros (say) to the right, i.e., to define $X_{i,n} = 0$ for $i > r(n)$; we end up with a a sequence of well-defined processes $X_1, X_2,\ldots$, where $X_n = \{X_{1,n}, X_{2,n}, \ldots\}$.  To show independence between the elements of the limiting process, it is enough to show that for each fixed $k$, the joint distribution of the first $k$ elements of the processes (properly normalized) converge in distribution to $k$ independent standard normal RVs.
But how do I go about this?  Simply proving that the limiting covariance is zero is not enough, I think, since zero covariance implies independence only in a multivariate normal setting, and I haven't proved that the limiting RVs have jointly a multivariate normal distribution.  The problem seems to me too elementary to be new, so I will be grateful for any references for an existing proof, or for ideas how to proceed.  Thanks.
 A: The simplest way is to couple your process with something obvious. It can be done in many ways. I would do the following.
Consider $r$ bins represented as disjoint unit intervals on $\mathbb R$. Run the Poisson process of intensity $\lambda=n/r$. Look at what we get in the bins. If the total number of balls $N$ is wrong, throw out or add $|N-n|$ balls at random from the whole configuration. You'll end up with the same distribution as if you threw in exactly $n$ balls in the original way. Now, before the correction procedure, the numbers of balls in bins are independent and have expectations and variances $n/r$ each. Also, since the total variance of $N$ is $n$, with high probability you'll need to use only $C\sqrt n$ balls during the correction. The claim now is that as long as the number of balls in the bin is not greater than $Cn/r$ (which has probability at least $(1-C^{-1}$), the correction can change the number of balls in the bin by only $C'r^{-1}\sqrt n$ on average, which is much less than the scaling size $\sqrt{n/r}$ if $r$ is large. Thus, the correction is invisible after scaling for each fixed finite set of bins and you can merely ignore it when passing to the limit.
