  Self-consistent strong plasma screening around light nuclei is implemented in
the Big Bang nucleosynthesis (BBN) epoch to determine the short-range screening
potential, $e\phi(r)/T \geq 1$, relevant for thermonuclear reactions. We
numerically solve the non-linear Poisson-Boltzmann equation incorporating
Fermi-Dirac statistics adopting a generalized screening mass to find the
electric potential in the cosmic BBN electron-positron plasma for finite-sized
$^4$He nuclei as an example. Although the plasma follows Boltzmann statistics
at large distances, Fermi-Dirac statistics is necessary when work performed by
ions on electrons is comparable to their rest mass energy. While strong
screening effects are generally minor due to the high BBN temperatures, they
can enhance the fusion rates of high-$Z>2$ elements while leaving fusion rates
of lower-$Z\le 2$ elements relatively unaffected. Our results also reveal a
pronounced spatial dependence of the strong screening potential near the
nuclear surface. These findings about the electron-positron plasma's role
refine BBN theory predictions and offer broader applications for studying
weakly coupled plasmas in diverse cosmic and laboratory settings.
