The catalog of gravitational-wave events is growing, and so are our hopes of
constraining the underlying astrophysics of stellar-mass black-hole mergers by
inferring the distributions of, e.g., masses and spins. While conventional
analyses parametrize this population with simple phenomenological models, we
propose an emulation-based approach that can compare astrophysical simulations
against gravitational-wave data. We combine state-of-the-art deep-learning
techniques with hierarchical Bayesian inference and exploit our approach to
constrain the properties of repeated black-hole mergers from the
gravitational-wave events in the most recent LIGO/Virgo catalog. Deep neural
networks allow us to (i) construct a flexible single-channel population model
that accurately emulates simple parametrized numerical simulations of
hierarchical mergers, (ii) estimate selection effects, and (iii) recover the
branching ratios of repeated-merger generations. Among our results, we find the
following: The distribution of host-environment escape speeds favors values
less than $100~\mathrm{km\,s^{-1}}$ but is relatively flat, with around $37\%$
of first-generation mergers retained in their host environments;
first-generation black holes are born with a maximum mass that is compatible
with current estimates from pair-instability supernovae; there is multimodal
substructure in both the mass and spin distributions, which, in our model, can
be explained by repeated mergers; and binaries with a higher-generation
component make up at least $14\%$ of the underlying population. Though these
results are inferred through emulation of a simplified model, the deep-learning
pipeline we present is readily applicable to realistic astrophysical
simulations
