The gran-canonical partition function for an ideal gas of bosons can be written as $$\mathcal{Z}=\prod_{(\vec{r},\vec{p})\in\mathbb{R}^{6}}\sum_{N=0}^\infty e^{-\beta E(\|\vec{p}\|)N+\mu\beta N}=\prod_{(\vec{r},\vec{p})\in\mathbb{R}^{6}}\frac{1}{1-e^{-\beta\left(E(\|\vec{p}\|)-\mu\right)}},$$ i.e. every mode $(\vec{r},\vec{p})$ is independent and the gran-canonical partition function of each mode is obtained by summing over all possible number of particles that can contain it. Using the usual trick of taking the exponential and the logarithm to transform the product into a sum, and then replacing that sum by the usual integral, one gets $$\mathcal{Z}=\exp\left(-\int\frac{\text{d}^3\vec{r}\text{d}^3\vec{p}}{h^3}\ln\left(1-e^{-\beta\left(E(\|\vec{p}\|)-\mu\right)}\right)\right)=\exp\left(-\frac{4\pi V}{h^3}\int_0^\infty\text{d}{p}\ln\left(1-e^{-\beta\left(E(p)-\mu\right)}\right)\right).$$ In particular, for photons the correct formula is $$\mathcal{Z}=\exp\left(-\frac{8\pi V}{h^3}\int_0^\infty\text{d}{p}p^2\ln\left(1-e^{-\beta\left(cp-\mu\right)}\right)\right)=\exp\left(-\frac{8\pi V}{(\beta ch)^3}\int_0^\infty\text{d}{u}u^2\ln\left(1-e^{-u+\beta\mu}\right)\right),$$ where we have added a factor of 2 to account for the fact that every position and momenta has two distinct modes associated to it.
Now, it is usually argued, in particular with respect to studies of black body radiation, that one should further set $\mu=0$ for photons. I don't understand why one should do this. Namely, with the result above, one could compute the average number of photons to be $$N=\left(\frac{\partial\ln\mathcal{Z}}{\partial(\beta\mu)}\right)_{\beta, V}=\frac{8\pi V}{(c\beta h)^3}\int_0^\infty\text{d}{u}\frac{u^2e^{-u+\beta\mu}}{1-e^{-u+\beta\mu}}=\frac{16\pi V}{(c\beta h)^3}\text{Li}_3(e^{\beta\mu})\rightarrow\frac{16\pi \zeta(3)V}{(c\beta h)^3}$$ as $\mu\rightarrow 0$. Even if $\mu=0$, one would not be able to recover this result if one sets $\mu=0$ from the start. Is this result meaningful? If it is meaningful, why should $\mu=0$? I've heard argument that it is because the photon is massless so there is no intrinsic energy cost associated to adding a photon. Would $\mu=0$ for a gas of massless scalar particles? I've also heard that $\mu=0$ because $N$ is not conserved. Then would $\mu=0$ for any QFT? Does the fact that $N$ is not conserved make its average ill-defined?
This has been bugging me for a while, any help would be much appreciated!
