Einstein solid: one or three dimensional quantum harmonic oscillator?

The Einstein model for solids assumes all atoms vibrate with the same frequency $$\omega$$, each atom being modeled as a quantum harmonic oscillator. The thing is: solids are three-dimensional objects, hence the atoms vibrate in three dimensions.

This is addressed by doing the following:

The one-dimensional quantum harmonic oscillator has energy:

$$E_n = \left(n + \frac{1}{2} \right) \hbar \omega$$

But since the quantum harmonic oscillator is three-dimensional, we may take the atoms as oscillating independently in each of the three dimensions:

$$E_{nx} = \left(n_x + \frac{1}{2} \right) \hbar \omega$$

$$E_{ny} = \left(n_y + \frac{1}{2} \right) \hbar \omega$$

$$E_{nz} = \left(n_z + \frac{1}{2} \right) \hbar \omega$$

We then find the total energy of the solid $$E(\omega, T)$$ by using Boltzmann statistics independently in each direction, thus finding some energy $$E_x$$ for a given direction and then we find the total energy $$E = 3E_x$$, thus we find a heat capacity of $$\frac{dE}{dT} = 3 \frac{dE_x}{dT}$$.

What I'm in doubt about is that I tried to solve this same problem but by considering the atom as the three dimensional quantum harmonic oscillator of energy:

$$E_{n} = \left(n + \frac{3}{2} \right) \hbar \omega$$

And then solving for the total energy. What we find for the heat capacity using Boltzmann statistics here is the same as the heat capacity if the oscillator was just one-dimensional. This of course is wrong, since a factor of 3 is missing.

What I want to understand is why we use three "one-dimensional quantum harmonic oscillators" instead of one three-dimensional oscillator. I mean, I believe that what I have done wrong is not considering the degeneracy, but I'm not sure if that's really it. Does it have anything to do with $$n = n_x + n_y + n_z$$? If not, then what is going on?

A problem of finding a partition function of three similar noninteracting oscillators is a simple one. Due to the absence of interaction we have $$Z_3 = (Z_1)^3,$$ where $$Z_1$$ is the partition function of one oscillator $$Z_1 = \sum_{n_x=0}^\infty e^{-\frac{\hbar\omega}{\theta}(n_x+1/2)}.$$ Energy levels of a three-dimensional oscillator are degenerate indeed. The degeneracy of a level with the energy $$E_n = \hbar\omega(n+3/2)$$ is equal to $$\Gamma_n = \sum_{n_x,n_y,n_z \geq 0} \Delta(n-n_x-n_y-n_z) = \frac{(n+1)(n+2)}{2},$$ where $$\Delta(x) = 1$$ if $$x=0$$ and $$\Delta(x) = 0$$ if $$x\neq 0$$. Due to this relation we still have $$Z_3 = \sum_{n=0}^\infty \Gamma_n e^{-\frac{\hbar\omega}{\theta}(n+3/2)} = (Z_1)^3,$$ but the latter equality is not so obvious now. It holds because of the discrete delta-function in the expression for $$\Gamma_n$$.