Partition function for an interacting bosonic mode Let us assume a single bosonic mode, in equilibrium with a reservoir.
For a non-interacting Bose gas, the partition function becomes
$\mathcal{Z_\text{nonint}}=\sum_{N=0}^\infty e^{-\beta(\epsilon-\mu)N}=\frac{1}{1-e^{-\beta(\epsilon-\mu)}}$.
from which it is easy to obtain a free energy $F=-T\log(\mathcal{Z})$ and particle number statistics $N=-\frac{\partial F}{\partial \mu}=\frac{1}{e^{\beta(\epsilon-\mu)}-1}$ yielding the Bose-Einstein distribution.
How is this modified for finite two-particle interaction, i.e, how to evaluate
$\mathcal{Z_\text{int}}=\sum_{N=0}^\infty e^{-\beta\left[(\epsilon-\mu)N+\frac{U}{2}N^2\right]}$ ?
 A: To evaluate the grand canonical partition function, as is done for no interaction, you need the microstates $E_i$. In the noninteracting case the energy of each microstate is simply the sum of the single particle energies, i.e. you could  write $E_i=\sum_j \epsilon_j n_{i_j}$, with $n_{i_j}$ the occupation of the single particle state $\epsilon_j$ in the microstate $i$. 
This is decomposition of many-body eigenstates into products of single particle eigenstates is no longer accurate when the interaction is nonzero. There is no way around finding all the exact eigenenergies $E_i$ of the many-body problem. Needless to say, this is extremely hard. 
Perturbation theory is therefore often used to approximate


*

*the energy of each microstate

*the spectrum itself, i.e. the types of possible excitations
The procedure is: compute a many-body ground state, do linear response on it, aka Bogolyubov, evaluate the partition function and hope for the best. 
Caveat: There are two main problems with this approach: 


*

*related to 1: the higher you go in energy, the worse your linear response approximation for the $E_i$ becomes.  

*related to 2: from exactly solvable models it is known that Bogolyubov excitations constitute only a fraction of the full excitation spectrum, see e.g. Exact Analysis of an Interacting Bose Gas. II. The Excitation Spectrum. So you do not evaluate (not even approximately) the full partition function when you rely on Bogolyubov theory. See also the plots of the excitation spectrum on Scholarpedia.
Nevertheless, if you are still interested, this question
Partition function of weakly interacting Bose gas discusses a Bogolyubov treatment. 
A: Replacing the sum by an integral, the partition function becomes
$\mathcal{Z}_{int}=\int_0^\infty e^{-\beta\left[(\epsilon-\mu)N+\frac{U}{2}N^2\right]}dN=\frac{e^{\frac{(\epsilon-\mu)^2\beta}{2U}}}{\sqrt{2U\beta/\pi}}Erfc\left[\frac{(\epsilon-\mu)\sqrt{\beta}}{\sqrt{2U}}\right]$,
at least for $U>0$. From this, a general equation of state can be derived as in the noninteracting case, with two particular limits.
In the weakly interacting limit $N^2U/T\rightarrow0$, the argument of the complementary error function goes to $+\infty$ so that an asymptotic expansion can be done, resulting to lowest order in $N\approx\frac{T}{\epsilon-\mu}\left(1-\frac{2TU}{(\epsilon-\mu)^2}\right)$ and $\frac{\sqrt{<N^2>-<N>^2}}{<N>}\approx1-\frac{U}{T}N^2$
In the strongly interacting limit $N^2U/T\rightarrow+\infty$, the argument of the complementary error fuction goes to $-\infty$ and $N\approx\frac{\mu-\epsilon}{U}$ and $\frac{\sqrt{<N^2>-<N>^2}}{<N>}\approx\sqrt{\frac{T}{UN^2}}$.
