Robust measurements of multiple black hole vibrational modes provide a unique opportunity to characterise gravity in extreme curvature and dynamical regimes, to better investigate the nature of compact objects and search for signs of new physics. We use a numerically-tuned quasicircular non-precessing ringdown model, TEOBPM, and the pyRing analysis infrastructure to perform a time-domain spectroscopic analysis of the third catalog of transient gravitational-wave signals, GWTC-3, searching for higher angular modes. The TEOBPM model effectively includes non-linearities in the early post-merger signal portion, and carries information about the progenitors parameters through time-dependent excitation amplitudes of the black hole quasinormal modes. Such a strategy allows us to accurately model the full post-merger emission, recovering higher signal-to-noise ratios compared to templates based on more agnostic superpositions of damped-sinusoids. We find weak evidence for the presence of (l,m)=(3,3)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$(l,m)=(3,3)$$\\end{document} [(l,m)=(2,1)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$(l,m)=(2,1)$$\\end{document}] mode in several events, with the largest Bayes factor in favour of this mode being B≃2.6\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${\\mathcal {B}} \\simeq 2.6$$\\end{document} [B≃1.2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${\\mathcal {B}} \\simeq 1.2$$\\end{document}] within the support of the peak time distribution. For GW190521, we observe B≃5.1\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$ {\\mathcal {B}} \\simeq 5.1$$\\end{document}, but only for times outside the peak time support reconstructed using the highly accurate NRSur7dq4 model, indicating significant systematics affecting such putative detection. Allowing for deviations from general relativity under the assumption of the presence of two modes, we find tentative support for the Kerr “final state conjecture”. Our work showcases a systematic methodology to robustly identify and characterise higher angular modes in ringdown-only signals, highlighting the significant impact of modelling assumptions and peak time uncertainty on spectroscopic measurements, at current signal-to-noise ratios.