According to general relativity, binary stellar systems emit gravitational radiation, and thereby lose energy and angular momentum. This dissipation causes the orbit to circularize and decay, with ever decreasing orbital period and orbital separation. Eventually, systems containing compact objects such as neutron stars or black holes, will decay to such tight, high-velocity orbits that a final burst of gravitational radiation will result, accompanied by a rapid inspiral toward coalescence. This final ``death-spiral'' of compact binaries is viewed as the most promising source of gravitational waves for kilometer-class laser-interferometric gravitational observatories which have recently begun science operations. These include LIGO (USA), VIRGO (France-Italy), and GEO600 (UK-Germany). The proposed space interferometer LISA will also open the low-frequency gravitational-wave window. Detection and study of the gravitational waves from such systems will initiate a new era of gravitational-wave astronomy and provide important new tests of general relativity, especially in its highly dynamical, strong-field regimes.

Extracting the characteristic inspiralling-binary signal from the noisy data in an array of detectors, and using the signal to obtain useful information about the sources, such as their directions, distances, masses, spins, and so on will involve a process of ``optimal matched filtering'' of a set of theoretical ``template waveforms'' against the output of the detectors. This method requires that the theoretical template be a sufficiently accurate representation of the waves from the hypothesized source that it match the observed signal to within a fraction of a phase over the hundreds to thousands of gravitational-wave cycles observed in the sensitive bandwidth.

For most of the binary inspiral phase, (v/c)^2 ~ (Gm/rc^2) is small, and so post-Newtonian (PN) approximations to general relativity can be used to describe the orbit and gravitational waveform. However, in order to calculate a template of sufficient accuracy, higher-order post-Newtonian corrections must be taken into account. Current conventional wisdom says that corrections through at least third PN (3PN) order should be calculated.

Most of our calculations of higher-order post-Newtonian effects have been carried out using a method called DIRE, Direct Integration of the Relaxed Einstein equations. It is based on the 1975 Epstein-Wagoner (EW) approach, in which Einstein's equations are cast in their so-called ``relaxed'' form, as a flat spacetime wave equation with a non-compact source consisting of matter and the stress-energy of the gravitational fields themselves. The solutions are expressed as integral equations, which can be iterated order-by-order. In the near zone, we evaluate the integrals by a slow-motion, multipole expansion. In the far zone, we evaluate the integrals over the past-light cone using a technique discovered by our group (Will and Wiseman 1996, Pati and Will 2000). In this new method, there are NO divergent integrals, at least through 3.5 PN order. Furthermore, the effects of gravitational-wave ``tails'' and of propagation of the radiation along the true cones of the spacetime are obtained.

Significant progress has been made in the past few years in calculating and studying the implications of higher-order post-Newtonian effects in the behavior of and gravitational waves from inspiralling binaries (see, for example, the review by Blanchet). Among the recent results achieved by members of WUGRAV are: