I was solving a physics exercise and I encountered this formula: $$\left< n_l \right>=\left[1+\sum_{k\neq l} \left(e^{bN(l-k)}\frac{\prod_{j\neq l} (1-e^{b(l-j)})}{\prod_{j\neq k} (1-e^{b(k-j)})}\right)\right]^{-1} \left( N+\sum_{h\neq l} \frac{1}{1-e^{b(h-l)}}\right) - \\\;\;\;\sum_{h\neq l} \left\{ \left[1+\sum_{k\neq h} \left(e^{bN(h-k)}\frac{\prod_{j\neq h} (1-e^{b(h-j)})}{\prod_{j\neq k} (1-e^{b(k-j)})}\right)\right]^{-1} \left( \frac{1}{1-e^{b(l-h)}} \right) \right\}$$
where $b$ is a positive real number, $N$ is a natural number, and all the sums and products are understood to run from $0$ to $M-1$, where $M$ is yet another natural number (all the $l$, $h$, $j$, $k$ indices are therefore bound to this interval, with the appropriate restrictions written under the sums and the products). Eventually I would also like to take the limit as $M$ goes to infinity, but I'm sure how well I can do that.
Do you think there is a hope to massage this formula to make it a little less ugly? For the moment, I've been trying to solve it numerically...
Additional info:
This formula arises in the context of an exercise I invented myself, and I haven't found any reference in the literature.
The problem is to find the mean occupation number of $M$ energy levels $\left\{ \varepsilon_l \mid l=0,... M-1\right\}$, populated by a Bose gas of $N$ indistinguishable quantum particles at temperature $T$. The system I consider is closed, so I used methods from statistical mechanics and calculated the canonical partition function (this approach is different from the one used classically in literature to derive the Bose-Einstein statistics, because the latter works in the gran-canonical ensemble, where the number of particles is not held fixed).
I so found, after some calculation, the canonical partition function to be
$$Z(N)=\sum_{l} \frac{e^{-\beta N \varepsilon_l}}{\prod_{k\neq l} \left( 1- e^{\beta (\varepsilon_l - \varepsilon_k)}\right)}$$
where $\beta=\frac{1}{k_BT}$.
One can then demonstrate that the mean occupation number $\left< n_l \right>$ of the energy level $\epsilon_l$ is given by
$$\left< n_l \right>=-\frac{1}{\beta} \frac{\partial}{\partial \varepsilon_l} \ln(Z(N))$$
So I evaluated this derivative and I considered it in a special case of this problem, where all the levels are equally spaced: $\varepsilon_l \triangleq l\varepsilon$, where $\varepsilon$ is a positive energy scale.
The result I found for $\left< n_l \right>$ is precisely the first formula of this post, where I renamed the product $\beta\varepsilon$ as $b$.
