In short, the function I am interested in is the following: $$I_n(p) = \sum_{k=0}^n \binom{n}{k} p^k (1-p)^{n-k} \left[h(p) - h\left(\tfrac{k}{n}\right)\right], \quad \text{with } h(x) \triangleq -x \log x - (1-x) \log x.$$ I am able to prove that for large $n$ this function attains a local maximum around $p \approx \frac{1.338}{n}$; assuming that $p = \frac{\alpha}{n}$ for $\alpha = O(1)$ and using a Poisson approximation of the binomial distribution, the value $\alpha \approx 1.338$ is a (numerically found) maximizing value of $$\max_{\alpha > 0} \left\{\sum_{z=1}^{\infty} \frac{\alpha^z}{z!} e^{-\alpha} (z \log z) - \alpha \log \alpha\right\}.$$ Note that the latter expression does not contain $n$ anymore; I am interested in large-$n$ asymptotics of $I_n(p)$ which leads to the above, $n$-independent expression. What I would like to show is that this is actually not just a local, but a global maximum, i.e., there are no values $p = \omega(\frac{1}{n})$ with even higher asymptotic values of $I_n(p)$ for large $n$. So my question can be formulated as:

Can we prove that for large $n$, $\max\limits_{p \in [0,1]} I_n(p) \to I_n(\frac{\alpha}{n})$ with $\alpha \approx 1.338184$?

Numerical inspection of $I_n(p)$ for say $n = 100\,000$ shows that $I_n(p)$ is nice and smooth, has two global maxima at $\frac{\alpha}{n}$ and $1 - \frac{\alpha}{n}$ and a local minimum at $p = \frac{1}{2}$ (and two trivial minima at $0$ and $1$). But clearly numerical inspection is not going to lead to a proof. I've tried various approaches using e.g. Taylor expansions of the entropy, but I am not able to prove that "order terms" are actually small for arbitrary $p$. Also, a bound given here using Jensen's inequality is not sharp enough; it's roughly a factor $2$ too loose to prove the above statement.

Any help would be appreciated!