A combinatorics problem related to Bose-Einstein statistics

Problem

Given $N$ and $E$, how many solutions $(n_0,n_1,...,n_E)$ are there to: $$\sum_{k=0}^En_k = N$$ $$\sum_{k=0}^Ek.n_k = E$$ where everything is non-negative integer ($N,E,n_k \in \mathbb{Z}_{\geq 0}$ )?

Context

$N$ is the number of indistinguishable particles that can occupy $Q+1$ distinguishable energy levels, with energies $E_k=0,1,2,...,Q$. In this case, $Q=E$, which is the total energy of the system.

Here's an example for $N=6$ and $E=9$:

Each rectangle represents a solution. There are 26 in total. The solution on the top left, for example, is $(5,0,0,0,0,0,0,0,0,1)$.

Source: HyperPhysics.

Attempts

The equivalent problem for distinguishable particles is trivial if solved by the stars and bars method. The solution is given by:

$$\binom{E+N-1}{E}$$

For $N=6$ and $E=9$ it would be:

$$\binom{9+6-1}{9} = \frac{14!}{9!\;5!} = 2002$$

I don't know if this is relevant, but I noticed that $2002 = 77\times 26$, or $77$ times the number of solutions to the problem I want to solve.

There might not even be a closed formula. Many simple combinatorics problems don't seem to have a simple solution.

Maybe the problem can be solved by generating functions, as suggested by the answer to this related question, but I have no idea how that works.

You're right about there not being a closed form, but there being a generating function solution.

In combinatorics, the objects you are counting are called partitions, which are ways of writing one number as a sum of positive integers.

Interpret $n_k$ as the number of times the integer $k$ appears in the sum. Then $$\sum_{k =1}^E k \cdot n_k = E$$ means that the numbers sum up to $E$, so we are interested in partitions of $E$ (note that we can ignore $k=0$ here as it does not affect this sum).

The first condition, $\sum_{k=0}^E n_k = N$, or equivalently $\sum_{k=1}^E n_k \le N$, means that we are interested in partitions of $E$ with at most $N$ parts.

Now unfortunately there is no closed solution to this problem. There is, however, a generating function. We first make a transformation.

We claim that the number of partitions of $E$ into at most $N$ parts is the same as the number of partitions of $E$ into parts of size at most $N$. To see this, let $m_j = \sum_{k=j}^E n_k$ count the number of parts of size at least $j$. Clearly $m_j=0$ for $j>N$, and $$\sum_j m_j = \sum_j \sum_{k \ge j} n_k = \sum_k \sum_{j \le k} n_k = \sum_k k \cdot n_k,$$ so given our partition $(k \cdot n_k)_k$ [that is, $n_k$ copies of $k$ as $k$ ranges from $1$ to $E$] of $E$ with at most $N$ parts, we formed a partition $(m_j)_j$ of $E$ with parts of size at most $N$. One can check that this is in fact a bijection, so we can count either type of partition.

Now for the generating function. A partition represents the (unordered) sum $m_1 + m_2 + ... + m_n = E$, where $n$ is the largest integer for which $m_j \ge 1$. One convenient way of representing such sums is through polynomial multiplication: $x^{m_1} \cdot x^{m_2} \cdot ... \cdot x^{m_n} = x^{E}$.

When multiplying out these monomials, we can group them together by their powers. So the terms with $m_j =1$ can be collected into the factor $x^{1 \cdot 0} + x^{1 \cdot 1} + x^{1 \cdot 2} + ... + x^{1 \cdot s_1} + ...$, where in $1 \cdot s_1$, the $1$ indicates we are dealing with parts of size $1$, and the $s_1$ keeps track of how many there are. The cool thing is that this is a geometric series, so we have $$x^{1 \cdot 0} + x^{1 \cdot 1} + x^{1 \cdot 2} + ... + x^{1 \cdot s_1} + ... = \frac {1}{1-x^1}.$$

We can do this for every other part size, and the factor corresponding to the parts of size $j$ will be $\frac {1}{1-x^j}$.

Finally, to get the generating function, we multiply these factors together, because remember that adding the parts corresponds to multiplying the monomials. Since we only want parts of size up to $N$, we take the product over $1 \le j \le N$, giving $$F_N (x) = \prod_{j=1}^N (1-x^j)^{-1}.$$

This gives a power series, and the coefficient of $x^E$ counts the partitions of $E$ into parts of size at most $N$ (or, equivalently, with at most $N$ parts), which is what you are looking for.

Even more finally, what good is a generating function? Well, it gives a systematic way to approach these problems. Note that the function $F_N$ doesn't just encode the solution to your problem, but to the problem for every possible value of $E$, all at once! In practical terms, it gives an algorithm for finding the solution, since the power series can be computed efficiently. We humans can also use analytic tools to extract the asymptotics of the solution, although that is not something I can do off the top of my head, so I'll leave that part to someone else.