  Heavy Neutral Leptons (HNLs) are strongly motivated by theory due to their
capability of simultaneously explaining the observed phenomena of dark matter,
neutrino oscillations and the baryon asymmetry of the Universe. The existence
of such particles would affect the expansion history of the Universe and the
synthesis of primordial abundances of light elements. In this work we review,
revise and extend the phenomenology of HNLs during the Big Bang Nucleosynthesis
(BBN) epoch for masses up to 1 GeV. This is of great importance, as BBN is able
to provide complementary bounds to those from upcoming and proposed laboratory
experiments. To this end we have developed a high-precision Boltzmann code that
simulates BBN in the presence of HNLs and takes into account all relevant HNL
decay channels, as well as subsequent interactions of decay products
(thermalization and decay showers), dilution due to QCD phase transition,
active neutrino oscillations and corrections to the weak reaction rates. We
present robust bounds on the lifetime and mixing angles of HNLs for masses
$3\,\mathrm{MeV}\leq m_N \leq 1\,\mathrm{GeV}$ and show that BBN is able to
constrain HNL lifetimes down to $0.03 - 0.05$ s, depending on the mixing
pattern. Moreover, combining our results with current experimental searches, we
can exclude HNLs that mix purely with electron neutrinos up to ${\sim}$450 MeV
and those that mix purely with muon neutrinos up to ${\sim}$360 MeV, in both
cases for lifetimes up to at least a few tens of seconds. Finally, we compare
the BBN constraints with those obtained from Cosmic Microwave Background
observations and explore how our results will be improved by a number of
upcoming and proposed laboratory experiments.
