  Gravitational-wave (GW) catalogs are rapidly increasing in number, allowing
for robust statistical analyses of the population of compact binaries.
Nonetheless, GW inference of cosmology has typically relied on additional
electromagnetic counterparts or galaxy catalogs. I present a new probe of
cosmological modifications of general relativity with GW data only. I focus on
deviations of the GW luminosity distance constrained with the astrophysical
population of binary black holes (BBHs). The three key observables are 1) the
number of events as a function of luminosity distance, 2) the stochastic GW
background of unresolved binaries and 3) the location of any feature in the
source mass distribution, such as the pair instability supernova (PISN) gap.
Despite a priori degeneracies between modified gravity and the unknown
evolution of the merger rate and source masses, a large damping of the GW
amplitude could be falsifiable since as redshift grows it reduces the events
and lowers the edges of the PISN gap, which is against standard astrophysical
expectations. Applying a hierarchical Bayesian analysis to the current
LIGO--Virgo catalog (GWTC-2), the strongest constraints to date are placed on
deviations from the GW luminosity distance, finding $c_{_M}=-3.2^{+3.4}_{-2.0}$
at $68\%$ C.L., which is $\sim10$ times better than multi-messenger GW170817
bounds. These modifications also affects the determination of the BBH masses,
which is crucial to accommodate the high-mass binary GW190521 away from the
PISN gap. In this analysis it is found that the maximum mass of $99\%$ of the
population shifts to lower masses with increased uncertainty,
$m_{99\%}=46.2^{+11.4}_{-9.1}M_\odot$ at $68\%$ C.L. Testing gravity at large
scales with the population of BBHs will become increasingly relevant with
future catalogs, providing an independent and self-contained test of the
standard cosmological model.
