  As gravitational-wave detectors become more sensitive, we will access a
greater variety of signals emitted by compact binary systems, shedding light on
their astrophysical origin and environment. A key physical effect that can
distinguish among formation scenarios is the misalignment of the spins with the
orbital angular momentum, causing the spins and the binary's orbital plane to
precess. To accurately model such systems, it is crucial to include multipoles
beyond the dominant quadrupole. Here, we develop the first multipolar
precessing waveform model in the effective-one-body (EOB) formalism for the
inspiral, merger and ringdown (IMR) of binary black holes: SEOBNRv4PHM. In the
nonprecessing limit, the model reduces to SEOBNRv4HM, which was calibrated to
numerical-relativity (NR) simulations, and waveforms from perturbation theory.
We validate SEOBNRv4PHM by comparing it to the public catalog of 1405
precessing NR waveforms of the Simulating eXtreme Spacetimes (SXS)
collaboration, and also to new 118 precessing NR waveforms, which span mass
ratios 1-4 and spins up to 0.9. We stress that SEOBNRv4PHM is not calibrated to
NR simulations in the precessing sector. We compute the unfaithfulness against
the 1523 SXS precessing NR waveforms, and find that, for $94\%$ ($57\%$) of the
cases, the maximum value, in the total mass range $20-200 M_\odot$, is below
$3\%$ ($1\%$). Those numbers become $83\%$ ($20\%$) when using the IMR,
multipolar, precessing phenomenological model IMRPhenomPv3HM. We investigate
the impact of such unfaithfulness values with two parameter-estimation studies
on synthetic signals. We also compute the unfaithfulness between those waveform
models and identify in which part of the parameter space they differ the most.
We validate them also against the multipolar, precessing NR surrogate model
NRSur7dq4, and find that the SEOBNRv4PHM model outperforms IMRPhenomPv3HM.
