Modeling relative peak times of gravitational wave harmonics

Accurate modeling of gravitational waves from binary black hole mergers is essential for extracting their rich physics. A key detail for understanding the physics of mergers is predicting the precise time when the amplitude of the gravitational wave strain peaks, which can differ significantly among the different harmonic modes. We propose two semianalytical methods to predict these differences using the same three inputs from numerical relativity (NR): the remnant mass and spin and the instantaneous frequency of each mode at its peak amplitude. The first method uses the frequency evolution predicted by the backwards one-body model, while the second models the motion of an equatorial timelike geodesic in the remnant black hole spacetime. We compare our models to the Simulating eXtreme Spacetimes Collaboration's waveform catalog for quasicircular, nonprecessing systems and find excellent agreement for $l=|m|$ modes up to $l=8$, with mean and median differences from NR below $1M$ in nearly all cases across the parameter space. We compare our results to the differences predicted by leading effective one-body and NR surrogate waveform models and find that in cases corresponding to the largest timing differences, our models can provide significant increases in accuracy.