  We calculate the evolution of the early universe through the epochs of weak
decoupling, weak freeze-out and big bang nucleosynthesis (BBN) by
simultaneously coupling a full strong, electromagnetic, and weak nuclear
reaction network with a multi-energy group Boltzmann neutrino energy transport
scheme. The modular structure of our code provides the ability to dissect the
relative contributions of each process responsible for evolving the dynamics of
the early universe in the absence of neutrino flavor oscillations. Such an
approach allows a detailed accounting of the evolution of the $\nu_e$,
$\bar\nu_e$, $\nu_\mu$, $\bar\nu_\mu$, $\nu_\tau$, $\bar\nu_\tau$ energy
distribution functions alongside and self-consistently with the nuclear
reactions and entropy/heat generation and flow between the neutrino and
photon/electron/positron/baryon plasma components. This calculation reveals
nonlinear feedback in the time evolution of neutrino distribution functions and
plasma thermodynamic conditions (e.g., electron-positron pair densities), with
implications for: the phasing between scale factor and plasma temperature; the
neutron-to-proton ratio; light-element abundance histories; and the
cosmological parameter \neff. We find that our approach of following the time
development of neutrino spectral distortions and concomitant entropy production
and extraction from the plasma results in changes in the computed value of the
BBN deuterium yield. For example, for particular implementations of quantum
corrections in plasma thermodynamics, our calculations show a $0.4\%$ increase
in deuterium. These changes are potentially significant in the context of
anticipated improvements in observational and nuclear physics uncertainties.
