Characterizing gravitational wave signals produced from black hole binaries requires accurately modeling the dynamics of those systems, to construct template waveforms for the inspiral, merger, and ringdown. I will describe recent progress in modeling the dynamics of black hole binaries in modified gravity theories. In particular, I will focus on the development of Numerical Relativity techniques to evolve black holes (and neutron stars) in Einstein scalar Gauss- Bonnet (ESGB) gravity. This modified gravity theory can be motivated by effective field theory reasoning and admits scalar "hairy" black hole solutions. I will discuss how gravitational wave observations can be used to constrain the amount of hair that could exist for black holes in our universe.