  Primordial black hole (PBH) formation during cosmic phase transitions and
annihilation periods, such as the QCD transition or the $e^+e^-$-annihilation,
is thought to be particularly efficient due to a softening of the equation of
state. We present a detailed numerical study of PBH formation during the QCD
epoch in order to derive an accurate PBH mass function. We also briefly
consider PBH formation during the $e^+e^-$-annihilation epoch. Our
investigation confirms that, for nearly scale-invariant spectra, PBH abundances
on the QCD scale are enhanced by a factor $\sim 10^3$ compared to a purely
radiation dominated Universe. For a power spectrum producing an (almost)
scale-invariant PBH mass function outside of the transition, we find a peak
mass of $M_{\rm pbh}\approx 1.9 M_{\odot}$ with a fraction $f\approx 1.5\times
10^{-2}$ of the PBHs having a mass of $M_{\rm pbh} > 10 M_{\odot}$, possibly
contributing to the LIGO-Virgo black hole merger detections. We point out that
the physics of PBH formation during the $e^+e^-$-annihilation epoch is more
complex as it is very close to the epoch of neutrino decoupling. We argue that
neutrinos free-streaming out of overdense regions may actually hinder PBH
formation.
