  The number of effective relativistic neutrino species represents a
fundamental probe of the thermal history of the early Universe, and as such of
the Standard Model of Particle Physics. Traditional approaches to the neutrino
decoupling are either very technical and computationally expensive, or assume
that neutrinos decouple instantaneously. In this work, we aim to fill the gap
between these two approaches by modeling the neutrino decoupling in terms of
two simple coupled differential equations for the electromagnetic and neutrino
sector temperatures, in which all the relevant interactions are taken into
account and which allows for a straightforward implementation of BSM species.
Upon including finite temperature QED corrections we reach an accuracy on
$N_{\rm eff}$ in the SM of $0.01$. We illustrate the usefulness of this
approach to the neutrino decoupling by considering, in a model independent
manner, the impact of MeV thermal dark matter on $N_{\rm eff}$. We show that
Planck rules out electrophilic and neutrinophilic thermal dark matter particles
of $m< 3.0\,\text{MeV}$ at 95\% CL regardless of their spin, and of their
annihilation being $s$-wave or $p$-wave. We point out that thermal dark matter
particles with non-negligible interactions with both electrons and neutrinos
are more elusive to CMB observations than purely electrophilic or
neutrinophilic ones. In addition, assisted by the accuracy of our approach, we
show that CMB Stage-IV experiments will generically test thermal dark matter
particles with $m \lesssim 15\,\text{MeV}$. We make publicly available the
codes developed for this study at https://github.com/MiguelEA/nudec_BSM .
