Characterizing Binary Black Hole Subpopulations in GWTC-4 with Binned Gaussian Processes: On the Origins of the $35M_{\odot }$ Peak

Understanding the astrophysical origins of binary black holes requires accurate and flexible modeling of multi-dimensional population properties. In this paper, using a data-driven framework based on binned Gaussian processes, we characterize the joint distribution of BBH primary masses, mass ratios, and effective inspiral spins. We identify three distinct subpopulations in the GWTC-4 sample of observations and investigate their astrophysical origins. We find that only one of the three subpopulations exhibits the $35M_{\odot }$ peak, which is characterized by a strong preference for equal mass systems and isotropic spin orientations. Our inferred distributions are consistent with a predominantly dynamical origin of this feature. By comparing with theoretical simulations, we further show that the subpopulation that exhibits the $35M_{\text{sun}}$ peak can exclusively comprise dynamically assembled systems in globular clusters, specifically if black hole birth spins are in the range~$(0.1-0.2)$, whereas the other two subpopulations require substantial contributions from alternative formation channels. We constrain the $\mathit{\text{lower}}\mathit{\text{bound}}$ on the merger rate of BBHs in globular clusters to be $0.69_{-0.33}^{+0.23}{\mathrm{G}\mathrm{p}\mathrm{c}}^{-3}{\mathrm{y}\mathrm{r}}^{-1}$, which is consistent with theoretical predictions. We conclude that dynamical formation in globular clusters remains a strong candidate for the origin of this excess near $30-40M_{\odot }$ and that more data and targeted parametric models are necessary to rigorously establish this interpretation.