Binomial probabilities at a fixed distance from the mode: size-biasing and the complete asymptotic expansion

Let $X\sim Bin(N,p)$ with $0<p<1$ and $q=1-p$, let $\nu=\lceil Np\rceil$ be the first lattice point not below the mean, and let $r$ be a fixed integer.

We determine the complete asymptotic expansion, in powers of $N^{-1}$, of the binomial mass $\Pr\{X=\nu+r\}$ -- a quotient of gamma functions -- uniformly for $p$ in compact subintervals of $(0,1)$ and for bounded $r$.

Because the mean $Np$ is not a lattice point, the coefficients cannot be constants: they are Bernoulli polynomials evaluated at the oscillating fractional displacement $h_N=\nu-Np\in[0,1)$, and are given in closed form to all orders.

Three structural features organise the answer.

The elementary tail of the naive Stirling expansion sums exactly to $-\log\bigl((\nu+r)/(Np)\bigr)$ and is removed by the binomial size-bias identity, leaving a pure Appell series in the Bernoulli polynomials.

Restoring the symmetry $p\leftrightarrow q$ replaces these by the even Appell sequence generated by $\tfrac z2\coth\tfrac z2$, the reciprocal Bernoulli polynomials, and returns the classical Stirling prefactor.

The integer shift $r$ contributes only Euler--Maclaurin power sums, the summation rule being fixed by the Appell sequence in force.

As consequences, the coefficient of $N^{-1}$ already yields the exact rule for the mode, with the skewness offset $\tfrac12(p-q)$; the case $p=\tfrac12$ recovers the central binomial coefficient; and the Cesàro means of the oscillating coefficients are obtained in closed form.