  In this work we present an extension of the time domain phenomenological
model IMRPhenomT for gravitational wave signals from binary black hole
coalescences to include subdominant harmonics, specifically the $(l=2, m=\pm
1)$, $(l=3, m=\pm 3)$, $(l=4, m=\pm 4)$ and $(l=5, m=\pm 5)$ spherical
harmonics. We also improve our model for the dominant $(l=2, m=\pm 2)$ mode and
discuss mode mixing for the $(l=3, m=\pm 2)$ mode. The model is calibrated to
numerical relativity solutions of the full Einstein equations up to mass ratio
18, and to numerical solutions of the Teukolsky equations for higher mass
ratios. This work complements the latest generation of traditional frequency
domain phenomenological models (IMRPhenomX), and provides new avenues to
develop computationally efficient models for gravitational wave signals from
generic compact binaries.
