Partition function of the single-particle vs Partition function of the system in the Canonical Ensemble In the Canonical Ensemble, given a quantum system with $N$ distingishable and non-interacting particles distributed amongst $r$ energy levels of energy $\epsilon _1,\epsilon _2,\epsilon _3,...,\epsilon _r$ and degeneracy $g_1,g_2,g_3,...,g_r$, the partition function of the single-particle is defined as
$$Z_{SP}=\sum_{i=1}^r g_ie^{-\epsilon_i\beta(T)}\tag{1}$$
with $\beta(T)=\frac{1}{K_BT}$, and the partition function of the whole system of $N$ particles is defined as
$$Z_{N}=\prod_{i=1} ^N (Z_{SP})_i =(Z_{SP})^N \tag{2}$$
(with the last equality holding if all particles are distinguishable and; in case they are identical and indistinguishable, $Z_N=Z_{SP}^N/N!$).

*

*Where does this definition, (2), come from? Why a product and not, let's say, a sum?


*On the other hand, would then be correct to define $Z_N$ also in this way (3)?
$$Z_N=\sum_{j=1}^Sg_je^{-Ej\beta(T)}\tag{3}$$
With $S$ the number of microstates of the whole system and $E_j$ the energy of the whole system at the microstate $j$.
 A: When considering a partition function of a system composed of several distinguishable subsystems you never add the separate partition functions up, and always multiply them. 
The reason is that the partition function covers the possible states of a system, and when for a system composed of subsystem we can set subsystem $A$ to a certain state and then we have to cover all of the states of subsystem $B$. Then change the state of subsystem $A$ and again sum over all states of subsystem $B$. This is multiplication
$$ Z_{AB} = \sum_{A,B} e^{-\beta(E_A+E_B)} = \sum_{A} e^{-\beta E_A} \sum_B e^{-\beta E_B} = Z_A Z_B$$
and the generalization to more than two subsystems is immediate.
Note that for this to be valid the subsystems must be separate and distinguishable. If they are interacting, for example, then you might have $E_{AB} \neq E_A + E_B$. If the particles are identical and indistinguishable, then they cannot be separated into subsystems $A$ and $B$ to begin with.
By the way - this rule of multiplication of partition functions is valid for classical systems as well as quantum ones.
A: Let us take the two-particle case explicitly: for one particle to have energy $\epsilon_1$ and the other to have energy  $\epsilon_2$, if they don't interact then the energy is $\epsilon_{(1,2)}=\epsilon_1+\epsilon_2$, further if they are distinguishable particles then the degeneracy is $g_{(1,2)}=g_1 g_2$. Both of these terms enter multiplicatively into the partition function and when you work it all out, the partition function (using your single-particle expression to sum over every pair $(i,j)$ of states) for Hamiltonian $H_A\otimes I+I\otimes H_B$ is $Z_A Z_B$ by the distributive rule of multiplication.
If the particles are indistinguishable, then this definition of the partition function, while it is still correct, turns out to be quite limited. It describes an ensemble with a fixed number of particles, but once they can be created or destroyed at will, it is useful to perform the same trick that we did with energy to get to the partition function in the first place: we connected it to a reservoir of energy and allowed energy to flow between the two systems. So when we have indistinguishable particles, the math is much easier if we connect the system to a larger system that is a reservoir of those particles, to allow particles to flow between the two systems. This is called the grand canonical ensemble, and it has a grand canonical partition function: now $i$ indexes a single particle state and it is helpful to renumber the states so that we give different numbers to each state with its different energy. Those states are occupied by $n_i$ particles and then the grand canonical partition function is
$$\mathcal Z = \sum_{i, n_i} \exp\big({-\beta (\epsilon_i(n_i) - \mu n_i)}\big).$$
Then if the particles are not interacting $\epsilon_i(n_i)=\epsilon_i(1)~n_i =\epsilon_i^\text{sp}~n_i$ and we can do the sum over $n_i$ either directly for fermions $1+\exp\left({-\beta(\epsilon_i^\text{sp}-\mu)}\right)$ or with a geometric series for bosons, $$1\over 1-\exp\left({-\beta(\epsilon_i^\text{sp}-\mu)}\right)$$
