Canonical ensemble near absolute zero In the canonical ensemble for an ideal gas of $N$ bosons, the partition function for $T\to 0$ scales like 
$$Z\sim e^{-\beta\epsilon_0N},$$
when $\epsilon_0$ is the lowest (non-degenerate) single particle energy niveau, $\beta = 1/k_B T$ and the Boltzmann constant $k_B$. The energy in this limit is $E=\epsilon_0 N$ and the entropy can be obtained by  
$$S=k_B(\log Z+\beta E),$$
which is independent of $T$ and identical zero in this case.
What is the next order approximation of the partition function with respect to $T$, such that there remains a temperature dependency of the entropy? 
EDIT: 
The general partition function of the $N$-particle system for arbitrary $T>0$ is the restricted sum
$$Z(\beta,N)=\sum_{\{N=n_0+n_1+...\}}\, x_0^{n_0}\,x_1^{n_1}\,...$$
with $x_i=e^{-\beta \epsilon_i}$, for $i=0,1,2,...$ Note, that $i$ is the quantum number and $\epsilon_i<\epsilon_{i+1}$, for all $i$. 
 A: First, let's consider the case of distinguishable particles. Consider a two-level system where each of the $N$ particles is either in a state with energy $\epsilon_0$ or $\epsilon_1$. The partition function is
$$
Z = \left(\text{e}^{-\beta \epsilon_0} + \text{e}^{-\beta \epsilon_1}\right)^N
$$ 
It follows that the free energy is
$$
F = - \beta^{-1} \ln Z = - \beta^{-1} N \ln  \left(\text{e}^{-\beta \epsilon_0} + \text{e}^{-\beta \epsilon_1}\right)
$$
and the average internal energy
$$
\bar{E} = - \frac{\partial \ln Z}{\partial \beta} = N \frac{\epsilon_0\text{e}^{\beta \epsilon_1} + \epsilon_1 \text{e}^{\beta \epsilon_0}}{\text{e}^{\beta \epsilon_0} + \text{e}^{\beta \epsilon_1}}
$$
The entropy is
$$
T S = \bar{E} - F
$$
Defining $\Delta = \epsilon_1 - \epsilon_0 > 0$ to be the spacing between energy levels we find after some rearrangements that
$$
\frac{S}{k_{\text{B}}} = N \left[ \frac{\beta \Delta}{1+\text{e}^{\beta \Delta} } + \ln \left( 1 + \text{e}^{-\beta \Delta}  \right) \right]
$$
As expected the entropy depends only on the difference $\Delta$ between the energy levels and is extensive (i.e., proportional to $N$). In the limit $\beta \Delta \rightarrow \infty$ (or $T \rightarrow 0$) the entropy goes to zero as before, while when $\beta \Delta \rightarrow 0$ (or $T \rightarrow \infty$) the entropy approaches $N \ln{2}$ as expected for a completely disordered two-state system.
Second, let's consider the case of indistinguishable bosons. The partition function is now
$$
Z = \sum_{n=0}^N x_0^n \, x_1^{N-n} = \frac{x_0^{N+1} - x_1^{N+1}}{x_0-x_1}
$$
where $x_i = \text{e}^{- \beta \epsilon_i}$ is the Boltzmann factor of the $i$-th energy level. The calculation of $F$, $\bar{E}$ and $S$ works the same as before. The expressions get longer now but they are not too complicated, which makes them perfect for something like Mathematica. For the entropy I get
$$
\frac{S}{k_{\text{B}}} =\frac{\beta \Delta  \left(-e^{\beta \Delta }+N \left(-e^{\beta \Delta 
   (N+2)}\right)+(N+1) e^{\beta \Delta  (N+1)}\right)+\left(e^{\beta \Delta }-1\right)
   \left(e^{\beta \Delta  (N+1)}-1\right) \log \left(e^{\beta \Delta 
   (N+1)}-1\right)}{\left(e^{\beta \Delta }-1\right) \left(e^{\beta \Delta 
   (N+1)}-1\right)}-\log \left(e^{\beta \Delta }-1\right)
$$
For $N=1$ this expression is the same as the obtained before, as expected. In the low-temperature limit $\beta\Delta \rightarrow \infty$ the entropy goes to zero; in the high-temperature limit $\beta\Delta \rightarrow 0$ the entropy is $\ln (N+1)$.
