 


When the black holes are within a few radii of each other, their event horizons become distorted, and they enter the merger phase. They fuse together in conditions of extreme gravity to form a single, highly perturbed black hole. Elucidating the dynamics of this process and calculating the accompanying waveforms cannot be accomplished analytically. Rather, we must use numerical relativity, in which the Einstein equations are solved on a computer, in three spatial dimensions plus time.  
Finally, the distorted remnant black hole evolves toward a symmetrically rotating state—known as a Kerr black hole—by emitting its perturbations in the form of gravitational waves. We say that the black hole "rings down," in analogy to the way a ringing bell emits its distortions as sound waves. We can solve the equations of general relativity during this ringdown phase analytically by considering the distortions to be perturbations of a Kerr black hole; the ringdown waveforms are sinusoids whose amplitudes drop off exponentially.  Detecting the gravitational waves of such mergers would be the scientific equivalent of striking gold. 

Simulating a black hole merger on a computer and calculating the
resulting gravitational waveforms was long considered the holy grail of
numerical relativity, both for the difficulty of this endeavor and the
importance of the results. From the perspective of fundamental physics,
black hole mergers are the ultimate "twobody" problem in general
relativity. The solution to the equations can teach us a great deal
about the effects of extreme gravity on the surrounding spacetime.
Black hole mergers also provide superb laboratories for testing general
relativity in very strong and dynamical fields. 

Seeking the Holy Grail In 1964 Susan Hahn and Richard Lindquist made the first attempt to solve Einstein's equations on a computer. They attempted to simulate the headon collision of two equalmass, nonspinning black holes. Since the term "black hole" had not yet been coined, they actually tried to model the headon collision of two "wormholes." While this valiant attempt was well ahead of its time, the computer code crashed before the wormholes could collide. 

In the 1970s, Larry Smarr and Kenneth Eppley first successfully simulated a headon collision of two black holes using supercomputers at Lawrence Livermore National Laboratory (LLNL). In addition to LLNL's computer power, Smarr and Eppley benefited from advances in understanding how to write Einstein's equations in a form suitable for integration on a computer. They solved the Einstein equations in two spatial dimensions (2D) plus time (that is, symmetric about an axis), and obtained information for the first time about the gravitational waves. They also produced the first numerical relativity scientific visualizations. 
Simulating a black hole merger on a computer and calculating the
resulting gravitational waveforms was long considered the holy grail of
numerical relativity, both for the difficulty of this endeavor and the
importance of the results. 

The next major step would be to simulate orbiting black holes, which requires solving Einstein's equations in three spatial dimensions (3D) plus time. But in the 1970s and 1980s, there were too many technical difficulties with the simulation codes and insufficient computer power available. While some work continued on evolving binary neutron stars, no further progress was made.  
The 1990s brought renewed interest in black hole mergers. As plans for LIGO developed, physicists recognized these mergers as major sources of strong gravitational waves. In response to the growing experimental effort to detect gravitational waves, the NSF funded a multiyear, multiinstitution Grand Challenge grant in the mid1990s aimed at simulating black hole mergers. The several groups involved in the Grand Challenge built largescale 3D codes, but found their evolutionary models beset by a host of instabilities. While they could simulate grazing black hole collisions, their codes crashed due to instabilities before any significant fraction of a binary orbit could be computed.  
After the Grand Challenge grant ended, the participating research groups continued working individually. A major focus was the formalisms, that is, the way the Einstein equations are expressed mathematically, including the choice of variables. They learned that instabilities from unphysical solutions to the equations were causing the codes to crash. Numerical simulations excite these modes, which grow exponentially and destroy the computation. Researchers devised strategies to rewrite Einstein's equations in different ways, to eliminate these unphysical solutions. Overall, progress was made, but incrementally.  
By the early 2000s, many in the community felt a sense of despair, with
some saying that "numerical relativity is impossible." This turned out
to be the darkest hour before the dawn, however, because a stunning
series of breakthroughs was about to begin. 

Stunning Advances In late 2003, a group led by Bernd Brügmann (now at Fredrich Schiller University in Jena, Germany) achieved the first complete orbit of a black hole binary. The team's code expressed the Einstein equations in a form known to be conducive to numerical stability. They represented the black holes as "punctures," a clever method for handling black hole singularities (where the spacetime curvature becomes infinite) on a numerical grid. 

Their specific approach, which was the basis for most work in numerical relativity at the time, required that the punctures remain fixed in the grid so that their singularities could be easily factored out. But since orbiting black holes actually move, Brügmann's group used a coordinate system in which the black holes are fixed in place, and the spacetime revolves around them. The team succeeded in simulating the binary for a little more than one orbit...before its code also crashed.  
In early 2005, Frans Pretorius (now at Princeton University) succeeded in simulating the first orbit, merger, and ringdown of equalmass black holes. Rather than using punctures, he handled the singularities by excising the regions inside the event horizons from the computational domain. These black holes were not fixed in place, but rather moved on their orbits through the numerical grid. His code was based on a very different way of writing the Einstein equations, called the generalized harmonic form. Pretorius' code did not crash, and he extracted the first gravitational waveform from a black hole merger.  
At this time, most of the numerical relativity community had codes based on the traditional techniques used by Brügmann's group, including puncture black holes. Pretorius' success led to the natural question of whether stable computations of black hole mergers required his very different techniques instead of the traditional puncture approach.  
The answer came later in 2005 with the moving puncture method, which was discovered simultaneously and independently by the numerical relativity group led by Manuela Campanelli (now at the Rochester Institute of Technology) and our group at NASA's Goddard Space Flight Center. This method is based on the traditional approach to evolving puncture black holes with Einstein's equations, but with a simple and crucial difference.  
When solving the Einstein equations, we are free to choose our coordinate systems; in fact, the coordinates we use can themselves change as the simulation proceeds, in order to keep the equations in a nice form for numerical computation. The moving puncture method is based on a specific choice of coordinates that enables the punctures to move across the grid. Using moving punctures, both the Campanelli and Goddard groups achieved accurate simulations of the orbit, merger, and ringdown of equal mass, nonspinning black holes and extracted the gravitational waveforms—with no code crashes.  
A
series of stunning advances resulted from this new moving puncture
method. Our Goddard group in particular was quickly able to accomplish
multiple binary orbits and uncovered the universal waveform signatures
of nonspinning, equalmass black holes. Since the moving puncture
method required only simple changes to existing numerical relativity
codes based on traditional puncture techniques, other groups worldwide
quickly adopted this method and got into the game. Within a few months,
these teams simulated mergers of unequal mass and of spinning black
holes, and their results were applied to astrophysics. A true golden
age of black hole simulations had opened, and the fun was only just
beginning! 

Anatomy of a Merger Simulation Now let's take an indepth look at how our group calculates black hole mergers on a computer. We start with the black holes on a computational grid at some initial time. The Einstein equations of general relativity are then used to evolve the binary forward in time. What could be simpler? 

A key issue in the simulation is the extreme spatial curvature inside the black holes, which approaches infinity at the centers and could thus cause computers to crash. Fortunately, things inside the black hole horizon cannot escape and influence the spacetime outside, so we do not have to include these physical infinities in our simulation. Instead, we represent each black hole as a puncture, with all the nastiness of the infinite curvature written in terms of factors taking the form (1/r)p, where p > 0. Here r is the distance from the center of the black hole, and these terms become infinite as the center is approached, r → 0. In the traditional puncture method, these (1/r)p quantities—the infinities—are factored out and dealt with separately, by hand. The computers then work only on finite quantities, and the infinite quantities are factored back in only when needed.  A series of stunning advances resulted from this new moving puncture method. 

Isolating the infinite curvatures in this manner requires an extra condition on the coordinate system to keep each puncture at a fixed location in the grid. But the black holes are orbiting each other and spiraling inward. If we evolve the binary with the black holes fixed in the coordinate system, the computational grid twists itself up and becomes highly stretched and warped. This causes various problems that can cause computer codes to crash, motivating researchers to try different strategies.  
In the summer of 2005, our group was experimenting with the puncture method and ran a simulation in which we did not isolate and factor out the infinities. Rather, we let the computer deal with the resulting large numbers. To our surprise, the code did not crash! Instead, the discrete nature of the computation cut off the infinities at the punctures while preserving the mathematical integrity of the black holes around the punctures. Campanelli's group independently came to a similar realization at about the same time.  


This discovery led to the development of a different approach. We found that if we did not isolate and factor out the infinities, we could allow the punctures to move across the grid. To ensure that the punctures move smoothly, and to produce a stable and accurate simulation, we needed a very careful choice of coordinate system. It took us several months to develop this system, but it allowed us to simulate multiple binary orbits and the final merger—our longsought breakthrough.  
Black hole merger simulations present other technical challenges. For example, the gravitational radiation produced by the merger has wavelengths typically 10 to 100 times larger than the radius of the horizons. For accurate models, we need to resolve both the strong gravitational fields near the black holes and the radiation in the wave zone. We accomplish this by using variable grid resolution in our runs, with finer grids near the black holes and coarser ones in the outer regions where the waves propagate (figure 6).  
Because the black holes move freely through the grid, this mesh spacing must adapt to the strength of the fields. In the region where the black holes are orbiting, we base our refinement criteria on a curvaturerelated quantity that takes large values near the black hole and then decreases away from the source. When the value of this quantity exceeds a chosen threshold, the code doubles the number of grid zones in each direction to increase the resolution in that region. In the wave zone and beyond we use progressively coarser, fixed mesh refinement; this also makes it computationally inexpensive to push the outer boundary of the simulated space far enough away that we need not worry about reflections of the waves from the outer edges of the grid.  
Despite these advances, the numerous orbits we wish to simulate
typically demand thousands of CPU hours. On the Columbia system at
NASA's Ames Research Center, as well as the Discover cluster here at
Goddard, we often run on as many as 500 processors for up to 70 hours.
As we push to binaries with higher mass ratios and spins, which require
greater accuracy, we expect the larger Jaguar system at Oak Ridge
National Laboratory to provide additional computational power. 

Gravitational Waveforms Once we had this moving puncture method working well, we first aimed to compute the definitive gravitational wave signal for the simplest case: the merger of equalmass nonspinning black holes. For simulations of astrophysically realistic binaries, the black holes should ideally start out on nearly circular orbits around their common center of mass. However, the types of initial conditions being used at that time were imperfect because they had some orbital eccentricity. 

To test the robustness of our models, we evolved equalmass, nonspinning binaries starting from several different separations, from roughly two to four orbits before merger. Despite the differing amounts of initial orbital eccentricity, we found that in every run the black holes locked on to the same trajectories roughly one orbit before merger.  


Figure 7 shows the gravitational waves emanating from the merger of an equalmass binary. We can calculate the waveform by averaging this radiation over an imaginary sphere centered on the source, and graphing the result as a function of time. When we compared the waveforms from our runs we found that the early stages of the waveforms differ by about 10%, corresponding to the different amounts of initial eccentricity. But for the final part of the waveform that corresponds to the last orbit and merger, the waveforms were nearly identical! We thus obtained a clear and universal waveform, independent of the initial starting point. This signal is shown as panel (a) in figure 8.  


So, despite our original expectations, the simulations showed that the uncertainty in the initial conditions did not significantly affect our ability to predict the gravitational waveforms expected from real mergers. Likewise, when we compared our waveforms with those produced by different groups using different techniques, we found good agreement. The field of gravitational waveform simulation was suddenly open for business.  
The waveforms were surprising in another way as well. There is an important difference in how gravitational radiation is generated in the early part of the waveform compared to its end. At first the waveform represents an inspiraling orbital system, while later it represents the "ringing" of the distorted single black hole produced in the merger. Are there dramatic features present in the waveforms that demark this seemingly important change in the underlying physics? For years, scientists had assumed that the complex and nonlinear interactions of strongfield sources near merger would leave strong imprints on the waveforms.  
Our simulations revealed quite a different result. The wave frequency initially "chirps" or increases smoothly and monotonically up to an expected peak, and after merger it matches the expected final frequency for the ringing of the final distorted black hole. Likewise, the wave amplitude increases smoothly to a peak before decaying as a damped sinusoid in the ringdown. As can be seen in figure 8 (a, p21), the waveforms are remarkable only in their simplicity, with the merger producing a smooth transition between the inspiral and ringdown. The nonlinearities of general relativity have conspired, it seems, to produce the simplest possible transition.  
Nevertheless, the gravitational waveforms from all mergers do not look
alike. Rather, the detailed shapes of the waveforms depend on the mass
of each black hole, their spins, and the orientations of their
rotational axes. Black hole binaries with unequal masses produce more
complicated distortions of the surrounding spacetime, resulting in more
complex waveforms from most viewing positions (figure 8, b, p21). And
mergers of spinning black holes bring additional interesting physics.
For example, compared to the nonspinning case, mergers of black holes
with their spins aligned with the orbital plane tend to "hang up,"
taking longer to plunge inward and merge. This leads to more densely
packed, higherfrequency wavefronts (figure 8, c, p21) . 

Recoil Kicks In the real Universe, black holes in binaries will almost always have unequal masses. Such a binary will emit gravitational radiation, and the momentum it carries, preferentially in the direction of the smaller hole. Early on, the black holes are widely separated and this asymmetry does not amount to much, when averaged over the course of an orbit. But in the final stages of inspiral, the radiation amplitude—and hence the amount of momentum loss—increases. When the binary enters its final plunge to merger, the system has run out of time to balance itself. 

The result is something we can predict from Newton's Third Law of Motion: since the binary emits more radiation—and hence momentum—in one direction, it must compensate by moving in the opposite direction. The more momentum it loses, the larger the "recoil kick" of the postmerger black hole. Because this effect depends so strongly on the endstage merger, it was not until the recent breakthroughs that the various research groups could provide a definitive answer to the size of such a recoil: the maximum kick is about 175 km s^{1}, for a binary mass ratio of about 3 to 1. 
And these advances are occurring across a broad front, with many
research groups carrying out calculations with independently written
codes. 

Now let's add spins to the mix. When the spins are perpendicular to the orbital plane and spinning in the same direction as the orbit, which is what astronomers expect when supermassive black holes are surrounded by gaseous accretion disks, the total kick can climb to about 500 km s^{1} for maximally spinning holes. In this case, the largest kick occurs for a binary with a mass ratio closer to 2 to 1.  


When the inspiraling holes spin in arbitrary directions, life becomes much more complicated: the spins precess, the orbits rise and fall, and the final recoil kick can be much larger. In fact, several research groups have found that some configurations may produce recoils up to 4,000 km s^{1}—easily enough to eject the merged hole from the center of even a large galaxy.  
Figure 9 shows the distribution of gravitationalwave energy emitted from the merger of two equalmass, spinning holes. The visible northsouth asymmetry gives rise to a "recoil kick" upwards.  
As you might imagine, this result has sparked the interest of
astronomers. If such kicks are commonplace, then many large galaxies
should no longer contain massive central holes. Given the apparent
robustness of the numerical results, and the ubiquity of massive
central black holes in galaxies, the focus has now shifted to
investigating the premerger configurations of the binary systems.
Physics has given way to astrophysics. 

Looking Ahead The recent progress in numerical relativity simulations of black hole mergers is both exciting and impressive by any measure. And these advances are occurring across a broad front, with many research groups carrying out calculations with independently written codes. 

For the simplest case of equalmass, nonspinning black holes, there is general agreement on the basic results. The waveform takes the simple shape shown in figure 8 (a, p21). The final black hole has a spin roughly 70% of the maximum value allowed by general relativity. And the total energy carried away by the gravitational waves during the merger is roughly 4% of the total energy of the binary. This yields an energy output of 10^{23} solar luminosities, with more massive black holes radiating this energy over longer timescales.  The mergers seen by LISA will be very energetic, and are expected to be detected with a very large signaltonoise ratio. 

Long runs of black hole mergers are now being carried out, starting roughly 10 orbits before the merger. This is a key development, because such waveforms are needed as templates to search for signals in the data from gravitational wave detectors. The groundbased detectors have already taken some early data, and efforts are underway to incorporate the numerical relativity results into their analysis. The availability of merger waveforms will increase LISA's sensitivity to black hole mergers reaching back very early in cosmic history.  
The mergers seen by LISA will be very energetic, and are expected to be detected with a very large signaltonoise ratio. Figure 10 shows how a "pure" waveform of the sort shown in figure 8 (a, p21) would appear in the LISA detector data stream.  


And
finally, there is an explosion of work on mergers of unequal mass and
spinning black holes and the resulting recoil kicks. Studies of the
waveforms and final spins of the black holes are underway. The kicks
especially have important astrophysical implications that are now under
study. We live in interesting and exciting times. Stay tuned! 

Contributors:
This article was written by Joan Centrella, Bernard Kelly, Jim Van
Meter, John Baker, and Robin Stebbins of the Gravitational Wave
Astrophysics Laboratory at NASA's Goddard Space Flight Center, and by
Goddard science writer Robert Naeye. The group acknowledges Chris Henze
of NASA's Ames Research Center for supplying images, and the use of the
supercomputers Columbia at NASA Ames and Discover at NASA Goddard.

