  The relation between the mass and spin of stellar-mass binary black holes
(BBHs) has been proposed to be a smoking gun for the presence of multiple
formation channels for compact objects. First-generation black holes (BHs)
formed by isolated binary stellar progenitors are expected to have nearly
aligned small spins, while nth-generation BBHs resulting from hierarchical
mergers are expected to have misaligned and higher spins. Leveraging data from
the third observing run O3 (GWTC-2.1 and GWTC-3), we employ hierarchical
Bayesian methods to conduct a comprehensive study of possible correlations
between the BBH masses and spins. We use parametric models that either
superpose independent BBH populations or explicitly model a mass-spin
correlation. We unveil strong evidence for a correlation between normalized
spin magnitudes and masses of BBHs. The correlation can be explained as a
transition from a BBH population with low spins at low masses and higher spins
for higher masses. Although the spin magnitude distribution at high masses
lacks robust constraints, we find strong evidence that a transition between two
BBH populations with different spin distributions should happen at 40-50
$M_{\odot}$. In particular, we find that the population of BBHs above 40-50
$M_{\odot}$ should compose the $\sim 2 \%$ of the overall population, with a
spin magnitude $\chi$ peaking around 0.7, consistently with the fraction of
nth-generation BBHs formed by hierarchical mergers in the latest
state-of-the-art BBH genesis simulations.
