  Hierarchical Bayesian inference is an essential tool for studying the
population properties of compact binaries with gravitational waves. The basic
premise is to infer the unknown prior distribution of binary black hole and/or
neutron star parameters such component masses, spin vectors, and redshift.
These distributions shed light on the fate of massive stars, how and where
binaries are assembled, and the evolution of the Universe over cosmic time.
Hierarchical analyses model the binary black hole population using a prior
distribution conditioned on hyper-parameters, which are inferred from the data.
However, a misspecified model can lead to faulty astrophysical inferences. In
this paper we answer the question: given some data, which prior
distribution--from the set of all possible prior distributions--produces the
largest possible population likelihood? This distribution (which is not a true
prior) is $\pistroke$ (pronounced "pi stroke"), and the associated
\textit{maximum population likelihood} is $\Lstroke$ (pronounced "L stroke").
The structure of $\pistroke$ is a linear superposition of delta functions, a
result which follows from Carath{\'e}odory's theorem. We show how $\pistroke$
and $\Lstroke$ can be used for model exploration/criticism. We apply this
$\Lstroke$ formalism to study the population of binary black hole mergers
observed in LIGO--Virgo--KAGRA's third Gravitational-Wave Transient Catalog.
Based on our results, we discuss possible improvements for gravitational-wave
population models.

Instrumentation and Methods for Astrophysics;General Relativity and Quantum Cosmology