  The Lab\'ik and Smith Monte Carlo simulation technique to implement the Widom
particle insertion method is applied using Molecular Dynamics (MD) instead to
calculate numerically the insertion probability, $P_0(\eta,\sigma_0)$, of
tracer hard-sphere (HS) particles of different diameters, $\sigma_0$, in a host
HS fluid of diameter $\sigma$ and packing fraction, $\eta$, up to $0.5$. It is
shown analytically that the only polynomial representation of $-\ln
P_0(\eta,\sigma_0)$ consistent with the limits $\sigma_0\to 0$ and
$\sigma_0\to\infty$ has necessarily a cubic form,
$c_0(\eta)+c_1(\eta)\sigma_0/\sigma+c_2(\eta)(\sigma_0/\sigma)^2+c_3(\eta)(\sigma_0/\sigma)^3$.
Our MD data for $-\ln P_0(\eta,\sigma_0)$ are fitted to such a cubic polynomial
and the functions $c_0(\eta)$ and $c_1(\eta)$ are found to be statistically
indistinguishable from their exact solution forms. Similarly, $c_2(\eta)$ and
$c_3(\eta)$ agree very well with the
Boubl\'ik-Mansoori-Carnahan-Starling-Leland and
Boubl\'ik-Carnahan-Starling-Kolafa formulas. The cubic polynomial is
extrapolated (high density) or interpolated (low density) to obtain the
chemical potential of the host fluid, or $\sigma_{0}\to\sigma$, as
$\beta\mu^{\text{ex}}=c_0+c_1+c_2+c_3$. Excellent agreement between the
Carnahan-Starling and Carnahan-Starling-Kolafa theories with our MD data is
evident.
