Balls and Bins analysis in the $m \gg n$ regime: gaps Suppose we are throwing $m$ balls into $n$ bins, where $m \gg n$. Let $X_i$ be the number of balls ending up in bin $i$, $X_\max$ be the heaviest bin, $X_\min$ be the lightest bin, and $X_{\mathrm{sec-max}}$ be the second heaviest bin. Roughly speaking, $X_i - X_j \sim N(0,2m/n)$, and so we expect $|X_i - X_j| = \Theta(\sqrt{m/n})$ for any two fixed $i,j$. Using a union bound, we expect $X_{\max} - X_{\min} = O(\sqrt{m\log n/n})$; presumably, we can get a matching lower bound by considering $n/2$ pairs of disjoint bins. This (not completely formal) argument leads us to expect that the gap between $X_{\max}$ and $X_{\min}$ is $\Theta(\sqrt{m\log n/n})$ with high probability.

I am interested in the gap between $X_\max$ and $X_{\mathrm{sec-max}}$. The argument outlined above shows that $X_\max - X_{\mathrm{sec-max}} = O(\sqrt{m\log n/n})$ with high probability, but the $\sqrt{\log n}$ factor seems extraneous. Is anything known about the distribution of $X_\max - X_{\mathrm{sec-max}}$?

More generally, suppose that each ball is associated with a non-negative score for each bin, and we are interested in the total score of each bin after throwing $m$ balls. The usual scenario corresponds to scores of the form $(0,\ldots,0,1,0,\ldots,0)$. Suppose that the probability distribution of the scores is invariant under permutation of the bins (in the usual scenario, this corresponds to the fact that all bins are equiprobable). Given the distribution of the scores, we can use the method of the first paragraph to get a good bound on $X_{\max} - X_{\min}$. The bound will contain a factor of $\sqrt{\log n}$ that comes from a union bound (via the tail probabilities of a normal variable). Can this factor be reduced if we're interested in bounding $X_{\max} - X_{\mathrm{sec-max}}$?
Answer: $\Theta\left(\sqrt{\frac{m}{n\log n}}\right)$.

Applying a multidimensional version of the Central Limit Theorem, we get that the vector $(X_1,\dots, X_n)$ has asymptotically multivariate Gaussian distribution with 
$$\mathrm{Var}[X_i] = m\left(\frac{1}{n} - \frac{1}{n^2}\right),$$
and
$$\mathrm{Cov}(X_i, X_j) = -m/n^2.$$
We will assume below that $X$ is a Gaussian vector (and not only approximately a Gaussian vector). Let us add a Gaussian random variable $Z$ with variance $m/n^2$ to all $X_i$ ($Z$ is independent from all $X_i$). That is, let
$$
\begin{pmatrix}
Y_1\\Y_2\\ \vdots\\Y_n 
\end{pmatrix}
=
\begin{pmatrix}
X_1+Z\\X_2+Z\\ \vdots\\X_n +Z
\end{pmatrix}.
$$
We get a Gaussian vector $(Y_1, \dots, Y_n)$. Now each $Y_i$ has variance $m/n$:
$$\mathrm{Var}[Y_i] = \mathrm{Var}[X_i] + \underbrace{2\mathrm{Cov}(X_i,Z)}_{=\, 0}+\mathrm{Var}[Z] = m/n,$$
and all $Y_i$ are independent:
$$\mathrm{Cov}(Y_i, Y_j) = \mathrm{Cov}(X_i, X_j) + \underbrace{\mathrm{Cov}(X_i,Z) + \mathrm{Cov}(X_j,Z)}_{=\, 0} +\mathrm{Cov}(Z, Z) = 0.$$

Note that $Y_i - Y_j = X_i - X_j$. Thus our original problem is equivalent to the problem of finding $Y_{\mathrm{max}} - Y_{\mathrm{sec-max}}$. Let us first for simplicity analyze the case when all $Y_i$ have variance $1$.


  Problem. We are given $n$ independent Gaussian r.v. $\gamma_1,\dots, \gamma_n$ with mean $\mu$ and variance $1$. Estimate the expectation of $\gamma_{\mathrm{max}} - \gamma_{\mathrm{sec-max}}$.


Answer: $\Theta\left(\frac{1}{\sqrt{\log n}}\right)$.

Informal Proof. 
Here is an informal solution to this problem (it's not hard to make it formal).
Since the answer does not depend on the mean, we assume that $\mu = 0$. Let $\bar\Phi(t) = \Pr[\gamma &gt; t]$, where $\gamma\sim{\cal N}(0,1)$. We have (for moderately large $t$),
$$\bar\Phi(t)\approx \frac{1}{\sqrt{2\pi}t} e^{-\frac{1}{2}t^2}.$$

Note that


$\Phi(\gamma_i)$ are uniformly and independently distributed on $[0,1]$, 
$\Phi(\gamma_{\mathrm{max}})$ is the smallest among $\Phi(\gamma_i)$,
$\Phi(\gamma_{\mathrm{sec-max})}$ is the second smallest among $\Phi(\gamma_i)$.


Thus $\Phi(\gamma_{\mathrm{max}})$ is close to $1/n$ and $\Phi(\gamma_{\mathrm{max}})$ is close to $2/n$ (there is no concentration but if we don't care about constants these estimates are good enough; in fact, they are even pretty good if we care about constants — but that needs a justification). Using the formula for $\bar\Phi(t)$, we get that
$$
2\approx \bar\Phi(\gamma_{\mathrm{sec-max}})\left/\bar\Phi(\gamma_{\mathrm{max}})\right. \approx e^{\frac{1}{2}\left(\gamma_{\mathrm{max}}^2 - \gamma_{\mathrm{sec-max}}^2\right)}.
$$

Thus $\gamma_{\mathrm{max}}^2 - \gamma_{\mathrm{sec-max}}^2$ is $\Theta(1)$ w.h.p. Note that $\gamma_{\mathrm{max}}\approx \gamma_{\mathrm{sec-max}} = \Theta(\sqrt{\log n})$. We have, 
$$\gamma_{\mathrm{max}} - \gamma_{\mathrm{sec-max}}\approx \frac{\Theta(1)}{\gamma_{\mathrm{max}} + \gamma_{\mathrm{sec-max}}} \approx \frac{\Theta(1)}{\sqrt{\log n}}.$$ 

QED

We get that 
\begin{align}
\mathbb{E}[{X_{\mathrm{max}} - X_{\mathrm{sec-max}}}] &amp;= \mathbb{E}[{Y_{\mathrm{max}} - Y_{\mathrm{sec-max}}}] \\
&amp;= \sqrt{\mathrm{Var}[Y_i]} \times\mathbb{E}[{\gamma_{\mathrm{max}} - \gamma_{\mathrm{sec-max}}}] = \Theta\left(\sqrt{\frac{m}{n\log n}}\right).
\end{align}


  The same argument goes through when we have arbitrary scores. It shows that 
  $$\mathbb{E}[X_{\mathrm{max}}- X_{\mathrm{sec-max}}] = c\, \left. \mathbb{E}[X_{\mathrm{max}}- X_{\mathrm{min}}]\right/\log n.$$