In physics, the following problem arises in the context of statistical mechanics. Let $L \gg 1$ be a fixed parameter. Let $\beta > 0$ be also fixed. Consider the following series: $$\rho = \frac{1}{L^{3}}\sum_{p \in \frac{2\pi}{L}\mathbb{Z}^{3}}\frac{1}{e^{\beta(|p|^{2}-\mu)}-1}. \tag{1}\label{1}$$ Here, $|p|^{2} = p_{1}^{2}+p_{2}^{2}+p_{3}^{2}$, as usual and $\mu$ is (another) parameter, which is supposed to be $-\infty < \mu \le 0$. One wants to take $L\to \infty$ and replace the sum by an integral: $$\frac{1}{L^{3}}\sum_{p\in \frac{2\pi}{L}\mathbb{Z}^{3}} \to \frac{1}{(2\pi)^{3}}\int dp$$
My question is: what are the conditions which allow such replacement or, putting in another words, what conditions guarantee convergence of the sum into an integral?
Let me just explain the motivation of the question. It turns out that, if one replaces the series by an integral without worrying about the convergence, it leads to a very strange physical consequence where the density of the system is bounded when $\mu \to 0^{-}$. The explanation in the physics literature is that $\mu$ is itself a function of $L^{3}$, implicitly by (\ref{1}). Hence, one has to be careful when taking the limit $L\to \infty$. So, if the parameters $\beta$ and $\rho$ are such that, in the limit $L\to \infty$, $\mu \to 1$, then the term $(e^{-\beta\mu}-1)^{-1}$ (the term with $p=0$) in the sum (\ref{1}) diverges, and one cannot replace the sum by an integral.
So, getting back to my question, I wonder if that ensuring that every term in the sum (\ref{1}) is not divergent for every $\mu$ (also in the limit $\mu \to 0$) is sufficient to ensure that the sum converges to an integral.