  Hierarchical analysis of the binary black hole (BBH) detections by the
Advanced LIGO and Virgo detectors has offered an increasingly clear picture of
their mass, spin, and redshift distributions. Fully understanding the formation
and evolution of BBH mergers will require not just the characterization of
these marginal distributions, though, but the discovery of any correlations
that exist between the properties of BBHs. Here, we hierarchically analyze the
ensemble of BBHs discovered by the LIGO and Virgo with a model that allows for
intrinsic correlations between their mass ratios $q$ and effective inspiral
spins $\chi_\mathrm{eff}$. At $98.7\%$ credibility, we find that the mean of
the $\chi_\mathrm{eff}$ distribution varies as a function of $q$, such that
more unequal-mass BBHs exhibit systematically larger $\chi_\mathrm{eff}$. We
find Bayesian odds ratio of $10.5$ in favor of a model that allows for such a
correlation over one that does not. Finally, we use simulated signals to verify
that our results are robust against degeneracies in the measurements of $q$ and
$\chi_\mathrm{eff}$ for individual events. While many proposed astrophysical
formation channels predict some degree correlation between spins and mass
ratio, these predicted correlations typically act in an opposite sense to the
trend we observationally identify in the data.
