Solution (source code)

= Solution

Write $\epsilon_k=E_0(k+\tfrac12)$ and use <Lagrange multipliers> $\alpha,\beta$ for the particle-number and energy constraints. The stationary equations are
$$
\frac{\partial}{\partial n_k}\left[-\sum_j(n_j\ln n_j-n_j)-\alpha\left(\sum_jn_j-N\right)-\beta\left(\sum_j\epsilon_jn_j-E\right)\right]
=-\ln n_k-\alpha-\beta\epsilon_k=0.
$$
Thus $n_k=e^{-\alpha}e^{-\beta\epsilon_k}$. Finite normalization over infinitely many levels requires $\beta>0$. Put $q=e^{-\beta E_0}$, so $0<q<1$. The <geometric series> gives
$$
N=e^{-\alpha}\frac{e^{-\beta E_0/2}}{1-q},\qquad
e^{-\alpha}=N(1-q)e^{\beta E_0/2}=2N\sinh\frac{\beta E_0}{2}.
$$
Therefore
$$
\boxed{n_k=2N\sinh\left(\frac{\beta E_0}{2}\right)e^{-\beta E_0(k+1/2)}.}
$$
For the energy constraint, use $n_k=N(1-q)q^k$ and $\sum kq^k=q/(1-q)^2$, obtained by differentiating the <geometric series>. This yields
$$
\frac E{NE_0}=(1-q)\sum_{k=0}^\infty\left(k+\frac12\right)q^k
=\frac12+\frac q{1-q}=\frac{1+q}{2(1-q)}.
$$
Recognizing the <hyperbolic cotangent> gives
$$
\boxed{E=\frac{NE_0}{2}\coth\left(\frac{\beta E_0}{2}\right).}
$$
The strict positivity of all occupations requires $E>NE_0/2$; the paper's positivity assumptions alone do not guarantee this feasibility condition. For admissible energy the multiplier is uniquely $\beta=E_0^{-1}\ln[(2E+NE_0)/(2E-NE_0)]>0$. Equality $E=NE_0/2$ is only the limiting state with all occupation in level zero, outside the stipulated strictly positive family; lower energy is impossible.

For rigor beyond formal multiplier equations, $p_k=n_k/N$ has fixed mean $E/(NE_0)-1/2$. The objective differs from $N$ times its <information entropy> only by the fixed constant $N-N\ln N$. The <maximum entropy distribution on the nonnegative integers> is uniquely geometric at that mean, confirming the stationary family and its maximum character.