$N$ kids each brought $k$ balls to a party. When they leave each kid brings back $k$ balls randomly. Let $X$ be the total number of balls brought back by their original owners. We fix $k$. Find the distribution of $X$ when $N \rightarrow \infty$
[EDIT, removed my earlier incorrect attempt] As the comment pointed out when $k=1$ it's a rather simple Poisson distribution. But how do we prove this case for arbitrary $k$?
The distribution you want should be Poisson(k) and this should follow from Stein-Chen method. (I will be referring to https://faculty.math.illinois.edu/~psdey/414CourseNotes.pdf)
Imagine each student kept $k$ buckets in front of them, and then after the reshuffling, checked whether the bucket contained a ball of theirs. We then have $kn$ Bernoulli random variables, each with probability of success 1/n. Now, we have
$$W = \sum_{i=1}^{kn} 1_i$$
We note that most of these random variables are positively correlated, since if we return someone's ball to them, buckets belonging to others are likelier to have their owners' balls. There are some of them that are negatively correlated though, which means we cannot use any of the results in that PDF directly; however, a small modification should still get the job done.
For $Y_j^i$, we consider the coupling that is induced by at first fixing bucket $i$ to have its owner's ball, and then swapping u.a.r. (analogous to Example 8.3 on the PDF linked above.)
We thus have
\begin{align} p_i E[|U_i - V_i|] &\leq E[|X_i + \sum_{j \neq i} - Y_j^i| \\ &\leq p_i E[X_i] + \sum_{O(j) \neq O(i)} p_i E(Y_j^i - X_j) + \sum_{j \neq i, O(j) = O(i)} p_i E[X_j - Y_j^i]\\ &\leq E[X_i]^2 + \sum_{O(j) \neq O(i)} E[X_i X_j] - E[X_i]E[X_j]) \\ &- \sum_{j \neq i, O(j) = O(i)} E[X_i X_j] - E[X_i] E[X_j] \end{align}
Now, we can evaluate both sums:
$$\sum_{O(j) \neq O(i)} E[X_i X_j] - E[X_i]E[X_j]) = k(n-1) (\frac{k^2}{(nk)(nk-1)} - \frac{1}{n^2} ) = \frac{k^2 (n-1)}{n^2 (k^2 n - k)} = O(n^2) $$
$$\sum_{j \neq i, O(j) = O(i)} E[X_i X_j] - E[X_i] E[X_j] = (k-1) (\frac{k(k-1)}{nk (nk-1) } - \frac{1}{n^2}) = O(n^3)$$
and then, summing it up, we see that the overall TV distance to a Poisson of rate $\lambda = \sum_{i}^{kn} 1/n = k$ can be bounded by order O(1/n).
