Consider a single particle system with states $\psi_j$ and corresponding energy levels ${\Large\varepsilon}_j$. The partition function for this single particle system is $$Z_{(1)}=\sum_je^{-\beta\,{\Large\varepsilon}_j}$$ where $$\beta=\frac{1}{k_B\,T}$$ with $T$ as the thermodynamic temperature and $k_B$ is Boltzmann's constant.
Now consider a system of two identical weakly interacting distinguishable particles. The microstates are now $\psi=\psi_{j_1}\,\psi_{j_2}$ for every combination of the two single-particle states and the energies are $${\Large\varepsilon}_{j_1j_2}={\Large\varepsilon}_{j_1}+{\Large\varepsilon}_{j_2}\tag{1}$$
The partition function for this two-particle system is
$$\begin{align}Z_{(2)}&=\sum_{j_1}\sum_{j_2}e^{-\beta\,{\Large\varepsilon}_{j_1j_2}}\\&=\sum_{j_1}\sum_{j_2}e^{-\beta\,({\Large\varepsilon}_{j_1}+{\Large\varepsilon}_{j_2})}\\&=\sum_{j_1}\sum_{j_2}e^{-\beta\,{\Large\varepsilon}_{j_1}}e^{-\beta\,{\Large\varepsilon}_{j_2}}\\&=\left(\sum_{\color{red}{j_1}}e^{-\beta\,{\Large\varepsilon}_{\color{red}{j_1}}}\right)\left(\sum_{\color{blue}{j_2}}e^{-\beta\,{\Large\varepsilon}_{\color{blue}{j_2}}}\right)\tag{2}\\&=Z_{(1)}\times Z_{(1)}\tag{3}\\&={Z_{(1)}}^2\end{align}$$
Where in going from equation $(2)$ to $(3)$ we used the fact that $j=\color{red}{j_1}=\color{blue}{j_2}$; since it is just a dummy index of notation. Mathematically I understand completely all of the above.
But physically and intuitively I cannot understand why $Z_{(2)}={Z_{(1)}}^2$.
Looking at equation $(1)$:
$${\Large\varepsilon}_{j_1j_2}={\Large\varepsilon}_{j_1}+{\Large\varepsilon}_{j_2}$$
by my logic $Z_{(2)}={Z_{(1)}}^2$ if and only if ${\Large\varepsilon}_{j_1}={\Large\varepsilon}_{j_2}$.
But I already know that in general ${\Large\varepsilon}_{j_1}\ne{\Large\varepsilon}_{j_2}$
My same confusion applies when we factorize the partition formula for composite systems that are weakly interacting. For example, consider a single diatomic molecule in a gas. The molecule has translational, vibrational and rotational components to its motion which can be treated as independent of each other giving states with the total energy
$$\begin{align}{\Large\varepsilon}&={\Large\varepsilon}_{\mathrm{translational}}+{\Large\varepsilon}_{\mathrm{vibrational}}+{\Large\varepsilon}_{\mathrm{rotational}}\\&={\Large\varepsilon}_{\mathrm{trans}}+{\Large\varepsilon}_{\mathrm{vib}}+{\Large\varepsilon}_{\mathrm{rot}}\qquad\qquad\qquad\qquad\text{(for simplicity)}\end{align}$$
The partition function then factorizes as before and we end up with $$Z=\sum_{j_{\mathrm{trans}}}\sum_{j_{\mathrm{vib}}}\sum_{j_{\mathrm{rot}}}e^{-\beta\,\left({\Large\varepsilon}_{{j}_{\mathrm{trans}}}+{\Large\varepsilon}_{{j}_{\mathrm{vib}}}+{\Large\varepsilon}_{{j}_{\mathrm{rot}}}\right)}=Z_{\mathrm{trans}}\times Z_{\mathrm{vib}}\times Z_{\mathrm{rot}}$$
and by my logic this must mean that ${\Large\varepsilon}_{\mathrm{trans}}={\Large\varepsilon}_{\mathrm{vib}}={\Large\varepsilon}_{\mathrm{rot}}$
Clearly I am missing the point and do not fully understand the physics that is taking place here.
I would be most grateful if someone could give me some hints or an explanation that will dispel my confusion.
EDIT:
One answer addresses the fact that the index of summations will not be equal outside the summation. This part I acknowledge and understand. But since $j_1\ne j_2$ outside the summation then this must mean that $${\Large\varepsilon}_{j_1}\ne{\Large\varepsilon}_{j_2}$$ so how can we have $$Z_{(1)}\times Z_{(1)}=\sum_je^{-\beta\,{\Large\varepsilon}_j}\sum_je^{-\beta\,{\Large\varepsilon}_j}$$ when the energies are different $({\Large\varepsilon}_{j_1}\ne{\Large\varepsilon}_{j_2})$?
It's almost as if the summation has taken away information about the energies.
