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.

In this paper we construct the first phenomenological waveform model which contains the "complete" $\ell=2$ spherical harmonic mode content for gravitational wave signals emitted by the coalescence of binary black holes with spin precession: The model contains the dominant part of the gravitational wave displacement memory, which manifests in the $(\ell=2, m=0)$ spherical harmonic in a co-precessing frame, as well as the oscillatory component of this mode. The model is constructed by twisting up the oscillatory contribution of the mode, as it was previously done for the rest of spherical harmonic modes in IMRPhenomTPHM and the Phenom family of waveform models. Regarding the displacement memory contribution present in the aligned spin (2,0) mode, we discuss a procedure to analytically compute the "precessing memory" in all the $\ell=2$ modes using the integration derived from the Bondi-Metzner-Sachs balance laws. The final waveform of the (2,0) mode is then obtained by summing together both contributions. We implement this as an extension of the computationally efficient IMRPhenomTPHM waveform model, and we test its accuracy by comparing against a set of Numerical Relativity simulations. Finally, we employ the model to perform a Bayesian parameter estimation injection analysis.