  We recently derived, using the density-of-states approximation, analytic
distribution functions for the outcomes of direct single-binary scatterings
(Stone & Leigh 2019). Using these outcome distribution functions, we present in
this paper a self-consistent statistical mechanics-based analytic model
obtained using the Fokker-Planck limit of the Boltzmann equation. Our model
quantifies the dominant gravitational physics, combining both strong and weak
single-binary interactions, that drives the time evolution of binary orbital
parameter distributions in dense stellar environments. We focus in particular
the distributions of binary orbital energies and eccentricities. We find a
novel steady state distribution of binary eccentricities, featuring strong
depletions of both the highest and the lowest eccentricity binaries. In energy
space, we compare the predictions of our analytic model to the results of
numerical N-body simulations, and find that the agreement is good for the
initial conditions considered here. This work is a first step toward the
development of a fully self-consistent semi-analytic model for dynamically
evolving binary star populations in dense stellar environments due to direct
few-body interactions.
