  Binary black holes (BBHs) in eccentric orbits produce distinct modulations
the emitted gravitational waves (GWs). The measurement of orbital eccentricity
can provide robust evidence for dynamical binary formation channels. We analyze
57 GW events from the first, second and third observing runs of the
LIGO-Virgo-KAGRA (LVK) Collaboration using a multipolar aligned-spin
inspiral-merger-ringdown waveform model with two eccentric parameters:
eccentricity and relativistic anomaly. This is made computationally feasible
with the machine-learning code DINGO which accelerates inference by 2-3 orders
of magnitude compared to traditional inference. First, we find eccentric
aligned-spin versus quasi-circular aligned-spin $\log_{10}$ Bayes factors of
1.84 to 4.75 (depending on the glitch mitigation) for GW200129, 3.0 for
GW190701 and 1.77 for GW200208_22. We measure $e_{\text{gw}, 10Hz}$ to be
$0.27_{-0.12}^{+0.10}$ to $0.17_{-0.13}^{+0.14}$ for GW200129,
$0.35_{-0.11}^{+0.32}$ for GW190701 and $0.35_{-0.21}^{+0.18}$ for GW200208_22.
Second, we find $\log_{10}$ Bayes factors between the eccentric aligned-spin
versus quasi-circular precessing-spin hypothesis between 1.43 and 4.92 for
GW200129, 2.61 for GW190701 and 1.23 for GW200208_22. Third, our analysis does
not show evidence for eccentricity in GW190521, which has an eccentric
aligned-spin against quasi-circular aligned-spin $\log_{10}$ Bayes factor of
0.04. Fourth, we estimate that if we neglect the spin-precession and use an
astrophysical prior, the probability of one out of the 57 events being
eccentric is greater than 99.5% or $(100 - 8.4 \times 10^{-4})$% (depending on
the glitch mitigation). Fifth, we study the impact on parameter estimation when
neglecting either eccentricity or higher modes in eccentric models. These
results underscore the importance of including eccentric parameters in the
characterization of BBHs for GW detectors.
