Computationally efficient models for the dominant and sub-dominant harmonic modes of precessing binary black holes.

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.