Ensemble associated with a system of Photons In the book "Fundamentals of Statistical and Thermal Physics" by F. Reif it is said that :

Photon Statistics is a special case of Bose-Einstein Statistics with no restriction on the total number of particles that the system comprises of, then the partition function in the case of Photon Statistics(Section-9.5), is written as : $$Z = \sum_{R=0} e^{-\beta(n_1e_1+n_2e_2+ ...)}$$
where, $$ n_1  e_1+n_2e_2 +... = E_R$$

And using this expression, the book goes on to derive the expression for Planck Distribution.
Now how is that possible? If I am not wrong then the expression used (in this case) for the partition function is derived for a system in the Canonical Ensemble, if the number of particles is not fixed, how can you use that expression?
 A: There are two perspectives, the first one is to consider this as an effective application of the Grand Canonical Ensemble, with chemical potential $\mu=0$. However, the specific choice the chemical potential seems arbitrarily specific and needs to be justified. The underlying for this is given in the second perspective.
The reason why you need to invoke the Canonical Ensemble is because photons are quantum excitation of modes of the EM field (see second quantization for more info). Each mode acts as an independent harmonic oscillator by construction, so you are actually applying the CE on uncoupled quantum harmonic oscillators. Since in this case the spectrum is regularly spaced out with a step $\hbar\omega$, you can interpret an energy eigenfunction as the creation of a finite number of particles at the same energy. Note that you typically set the zero point energy to $0$, which does not matter in this case since it will only add an overall factor to $Z$.
Hope this helps and tell me if you find some mistakes.
A: The formula looks like a formula in the Canonical Ensemble but is the genuine result of a Grand Canonical Ensemble calculation at the $\mu=0$ condition.
The way to see it is as follows.
The grand canonical partition function, in general, is
$$
\Xi =\sum_{N\ge 0} e^{\mu N} Q_N
$$
Where $\mu$ is the (fixed) chemical potential, and $Q_N$ is the Canonical partition function at a fixed ($N$) number of degrees of freedom.
For a separable system, the energy of a $N$-degrees of freedom state (independently whether they represent quasi-particles or else)  can be written as
$$
E(\{n_i\},N)=  \sum_i  n_i \varepsilon_i,
$$
where $\varepsilon_i$ is the eigenvalue of the $i$-th one particle state, and the values $n_i$ satisfy the constraint
$$
\sum_i n_i = N. \tag{1}
$$
Therefore, in the case $\mu=0$, we have
$$
\Xi=\sum_{N\ge 0} Q_N = \sum_{N\ge 0} \left. \sum_{\{ n_i\}} \right.^{\prime} e^{ - \beta \sum_i  n_i \varepsilon_i}.
$$
The prime on the inner summation is to recall that the configurations $\{ n_i\}$ must satisfy the constraint ($1$).
It is pretty straightforward to rewrite this last expression as
$$
\sum_{N\ge 0} \left. \sum_{\{ n_i\}} \right.^{\prime} e^{ -\beta \sum_i  n_i \varepsilon_i} = \sum_{N\ge 0} \sum_{\{ n_i\}} e^{ -\beta \sum_i  n_i \varepsilon_i} \delta_{N,\sum_i  n_i },
$$
by summing over all possible (unconstrained) values of the $\{ n_i \}$ and introducing a Kronecker's delta to take care of the constraint. In this way, we realize that by interchanging the order of summation, for each unconstrained choice of the set  $\{ n_i \}$, there is only one value of $N$ contributing to the sum. Said in another way, the sum over $N$ and the delta disappear, and we end up with
$$
\Xi =  \sum_{\{ n_i\}} e^{ -\beta \sum_i  n_i \varepsilon_i}.
$$
where each $n_i$ is independent of the value of the others and goes from $0$ to $\infty$.
