  In this work we study the dynamics of spinning binary black hole systems in
the strong field regime. For this purpose we extract from numerical relativity
simulations the binding energy, specific orbital angular momentum, and
gauge-invariant orbital frequency. The goal of our work is threefold: First, we
extract the individual spin contributions to the binding energy, in particular
the spin-orbit, spin-spin, and cubic-in-spin terms. Second, we compare our
results with predictions from waveform models and find that while
post-Newtonian approximants are not capable of representing the dynamics during
the last few orbits before merger, there is good agreement between our data and
effective-one-body approximants as well as the numerical relativity surrogate
models. Finally, we present phenomenological representations for the binding
energy for non-spinning systems with mass ratios up to $q = 10$ and for the
spin-orbit interaction for mass ratios up to $q = 8$ obtaining accuracies of
$\lesssim 0.1\%$ and $\lesssim 6\%$, respectively.
