Quasi-normal mode (QNM) modeling is an invaluable tool for studying strong gravity, characterizing remnant black holes, and testing general relativity. To date, most studies have focused on the dominant (2,2) mode, and have fit to standard strain waveforms from numerical relativity. But, as gravitational wave observatories become more sensitive, they can resolve higher-order modes. Multimode fits will be critically important, and in turn require a more thorough treatment of the asymptotic frame at ℐ⁺. The first main result of this work is a method for systematically fitting a QNM model containing many modes to a numerical waveform produced using Cauchy-characteristic extraction, which is known to exhibit memory effects. We choose the modes to model based on their power contribution to the residual between numerical and model waveforms. We show that the all-angles mismatch improves by a factor of ∼10⁵ when using multimode fitting as opposed to only fitting (2,±2,n) modes. Our second main result addresses a critical point that has been overlooked in the literature: the importance of matching the Bondi-van der Burg-Metzner-Sachs (BMS) frame of the simulated gravitational wave to that of the QNM model. We show that by mapping the numerical relativity waveforms to the super rest frame, there is an improvement of ∼10⁵ in the all-angles strain mismatch, compared to using the strain whose BMS frame is not fixed. We illustrate the effectiveness of these modeling enhancements by applying them to families of waveforms produced by numerical relativity, and comparing our results to previous studies.