  Recent population studies have searched for a subpopulation of primordial
black holes (PBHs) in the gravitational-wave (GW) events so far detected by
LIGO/Virgo/KAGRA (LVK), in most cases adopting a phenomenological PBH mass
distribution. When deriving such population from first principles in the
standard scenario, however, the equation of state of the Universe at the time
of PBH formation may strongly affect the PBH abundance and mass distribution,
which ultimately depend on the power spectrum of cosmological perturbations.
Here we improve on previous population studies on several aspects: (i) we adopt
state-of-the-art PBH formation models describing the collapse of cosmological
perturbations across the QCD epoch; (ii) we perform the first Bayesian
multi-population inference on GW data including PBHs and directly using power
spectrum parameters instead of phenomenological distributions; (iii) we
critically confront the PBH scenario with LVK phenomenological models
describing the GWTC-3 catalog both in the neutron-star and in the BH mass
ranges, also considering PBHs as subpopulation of the total events. Our results
confirm that LVK observations prevent the majority of the dark matter to be in
the form of stellar mass PBHs. We find that the best fit PBH model can comprise
a small fraction of the total events, in particular it can naturally explain
events in the mass gaps. If the lower mass-gap event GW190814 is interpreted as
a PBH binary, we predict that LVK should detect up to a few subsolar mergers
and one to $\approx 30$ lower mass gap events during the upcoming O4 and O5
runs. Finally, mapping back the best-fit power spectrum into an ultra slow-roll
inflationary scenario, we show that the latter predicts detectable PBH mergers
in the LVK band, a stochastic GW background detectable by current and future
instruments, and may include the entirety of dark matter in asteroid-mass PBHs.
