In this paper we present an extensive analysis of the GW190521 gravitational wave event with the current (fourth) generation of phenomenological waveform models for binary black hole coalescences. GW190521 stands out from other events since only a few wave cycles are observable. This leads to a number of challenges, one being that such short signals are prone to not resolve approximate waveform degeneracies, which may result in multi-modal posterior distributions. The family of waveform models we use includes a new fast time-domain model IMRPhenomTPHM, which allows us extensive tests of different priors and robustness with respect to variations in the waveform model, including the content of spherical harmonic modes. We clarify some issues raised in a recent paper [Nitz&Capano], associated with possible support for a high-mass ratio source, but confirm their finding of a multi-modal posterior distribution, albeit with important differences in the statistical significance of the peaks. In particular, we find that the support for both masses being outside the PISN mass-gap, and the support for an intermediate mass ratio binary are drastically reduced with respect to what Nitz&Capano found. We also provide updated probabilities for associating GW190521 to the potential electromagnetic counterpart from ZTF.