  Astrophysically motivated population models for binary black hole observables
are often insufficient to capture the imprints of multiple formation channels.
This is mainly due to the strongly parametrized nature of such investigations.
Using a non-parametric model for the joint population-level distributions of
binary black hole component masses and effective inspiral spins, we find hints
of multiple subpopulations in the third gravitational wave transient catalog.
The higher (more positive) spin subpopulation is found to have a mass-spectrum
without any feature at the $30-40M_{\odot}$, which is consistent with the
predictions of isolated stellar binary evolution, simulations for which place
the pile up due to pulsational pair-instability supernovae near $50M_{\odot}$
or higher. The other sub-population with effective spins closer to zero shows a
feature at $30-40M_{\odot}$ and is consistent with binary black holes formed
dynamically in globular clusters, which are expected to peak around
$30M_{\odot}$. We also compute merger rates for these two subpopulations and
find that they are consistent with the theoretical predictions of the
corresponding formation channels. We validate our results by checking their
robustness against variations of several model configurations and by analyzing
large simulated catalogs with the same model.
