If we have a large number $n$ of bins and throw $\lambda n$ balls into them independently with uniform probability then we expect $X_r$, the number of bins that end up with exactly $r$ balls is $E(X_r)=\frac{\lambda^r}{r!}e^{-\lambda}.$
With help from André Nicolas and Robert Israel I can compute the covariance matrix for these counts as
$\Sigma = n(D - p'p - q'q)$
where $p = p_0,p_1,\dots$, and $q = q_0,q_1,\dots$ where $p_i = \frac{\lambda^i}{i!}e^{-\lambda},$ and $q_i = \frac{i-\lambda}{\sqrt{\lambda}}\frac{\lambda^i}{i!}e^{-\lambda},$ and $D$ is the diagonal matrix made of the elements of $p$.
Is there a simple explanation for the simple form of $\Sigma$? It looks vaguely related to the covariance of a multinomial (which I'd expect) with an extra perturbation by $q'q$.