Suppose that I have a random variable $X \sim \text{Multinomial}(N, M, \mathbf{p})$, where $\mathbf{p} = \frac{1}{M} \mathbf{1}$, $N$ is the number of trials, and $M$ is the number of bins.
I want to find out:
$$ P\left(\text{max}\left\{ x_{2}, x_{3}, \ldots, x_{M} \right\} \leq \gamma \right) $$
In other words, I want to find out the probability that the maximum of all but one of the $x_{i}$ are less than a threshold.
My immediate thought is that because $x_{i} \sim \text{Binomial}\left( N, p \right)$, that the answer would simply be:
$$ P\left(\text{max}\left\{ x_{2}, x_{3}, \ldots, x_{M} \right\} \leq \gamma \right) = \prod_{i=2}^{M} P(x_{i} \leq \gamma) = P(x_{i} \leq \gamma)^{M-1} $$
Therefore, if I wanted the maximum of all but $k$ of the $x_{i}$, I would have:
$$ P\left(\text{max}\left\{ x_{k+1}, x_{k+2}, \ldots, x_{M} \right\} \leq \gamma \right) = \prod_{i=k+1}^{M} P(x_{i} \leq \gamma) = P(x_{i} \leq \gamma)^{M-k} $$
However, this doesn't seem to be the case when I compare against a numerical experiment. The figure shows $P(\text{max}\left\{ x_{2}, x_{3} \right\} \leq \gamma)$ when $X \sim \text{Multinomial}\left(N=1000, M=3, p = [\frac{1}{3}, \frac{1}{3}, \frac{1}{3}]^{T} \right)$. The red curve represents my thoughts above and the blue curve represents the hard data.
So, what am I missing?
What you're missing is that the $x_i$ are not independent binomials. The sum $x_1 + x_2 + \dots + x_n$ is fixed to $N$, so if one of them is small, that makes the rest likely to be larger.
Consider the $\gamma=0$ case of your example, chosen because it's the easiest to compute. We have $\max\{x_2, x_3\} \le \gamma$ precisely when $X = (N,0,0)$, which happens with probability $(\frac13)^N$. But the individual probability that $x_i \le \gamma$ is $(\frac23)^N$, so your method gives $(\frac23)^{2N} = (\frac49)^N$ as the answer instead.
I don't expect a closed form to exist for your expression, because even $\Pr[x_i \le \gamma]$ has no closed form.