$N$ bosons in two energy levels

For a system of $N$ indistinguishable bosons and two single-particle states (e.g. lattice sites), we have $N+1$ basis states that we can enumerate as follows: ${|N-m,mrangle}_m$, where $m=0,ldots,N$. In the basis state $|N-m,mrangle$ there are $N-m$ bosons in the single-particle state $1$ and $m$ bosons in the single-particle state $2$. The Hamiltonian $$hat{H} = E_1, hat{n}_1 + E_2, hat{n}_2 quad,$$ describes a system of non-interacting bosons. Here, the operator $hat{n}_i$ counts the number of bosons in the single-particle state $i=1,2$. Note that in particular the basis introduced above is an eigenbasis of this Hamiltonian. We define the partition function associated with $hat{H}$ as $$ Zequiv mathrm{Tr},e^{-betahat{H}} quad .$$ Using the above considerations, we obtain $$ Z= e^{-beta N E_1},sumlimits_{m=0}^{N}left(e^{-beta (E_2-E_1)}right)^m quad .tag{1}label{1}$$ We note that the $n$-th partial sum of the geometric series for $qneq1$ is given by: $$ s_n equiv sumlimits_{k=0}^n q^k= frac{1-q^{n+1}}{1-q} quad .$$

Substituting $nrightarrow N$, $qrightarrow e^{-beta (E_2-E_1)}$ and $krightarrow m$ eventually yields $$Z = e^{-beta N E_1} , left(frac{1-e^{-beta (E_2-E_1)(N+1)}}{1-e^{-beta (E_2-E_1)}}right) quad. $$

Edit: Regarding your comment: How do we obtain equation $(1)$?

First, note that the energy of a basis state $|N-m,mrangle$ is given by $$ epsilon_mequiv langle N-m,m|hat{H}|N-m,mrangle quad.$$ This yields $epsilon_m = E_1, (N-m) + E_2 ,m$. These are the eigenvalues of the Hamiltonian.

Second, the partition function is defined as the trace of the exponentiated Hamiltonian. Roughly speaking, this means that we must sum over all expectation values of this exponentiated Hamiltonian, i.e. we have to calculate:

$$Z = sumlimits_{m=0}^N langle N-m,m|e^{-betahat{H}}|N-m,mrangle quad. $$ This can be rewritten to (this should be covered in any textbook of quantum mechanics or statistical physics) $$ Z = sumlimits_{m=0}^N e^{-beta epsilon_m} quad .$$ In other words, the partition function is given by the sum of the exponentiated eigenvalues of the Hamiltonian.

As a final step, replace the expression of $epsilon_m$ with the value we have derived earlier. Then use simple rules of the exponential function to obtain equation $(1)$.

A remark: Actually, the trace of an operator is independent of the basis we choose. But the basis chosen is of course the most convenient one.