[eng] We present the IMRPhenomXHM frequency domain phenomenological waveform model for the inspiral, merger and ringdown of quasi-circular non-precessing black hole binaries. The model extends the IMRPhenomXAS waveform model [1], which describes the dominant quadrupole modes ` = |m| = 2, to the harmonics (`, |m|) = (2, 1), (3, 3), (3, 2), (4, 4), and includes mode mixing effects for the (3, 2) spherical harmonic. IMRPhenomXHM is calibrated against hybrid waveforms, which match an inspiral phase described by the effective-onebody model and post-Newtonian amplitudes for the subdominant harmonics to numerical relativity waveforms and numerical solutions to the perturbative Teukolsky equation for large mass ratios up to 1000. A computationally efficient implementation of the model is available as part of the LSC Algorithm Library Suite [2].