Searching for binary black hole sub-populations in gravitational wave data using binned Gaussian processes

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_, 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_ or higher. The other sub-population with effective spins closer to zero shows a feature at 30-40M_ and is consistent with binary black holes formed dynamically in globular clusters, which are expected to peak around 30M_. 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.