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_{\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 \textit{lower bound} on the merger rate of BBHs in globular clusters to be $0.69^{+0.23}_{-0.33} \rm{Gpc}^{-3}\rm{yr}^{-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.