$Z_{N,M}$ can be recovered from the grand partition function by $Z_{N,M}=\oint {dz \over 2\pi i} z^{-M-1} F(z)^N$, where the contour is within the radius of convergence, $|z|\leq 1$ [@0305-4470-38-19-R01]. Using Eq. (\[eq:Fz\]), one finds that the dominant contribution is made around $z_*=e^{-\mu}$ with $$\mu = e^{-\rho \zeta(2)},$$ and thus the integral is approximated in the limit $\mu\to 0$, $N\to\infty$, and $N\mu$ finite by employing the steepest descent path $z=e^{-\mu y}$ with $y$ running from $1-i\infty$ to $1+i\infty$ as $$Z_{N,M}=\mu\, \phi\left({N\mu \over \zeta(2)}\right), \ \phi(\eta)=\int_{1-i\infty}^{1+i\infty} {dy \over 2\pi i} e^{\eta (y \log y - y)}.
\label{eq:ZNM}$$ Using Eq. (\[eq:ZNM\]) into Eq. (\[eq:pm\]), we find that $$p_\infty(m) = {1\over \zeta(2) m^2} e^{{m\over N}\zeta(2)} {\phi\left({N\mu \over \zeta(2)}e^{{m\over N}\zeta(2)}\right)\over \phi\left({N\mu \over \zeta(2)}\right)}.
\label{eq:pmexact}$$

For $\phi(\eta)$ defined in Eq. (\[eq:ZNM\]), when $\eta$ is large, the contribution near $y=0$ to $\phi(\eta)$ is dominant, allowing us to approximate it as $$\phi(\eta) \simeq \int_{-\infty}^\infty {d\beta \over 2\pi} e^{-{\eta} ( 1 + {\beta^2 \over 2})} \simeq  \sqrt{1\over 2\pi \eta} e^{-{\eta}}.
\label{eq:phi1}$$ According to our numerical computation, this approximation is good even for $\eta\simeq 1$. For $\eta\ll 1$, we consider the path surrounding the branch cut of $\log y$, $y=x \pm i\epsilon$ with $x\in (-\infty,0)$ and $\epsilon\to 0$ such that 