Why do gravitational waves circularize a binary?

The other answers have given good “rigorous mathematical” arguments for why this happens, but I'd like to add a simpler hand-waving one.

Gravitational waves are emitted when massive bodies accelerate. The acceleration is strongest at periapsis (i.e. when the bodies are closest). The GW emission removes energy. As a result, the orbiting body has less kinetic energy left to rise out of the periapsis, i.e. the next apoapsis will be lower.

By comparison, at the apoapsis there's not so much acceleration, thus the height of the periapsis doesn't change very much. But of course the next time you are at that periapsis, more energy is lost. Repeat, until eventually the apoapsis isn't really higher than the periapsis anymore: you have a circular orbit.

Note that the emission of gravitational waves does not necessarily makes an orbit more circular.

This happens to be the case in the weak field (as detailed in G. Smith's answer). However, for binaries with sufficiently different masses, emission of gravitational waves can actually increase the eccentricity in the strong field regime. This effect was discovered by Glampedakis and Kennefick in gr-qc/0203086, and has since been confirmed by many independent calculations.

Since whether gravitational wave increase or decrease eccentricity, apparently depends on the precise properties of the binary, we should not expect a simple qualitative argument explaining why its does so one way or the other.

The shape of a Keplerian orbit can be characterized geometrically by its semimajor axis $a$ and eccentricity $e$ or dynamically by its energy $E$ and angular momentum $L$.

The latter are expressible in terms of the former as

$$E=-\frac{Gm_1m_2}{2}\frac{1}{a}\tag{1}$$

and

$$L^2=\frac{Gm_1^2m_2^2}{m_1+m_2}a(1-e^2)\tag{2}.$$

The former are expressible in terms of the latter as

$$a=-\frac{Gm_1m_2}{2}\frac{1}{E}\tag{3}$$

and

$$e^2=1+2\frac{m_1+m_2}{G^2m_1^3m_2^3}EL^2\tag{4}.$$

Note that $E$ is negative for a bound orbit.

Gravitational radiation carries energy and angular momentum (and also linear momentum) off to infinity. $E$ decreases and becomes more negative, so its absolute value increases; $L^2$ decreases and becomes less positive, so its absolute value decreases. Whether the absolute value of the negative number $EL^2$ increases or decreases, and thus what happens to the eccentricity, is not obvious.

One must do the calculation! Here is what Peters did.

He first derives/rederives the formulas

$$\frac{dE^\text{rad}}{dt}=\frac{G}{c^5}\left(\frac15\dddot{Q}_{ij}\dddot{Q}_{ij}\right)\tag{5}$$

and

$$\frac{dL_i^\text{rad}}{dt}=\frac{G}{c^5}\left(\frac25\epsilon_{ijk}\ddot{Q}_{jl}\dddot{Q}_{kl})\right)\tag{6}$$

for the rate at which energy and angular momentum are carried to infinity by gravitational waves, in the leading order of a multipole expansion. Here

$$Q_{ij}=\sum_n m^{(n)}\left(x_i^{(n)}x_j^{(n)}-\frac13\delta_{ij}x_k^{(n)}x_k^{(n)}\right)\tag{7}$$

is the system's traceless mass quadrupole moment tensor when the system is considered as $n$ point masses.

He then applies this to a Keplerian binary and averages over one elliptical orbit. Using the conservation of energy and angular momentum, he finds that the binary's energy and momentum decrease at the average rate

$$\left\langle\frac{dE}{dt}\right\rangle=-\frac{32}{5}\frac{G^4}{c^5}\frac{m_1^2m_2^2(m_1+m_2)}{a^5}\frac{1+\frac{73}{24}e^2+\frac{37}{96}e^4}{(1-e^2)^{7/2}}\tag{8}$$

and

$$\left\langle\frac{dL}{dt}\right\rangle=-\frac{32}{5}\frac{G^{7/2}}{c^5}\frac{m_1^2m_2^2(m_1+m_2)^{1/2}}{a^{7/2}}\frac{1+\frac{7}{8}e^2}{(1-e^2)^2}\tag{9}.$$

Differentiating (3) gives

$$\frac{da}{dt}=\frac{Gm_1m_2}{2}\frac{1}{E^2}\frac{dE}{dt}\tag{10}$$

and differentiating (4) gives

$$e\frac{de}{dt}=\frac{m_1+m_2}{G^2m_1^3m_2^3}\left(L^2\frac{dE}{dt}+2EL\frac{dL}{dt}\right)\tag{11}.$$

Substituting (1), (2), (8), and (9) into (10) and (11) gives

$$\left\langle\frac{da}{dt}\right\rangle=-\frac{64}{5}\frac{G^3}{c^5}\frac{m_1m_2(m_1+m_2)}{a^3}\frac{1+\frac{73}{24}e^2+\frac{37}{96}e^4}{(1-e^2)^{7/2}}\tag{12}$$

and

$$\left\langle\frac{de}{dt}\right\rangle=-\frac{304}{15}\frac{G^3}{c^5}\frac{m_1m_2(m_1+m_2)}{a^4}\frac{e(1+\frac{121}{304}e^2)}{(1-e^2)^{5/2}}\tag{13}.$$

You can see that the rate of decrease of the eccentricity is very rapid for a highly eccentric orbit with $e$ near $1$, due to the $(1-e^2)^{5/2}$ in the denominator. In others words, the orbit rapidly circularizes.

From these equations, Peters proceeds to find $a$ as a function of $e$ (with two unusual exponents, $12/19$ and $870/2299$) and a differential equation for $e(t)$ from which the lifetime of the binary can be found.