  Observations of exoplanets and protoplanetary disks show that binary stellar
systems can host planets in stable orbits. Given the high binary fraction among
stars, the contribution of binary systems to Galactic habitability should be
quantified. Therefore, we have designed a suite of Monte Carlo experiments
aimed at generating large (up to $10^6$) samples of binary systems. For each
system randomly extracted we calculate the intersection between the radiative
habitable zones and the regions of dynamical stability using published
empirical formulations that account for the dynamical and radiative parameters
of both stars of the system. We also consider constraints on planetary
formation in binary systems. We find that the habitability properties of
circumstellar and circumbinary regions are quite different and complementary
with respect to the binary system parameters. Circumbinary HZs are, generally,
rare ($\simeq 4\%$) in the global population of binary systems, even if they
are common for stellar separations $\lesssim 0.2$ AU. Conversely, circumstellar
HZs are frequent ($\ge 80\%$) in the global population, but are rare for
stellar separations $\lesssim 1$ AU. These results are robust against
variations of poorly constrained binary systems parameters. We derive ranges of
stellar separations and stellar masses for which HZs in binary systems can be
wider than the HZs around single stars; the widening can be particularly strong
(up to one order of magnitude) for circumstellar regions around M-type
secondary stars. The comparison of our statistical predictions with
observational surveys shows the impact of selection effects on the habitability
properties of detected exoplanets in binary systems.
