Searching for additional structure and redshift evolution in the observed binary black hole population with a parametric time-dependent mass distribution

The population of the observed gravitational wave events encodes unique information on the formation and evolution of stellar-mass black holes, from the underlying astrophysical processes to the large-scale dynamics of the Universe. We use the ICAROGW analysis infrastructure to perform hierarchical Bayesian inference on the gravitational wave signals from the LIGO-Virgo-KAGRA third observing run, O3. Searching for additional structure and redshift evolution in the primary mass distribution, we explore the dependence of the mass spectrum reconstruction on different parametrizations and prior choices. For the stationary case, we find strong evidence (Bayes factor $B \simeq 180$) that the results obtained using a power-law model with a peak (Powerlaw-Gaussian)--the model preferred so far in the literature--are sensitive to prior bounds, affecting the resolvability of the $\sim 35 M_{\odot}$ peak. This behaviour is reproduced by simulated data, indicating a bimodal structure in the likelihood. Models with three mass features simultaneously capture a sharp $\sim 10M_{\odot}$ peak, a $\sim 35 M_{\odot}$ overdensity, and support for a $\sim 20 M_{\odot}$ overdensity preceded by a dip. Among these, a model with three power-law peaks (Powerlaw-Powerlaw-Powerlaw) is equally favored, in terms of evidence, over the Powerlaw-Gaussian model with wide priors. We find no statistical support for redshift evolution in the current data and provide constraints on the parameters governing this evolution, showing consistency with stationarity. We highlight possible limitations of the hierarchical Bayesian inference framework in reconstructing evolving features outside the detector horizon. Our work lays the foundations for a robust characterization of time-dependent population distributions, with significant implications for black hole astrophysics and gravitational wave cosmology.