  We present IMRPhenomXPHM, a phenomenological frequency-domain model for the
gravitational-wave signal emitted by quasi-circular precessing binary black
holes, which incorporates multipoles beyond the dominant quadrupole in the
precessing frame. The model is a precessing extension of IMRPhenomXHM
(Garc\'ia-Quir\'os 2020), based on approximate maps between aligned-spin
waveform modes in the co-precessing frame and precessing waveform modes in the
inertial frame, which is commonly referred to as "twisting up" the
non-precessing waveforms. IMRPhenomXPHM includes IMRPhenomXP as a special case,
the restriction to the dominant quadrupole contribution in the co-precessing
frame. We implement two alternative mappings, one based on a single-spin PN
approximation, as used in IMRPhenomPv2 (Hannam 2013), and one based on the
double-spin MSA approach (Chatziioannou 2017). We include a detailed discussion
of conventions used in the description of precessing binaries and of all
choices made in constructing the model. The computational cost of \phXPHM is
further reduced by extending the interpolation technique of (C.
Garc\'ia-Quir\'os 2020) to the Euler angles. The accuracy, speed, robustness
and modularity of the IMRPhenomX family will make these models productive tools
for gravitational wave astronomy in the current era of greatly increased number
and diversity of detected events.
