Dongze Sun, Béatrice Bonga, Leo C. Stein, Guido Da Re

Post-Newtonian (PN) theory provides the analytic foundation for modeling the early inspiral of binary black holes. However, as an asymptotic series, successive PN orders do not necessarily improve agreement with the full nonlinear dynamics. While this has been explored in the extreme-mass-ratio limit, comparable-mass systems most relevant to current observations have not been benchmarked as systematically at high PN order. We study the convergence of the PN series for non-spinning and quasi-circular systems by comparing the PN energy flux at future null infinity to a long, high-accuracy numerical relativity (NR) simulation. To enable a gauge-consistent comparison, we place both descriptions in the same BMS frame and calibrate the intrinsic PN parameters by fitting to the NR waveform in the early inspiral. We find that for orbital velocities $v\lesssim0.45$, higher PN orders continue to reduce the PN--NR flux discrepancy, with (incomplete) 6PN providing the best agreement among the orders considered. The improvement with PN order is non-monotonic with local extrema around 2.5PN and 4PN. This implies that the optimal truncation order of the PN series cannot be identified from the first local minimum in the energy flux residuals, contrary to suggestions in earlier work. As $v$ approaches $\sim 0.5$ near the innermost circular orbit, higher PN orders no longer improve the agreement between NR and PN, indicating a loss of convergence. These results motivate continued high-order PN calculations and clarify the NR accuracy needed to validate them.