Derivation of partition function for $N$ identical quantum harmonic oscillators What is the partition function
$$\mathcal Z^{(N)}_\beta(H) : =\mathrm{Tr}\exp(-\beta H) \tag{Z} $$
$\left(\beta >0\right)$ for a system of $N$ indistinguishable and non-interacting bosons (e.g. harmonic oscillators) with Hamiltonian $H$?
Why is $$\mathscr Z^{(N)}_\beta(H) := \frac{\mathcal Z_\beta^{(1)}(h)^N}{N!}\tag{Z$^\prime$}$$
with the corresponding single particle Hamiltonian $h$ and the partition function for a single particle $\mathcal Z_\beta^{(1)}$ not the correct partition function? Is it an approximation? If so, under which (physical) circumstances is this a good approximation? How does the factor $N!$ arise?
 A: The problem with the partition function in $(\mathrm Z^\prime)$ is that there the physical states are not counted correctly (cf. the answer by @SolubleFish). However, this partition function can be obtained in a certain limit, as shown in the end. To obtain the correct expression for the partition function, we should start more or less from scratch:
The Hilbert space of $N$ identical bosons is given by $\mathcal H := \vee^N \mathcal H_1$, where $\mathcal H_1$ is the single-particle Hilbert space. If $h$ denotes the Hamiltonian for a single particle (e.g. harmonic oscillator), then the Hamiltonian for the system of interest is given by
$$H:= \sum\limits_{i=1}^N h_i \quad ,  \tag{1}$$
where $h_i:=  \mathbb I \otimes \ldots \otimes h \otimes \ldots \otimes \mathbb I$ ($h$ is at the $i$-th position and the total number of factors is $N$).
Let $\{|k\rangle\}_{k \in \mathbb N_0} \subset \mathcal H_1$ denote the eigenbasis of $h$ with $h |k\rangle = \epsilon_k |k\rangle$ and
$$|k_1,k_2,\ldots , k_N\rangle := \sqrt{\frac{N!}{n_{k_1}!n_{k_2}!\cdots n_{k_N}!}} S \left(|{k_1}\rangle \otimes |{k_2}\rangle \otimes \ldots \otimes |{k_N}\rangle \right) \quad , \tag{2} $$
with the symmetrization operator $S:=\frac{1}{N!} \sum\limits_p P$ and the permutation operator for $N$ particles $P$. Then a basis in $\mathcal H$ is given by (we employ some standard ordering):
$$\{|k_1,k_2,\ldots ,k_N\rangle\}_{k_1\leq k_2 \leq \ldots \leq k_N} \quad . \tag{3}$$
It holds that
$$\mathbb I_{\mathcal H} = \sum\limits_{k_1=0}^{\infty}\sum\limits_{k_2=k_1}^{\infty}  \cdots \sum\limits_{k_N=k_{N-1}}^\infty |k_1,k_2,\ldots, k_N\rangle \langle k_1,k_2,\ldots ,k_N|  \tag{4}$$
and $$\langle k_1,k_2,\ldots ,k_N|k_1^\prime,k_2^\prime,\ldots k^\prime_N\rangle = \delta_{k_1,k^\prime_1}\, \delta_{k_2,k^\prime_2} \cdots\, \delta_{k_N,k^\prime_N} \quad . \tag{5}$$
Moreover, we have $$H|k_1,k_2,\ldots, k_N\rangle = \left(\epsilon_{k_1} +\epsilon_{k_2} + \ldots  +\epsilon_{k_N}  \right)|k_1,k_2,\ldots, k_N\rangle \quad . \tag{6}$$
The partition function in equation $(\mathrm{Z})$ can now be simplified. Indeed, by making use of $(4)$-$(6)$ we find
\begin{align}
\mathcal Z^{(N)}_\beta (H) &= \sum\limits_{k_1=0}^{\infty}\sum\limits_{k_2=k_1}^{\infty}  \cdots \sum\limits_{k_N=k_{N-1}}^\infty  \exp-\beta\,\left(\epsilon_{k_1} +\epsilon_{k_2} + \ldots  +\epsilon_{k_N}\right)\\
 &= \sum\limits_{k_1=0}^{\infty} e^{-\beta\epsilon_{k_1}}\sum\limits_{k_2=k_1}^{\infty} e^{-\beta\epsilon_{k_2}}\, \cdots \sum\limits_{k_N=k_{N-1}}^\infty e^{-\beta\epsilon_{k_N}} \quad . \tag{7}
\end{align}

For a quantum harmonic oscillator, we have $\epsilon_k := \hbar \omega \left(k +\frac{1}{2}\right)$ for some $\omega >0$. By defining $q:= e^{-\beta \hbar \omega}$ we can hence write
$$\mathcal Z^{(N)}_\beta(H) = q^{N/2} \sum\limits_{k_1=0}^{\infty} q^{k_1}\sum\limits_{k_2=k_1}^{\infty} q^{k_2} \cdots \sum\limits_{k_N=k_{N-1}}^\infty q^{k_N}\quad . \tag{8}$$
To proceed, note that here $0<|q|<1$ and thus
$$\sum\limits_{k=m}^\infty q^k = \frac{q^m}{1-q} \tag{9} \quad . $$
Consequently, we find
\begin{align}
\mathcal Z^{(N)}_\beta (H) &= q^{N/2} \sum\limits_{k_1=0}^{\infty} q^{k_1}\sum\limits_{k_2=k_1}^{\infty} q^{k_2} \cdots \sum\limits_{k_{N-1}=k_{N-2}}^\infty q^{k_{N-1}}\sum\limits_{k_N=k_{N-1}}^\infty q^{k_N}\\
  &= q^{N/2} \sum\limits_{k_1=0}^{\infty} q^{k_1}\sum\limits_{k_2=k_1}^{\infty} q^{k_2} \cdots \sum\limits_{k_{N-1}=k_{N-2}}^\infty q^{k_{N-1}} \frac{q^{k_{N-1}}}{1-q} \tag{10}\\
 &= \frac{q^{N/2}}{1-q} \sum\limits_{k_1=0}^{\infty} q^{k_1}\sum\limits_{k_2=k_1}^{\infty} q^{k_2} \cdots \sum\limits_{k_{N-1}=k_{N-2}}^\infty (q^2)^{k_{N-1}} \\&= \ldots 
\end{align}
which eventually yields
$$\mathcal Z^{(N)}_\beta(H) = q^{N/2} \prod\limits_{i=1}^N \frac{1}{1-q^i} \quad . \tag{11} $$
The reader may formalize this argument as an exercise.
The result coincides with equation $(14)$ of Investigations on finite ideal quantum gases and $(9)$ of Statistical mechanics and the partitions of numbers.

The second linked paper shows/ states that
$$\mathcal Z^{(N)}_\beta(H) \underbrace{\longrightarrow}_{N^2 \mu\to 0} \mathscr Z_\beta^{(N)}(H) \quad , \tag{12}  $$
with $\mu:=\beta\hbar\omega$. Alternatively,  we can use the following relation, which may be proven by induction:
$$ N!= \lim\limits_{\mu \to 0} \frac{\prod\limits_{i=1}^N 1-e^{-\mu i}}{(1-e^{-\mu})^N}  = \lim\limits_{\mu \to 0} \frac{\mathcal Z_\beta^{(1)}(h)^N}{\mathcal Z^{(N)}_\beta(H)}\quad. \tag{13} $$
This shows that for a sufficiently small $\mu^*(N)\gtrsim 0$ it holds that
$$\mathscr Z_{\beta^*}^{(N)}(H) \approx  \mathcal Z_{\beta^*}^{(1)}(h)^N \, \frac{\mathcal Z^{(N)}_{\beta^*}(H) }{ \mathcal Z_{\beta^*}^{(1)}(h)^N}  = \mathcal Z^{(N)}_{\beta^*}(H) \quad .  \tag{14} $$
The reader is encouraged to test this approximation for some finite $\mu$ as a function of $N$. We find that the larger $N$ is, the smaller $\beta$, i.e. the higher the temperature, must be, in order to be a good approximation. This is in a qualitatively agreement with the limit taken in $(12)$.
All in all, in the high temperature limit, the partition function $\mathscr Z$ should be a good approximation for $\mathcal Z$, at least in our example.
A: I would like to add something to @JasonFunderberker's answer :
$Z'(h,N,\beta) = \mathcal Z^{(1)}_\beta(h)^N$ is the partition function for $N$ distinguishable particles. The Hilbert space for this system is just the $n$-fold tensor product of the $1$-particle states.
When, the particle are indistinguishable, we have to make sure we are not over-counting states. Dividing the whole partition function by $N!$ would be correct if every state was overcounted $N!$ times in $Z'$ above, or equivalently in the Hilbert space of symmetric states.
But this is not correct : while a generic states is indeed overcounted by a factor of $N!$, some states are not. For example, if all the particles are in the same state, there is no overcounting at all.
