We study the distortions of equilibrium spectra of relic neutrinos due to the interactions with electrons, positrons, and neutrinos in the early Universe. We solve the integro-differential kinetic equations for the neutrino density matrix, including three-flavor oscillations and finite temperature corrections from QED up to the next-to-leading order $mathcal{O}(e^3)$ for the first time. In addition, the equivalent kinetic equations in the mass basis of neutrinos are directly solved, and we numerically evaluate the distortions of the neutrino spectra in the mass basis as well, which can be easily extrapolated into those for non-relativistic neutrinos in the current Universe. In both bases, we find the same value of the effective number of neutrinos, $N_{rm eff} = 3.044$, which parameterizes the total neutrino energy density. The estimated error for the value of $N_{rm eff}$ due to the numerical calculations and the choice of neutrino mixing parameters would be at most 0.0005.