  We present the first systematic comparison between gravitational waveforms
emitted by inspiralling, quasi-circular and nonspinning black hole binaries
computed with three different approaches: second-order gravitational self-force
(2GSF) theory, as implemented in the 1PAT1 model; numerical relativity (NR), as
implemented by the SXS collaboration; and the effective one body (EOB)
formalism, as implemented in the TEOBResumS waveform model. To compare the
models we use both a standard, time-domain waveform alignment and a
gauge-invariant analysis based on the dimensionless function
$Q_\omega(\omega)\equiv \omega^2/\dot{\omega}$, where $\omega$ is the
gravitational wave frequency. We analyse the domain of validity of the 1PAT1
model, deriving error estimates and showing that the effects of the final
transition to plunge, which the model neglects, extend over a significantly
larger frequency interval than one might expect. Restricting to the inspiral
regime, we find that, while for mass ratios $q = m_1/m_2\le 10$ TEOBResumS is
largely indistinguishable from NR, 1PAT1 has a significant dephasing $\gtrsim
1$rad; conversely, for $q\gtrsim 100$, 1PAT1 is estimated to have phase errors
$<0.1$rad on a large frequency interval, while TEOBResumS develops phase
differences $\gtrsim1$rad with it. Most crucially, on that same large frequency
interval we find good agreement between TEOBResumS and 1PAT1 in the
intermediate regime $15\lesssim q\lesssim 64$, with $<0.5$rad dephasing between
them. A simple modification to the TEOBResumS flux further improves this
agreement for $q\gtrsim 30$, reducing the dephasing to $\approx0.27$rad even at
$q=128$. Our results pave the way for the construction of GSF-informed EOB
models for both intermediate and extreme mass ratio inspirals for the next
generation of gravitational wave detectors.
