{'fraction_non_alphanumeric': 0.02946078348849346, 'fraction_numerical': 0.008669887712328076, 'mean_word_length': 4.531427507216687, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulations of Binary Black Hole Mergers Using Spectral Methods

Several improvements in numerical methods and gauge choice are presented that make it possible now to perform simulations of the merger and ringdown phases of"generic"binary black-hole evolutions using the pseudo-spectral evolution code SpEC. These improvements include the use of a new damped-wave gauge condition, a new grid structure with appropriate filtering that improves stability, and better adaptivity in conforming the grid structures to the shapes and sizes of the black holes. Simulations illustrating the success of these new methods are presented for a variety of binary black-hole systems. These include fairly ``generic'' systems with unequal masses (up to 2:1 mass ratios), and spins (with magnitudes up to 0.4 M^2) pointing in various directions.

These first results on binary black-hole systems were obtained by several different groups using different codes based on two different formulations of the Einstein equations (BSSN and generalized harmonic), using two different methods for treating the black-hole interiors (moving puncture and excision), and using rather different gauge conditions to fix the spacetime coordinates. All of these early results, however, were obtained with codes based on finite difference numerical methods and adaptive mesh refinement for computational efficiency. The Caltech/Cornell collaboration decided a number of years ago to follow a different path by developing an Einstein evolution code, called SpEC, based on spectral methods. The advantages of spectral methods are their superior accuracy and computational efficiency, and their extremely low numerical dissipation that is ideal for solving wave propagation problems with high precision. The disadvantages are the relative complexity of spectral codes, the lack of any appropriate pre-existing spectral code infrastructure (analagous to CACTUS [43] for example), and the extreme sensitivity of spectral algorithms to develop-ing instabilities when any aspect of the solution method is mathematicaly ill-posed. Consequently it has taken our group somewhat longer to bring SpEC up to the level of the state-of-the-art codes in the field.
We, along with our Caltech/Cornell collaborators, have developed a large number of numerical and analytical tools over the years that make it possible for SpEC to evolve binary black-hole systems with greater precision than any other code (at the present time) [44]. These technical developments include the derivation and implementation of constraint preserving and physical (no incoming gravitational wave) boundary conditions [45,46], dual frame evolution methods and feedback and control systems to lock the computational grids onto the location of the black holes [47], and special angular filtering methods needed to cure an instability that occurs when tensor fields are evolved [48]. Using these methods we and our collaborators have performed a number of high precision evolutions of the inspiral portions of binary black-hole systems, and have used the gravitational waveforms from these evolutions to calibrate the accuracy of the approximate post-Newtonian waveforms that are widely used in gravitational wave data analysis [49,50,51,52].
During the past year our group has begun to have some success in performing the more dynamical and difficult (for spectral codes) merger and ringdown portions of binary black-hole evolutions [50,52,53,54,55,56]. The techniques we used to achieve these first merger calculations turned out to be rather non-robust however, requiring a great deal of hands-on adjustment and fine tuning for each new case we attempted. In this paper we report on several new numerical and analytical developments that allow us now to perform stable and accurate binary merger and ringdown simulations of a fairly wide range of binary black-hole systems. These new methods appear to be quite robust, and no fine tuning is required: the same basic method works on each of the cases we have attempted.
The three new technical breakthroughs that allow us now to perform successful binary black-hole merger and ringdown simulations are described in Sec. II. These major advances include the development of a new gauge condition that works extremely well with the generalized harmonic formulation of the Einstein system used in our code. This new gauge chooses spatial coordinates that are solutions of a damped wave equation, and chooses the time coordinate in a way that limits the growth of √ g/N , here g is the determinant of the spatial metric and N is the lapse function. This new gauge condition is described in some detail in Sec. II A. Another important development, reported in Sec. II B, is the construction of a new grid structure on which our binary black-hole evolutions are performed. This new non-overlapping grid structure, and a certain type of spectral filtering used with the new grid structure, removes a class of numerical instabilities that were a limiting factor in our ability to perform the highly dynamical portions of merger calculations, and allows us to increase resolution near the black holes in a more targeted and computationally efficient way. The third technical development is a new more efficient and robust method to conform the structure of the grid to the shape and size of the black holes in an automatic dynamical way. These developments are described in Sec. II C and the Appendix .
Using these new technical tools we have now performed successful merger and the succeeding ringdown simulations of a fairly wide range of binary black-hole systems. We describe in Sec. III merger simulations for six reasonably diverse cases. These include two systems in which the black holes are non-spinning: one has equal mass black holes, the other has holes with a 2:1 mass ratio. We also describe successful merger and ringdown calculations with two different equal-mass binary systems in which the holes have identical intrinsic spins of magnitude 0.4M 2 : aligned to the orbital angular momentum in one system, and anti-aligned in the other. We have also performed successful merger and ringdown simulations on two fairly generic binary systems. These two systems have black holes with the same mass ratio M 2 /M 1 = 2, and the same "randomly oriented" initial spins of magnitudes 0.2M 2 1 and 0.4M 2 2 respectively. But the initial separations of the holes are different in the two cases: one starts about 8.5 orbits before merger, and another (used to perform careful convergence studies) starts about 1.5 orbits before merger. The same methods were used to perform all of these merger and ringdown simulations, and no particular fine tuning was required.

II. TECHNICAL DEVELOPMENTS
This section describes the three new technical breakthroughs that allow us now to perform successful binary black-hole merger and ringdown simulations. These major advances include the development of a new gauge condition, described in Sec. II A, in which the spatial coordinates satisfy a damped wave equation and the time coordinate is chosen in a way that controls the growth of the spatial volume element. Another important development, described in Sec. II B, is a new non-overlapping grid structure for our binary black-hole simulations, and a type of spectral filtering for these new grid structures that improves their accuracy and stability. The third new technical development, described in Sec. II C, is a more efficient and robust method of conforming the structure of the grid to the shape and size of the black holes in a dynamical and automatic way.

A. Damped-Wave Gauge
Harmonic gauge is defined by the condition that each coordinate x a satisfies the co-variant scalar wave equation: Harmonic coordinates have proven to be extremely useful for analytical studies of the Einstein equations, but have found only limited success in numerical problems like simulations of complicated highly dynamical blackhole mergers. One reason for some of these difficulties is the wealth of "interesting" dynamical solutions to the harmonic gauge condition itself, Eq. (1). Since all "physical" dynamical fields are expressed in terms of the coordinates, an ideal gauge condition would limit coordinates to those that are simple, straightforward, dependable, and non-singular; having "interesting" dynamics of their own is not a desirable feature for coordinates. The dynamical range available to harmonic coordinates can be reduced by adding a damping term to the equation: [57,58] where t a is the future directed unit normal to the constant-t hypersurfaces. Adding such a damping term to the equations for the spatial coordinates x i tends to remove extraneous gauge dynamics and drives the coordinates toward solutions of the co-variant spatial Laplace equation on the timescale 1/µ. Choosing 1/µ to be comparable to (or smaller than) the characteristic timescale of a particular problem should remove any extraneous coordinate dynamics on timescales shorter than the physical timescale. The addition of such a damping term in the time-coordinate equation is not appropriate however. Such a damped-wave time coordinate is driven toward a constant value, and therefore toward a state in which it fails to be a useful time coordinate at all. It makes sense then to use the damped-wave gauge condition only for the spatial coordinates: where N i is the shift, and N is the lapse. The appropriate contra-variant version of this damped-wave gauge condition is therefore where g ab is the spatial metric of the constant-t hypersurfaces. (Note that we use Latin letters from the beginning of the alphabet, a, b, c, . . ., to denote four-dimensional spacetime indices, and letters from the middle of the alphabet, i, j, k, . . ., for three-dimensional spatial indices.) While the damped-wave gauge is a poor choice for the time coordinate, the idea of imposing a gauge that adds dissipation to the gauge dynamics of the time coordinate is attractive. To find the appropriate expression for t a H a , the component of H a not fixed by Eq. (4), we note that the gauge constraint H a + Γ a = 0 implies that t a H a is given by where g = det g ij is the spatial volume element. In our experience, a frequent symptom of the failure of simpler gauge conditions in binary black-hole simulations is an explosive growth of g in the spacetime region near the black-hole horizons. This suggests that a good use of the remaining gauge freedom would be to attempt to control the growth of the spatial volume element, g. Choosing the gauge condition, together with Eq. (5) implies the following evolution equation for √ g/N : whose solutions tend to suppress growth in √ g/N . The discussion of this gauge condition in Ref. [57] shows that it also implies that the lapse N satisfys a damped wave equation, with the damping factor µ L . So in this sense, the gauge condition on the time coordinate, Eq. (6), is the natural extension of the spatial-coordinate dampedwave gauge condition, Eq. (4).
Combining this new lapse condition, Eq. (6), with the damped-wave spatial coordinate condition, Eq. (4), gives the gauge-source function for our full damped-wave gauge condition: The damping factors µ L ≥ 0 and µ S ≥ 0 can be chosen quite arbitrarily as functions of spacetime coordinates x a , or even as functions of the spacetime metric ψ ab . The gauge source function H a depends only on coordinates and the spacetime metric ψ ab in this case, so these gauge conditions can be implemented directly in the generalized harmonic Einstein system without the need for a gauge driver. Previous studies of this condition (developed initially as a test of the first-order gauge-driver system [57]) showed it to be quite useful for evolving single black-hole spacetimes. In those tests, which included several different evolutions of maximal-slice Schwarzschild initial data with large non-spherical gauge perturbations, the black hole always evolved quickly toward a non-singular timeindependent equilibrium state. The solutions to the lapse gauge condition, Eq. (7), can be thought of as equating log( √ g/N ) to a certain weighted time average of ∂ k N k /N . The timescale associated with this time averaging is set by µ L , which determines for example the rate at which √ g/N is driven toward an asymptotic equilibrium state. When µ L is constant, it is easy to show that √ g/N is driven exponentially toward this asymptotic state. In highly dynamical spacetimes, however, we find that √ g/N must be driven even faster than exponential in order to prevent the formation of singularities in g. This can be accomplished by making µ L larger whenever g becomes large. In practice we find that the choice is a constant or perhaps a function of time) is very effective at suppressing the growth of these singularities. We find that choosing the same damping factor µ S , for the spatial part of the gauge condition is also quite effective. The binary black-hole merger and ringdown simulations described in Sec. III use these µ L and µ S , with µ 0 taken to be an order-unity function of time (to accomodate starting up evolutions from initial data satisfing a different gauge condition).

B. Grid Structure and Filtering
The pseudo-spectral numerical methods used by our code, the Spectral Einstein Code (SpEC) [45,47,59], represent dynamical fields on a spatial grid structure that is specially constructed for each problem. For binary blackhole simulations, our group has been using grid structures constructed from layers of spherical-shell subdomains centered on each black hole, surrounded and connected by cylindrical, cylindrical-shell, and/or rectangular-block subdomains, all surrounded by spherical-shell subdomains that extend to the outer boundary of the computational domain (located far from the holes). The intermediate-zone cylindrical-shell and/or rectangularblock subdomains must overlap both the inner and outer spherical-shell subdomains in these grid structures to cover the computational domain completely. These grid structures are quite efficient, and have allowed our group to perform long stable inspiral calculations for a variety of simple binary black-hole systems [47,49,52,60], and also simulations of the merger and ringdown phases of a few of these simple cases [53,56]. The overlap regions in these grid structures are well behaved in these successful cases. But in more generic inspiral and merger simulations, these overlap regions become significant sources of numerical error and instability.
The overlap-related instability discussed above is most likely caused by the method we use to exchange information between adjacent subdomains. We set the incoming characteristic fields, whose values are determined by interpolation from the adjacent subdomain, using a penalty method [61,62,63,64]. This penalty method for imposing boundary conditions was derived explicitly for the case where adjacent subdomain boundaries touch but do not overlap. One possibility for resolving this instability would be to re-derive the appropriate penalty boundary terms for the case of overlapping subdomains. We have chosen the simpler option of constructing new grid structures without overlapping subdomains. We do this by changing the intermediate-zone grid structure from the cylindrical-shells and/or rectangular blocks used previously into a series of deformed rectangular blocks whose boundaries conform to the spherical shells used near the inner and outer boundaries. The new type of grid element needed for this construction is illustrated in Fig. 1, which shows how the 3D volume between a sphere and a concentric cube can be mapped into six deformed cubes. This construction is based on the "cubed-sphere" coordinate representations of the sphere (i.e. the mapping of the six faces of the cube onto the sphere as illustrated in Fig. 1), that is widely used in atmospheric and geophysical modeling [65,66,67,68]. Three-dimensional grid structures based on cubed spheres have also been used in other types of simulations, including a few black-hole simulations [69,70].
Our new grid structure is illustrated in Figs. 2 and 3. It consists of a series of spherical-shell subdomains and several layers of cubed-sphere subdomains surrounding each black hole, as illustrated in Fig. 2. The two cubic blocks containing the black holes are surrounded by a series of rectangular-block subdomains to fill out the volume of a large cube. The center of this large cube is placed at the center-of-mass of the binary system, as illustrated in Fig. 3 for a system having a 2:1 mass ratio. This large cube is connected to the outer spherical-shell This collection of blocks filles up a cube which is surrounded by a series of cubed-sphere layers (two shown) all surrounded by a series of spherical shells (one shell shown) that are centered on the center of mass of the binary system. subdomains by several additional layers of cubed-sphere subdomains, also illustrated in Fig. 3.
Previous attempts to use this type of cubed-sphere grid structure in SpEC have always proven unsuccessful, due to numerical instabilities at the interdomain boundaries. These previous attempts used no spectral filtering, or a two-thirds anti-aliasing filter, in the cubed-sphere subdomains. It has been shown, however, that an appropriate filter is required for stability (and improved accuracy) of Chebyshev polynomial spectral expansions, such as those used in these cubed-sphere subdomains. Such a filter must satisfy two important conditions. The first criterion is that when represented in physical space, this filter must act like a dissipation term in the evolution equations that vanishes on the subdomain boundaries (to avoid conflicting with the boundary conditions). The second criterion is that it must set the highest spectral coefficient to zero. Adding this type of spectral filtering turns out to be one of the key elements in improving our numerical method enough to make possible the binary black-hole simulations described in Sec. III. The needed spectral filter is applied to each field u(x, t) that is expanded as a sum of Chebyshev polynomials: These are used for the radial-coordinate expansions in the spherical-shell subdomains, and for all three spectralcoordinate expansions in the rectangular-block and the cubed-sphere subdomains. Each spectral coefficient u k in these expansions is filtered according to the expression, after each time step. This filter satisfies the first necessary filter criterion by adding a 2p th -order dissipation term to the evolution equations which vanishes on each subdomain boundary [64]. For the binary black-hole simulations presented here, we use the filter parameters α = 36 and p = 32. The choice α = 36 guarantees that this filter satisfies the second necessary filter criterion by ensuring that the highest spectral coefficient is set to (double precision) zero, while the choice p = 32 insures that this is a very mild filter that only influences the largest few spectral coefficients. This new filter is applied wherever Chebyshev expansions, Eq. (10), are used. We use the same filter we have used in previous binary black-hole simulations in the angular directions in the spherical-shell subdomains: we set to zero the top four ℓcoefficients in the tensor spherical-harmonic represention of each dynamical field [47,48].

C. Adaptive Conforming Grid
One of the important developments that allowed our group (several years ago) to begin performing successful binary black-hole inspiral simulations was the introduction of the dual-coordinate-frame evolution method [47]. This method uses two distinct coordinate frames: One is a non-rotating and asymptotically Cartesian coordinate system used to construct the tensor basis for the components of the various dynamical fields evolved by our code. The second is a coordinate system chosen to follow (approximately) the motions of the black holes. We fix the computational grid to this second coordinate frame, and solve the evolution equations for the "inertial-frame" tensor field components as functions of these "grid-frame" coordinates. The map M connecting these grid frame coordinatesx i to the inertial frame x i , can be written as the composition of more elementary maps: M K represents a "kinematical" map that keeps the centers of the black holes located (approximately) at the centers of the excised holes in our grid structures (cf. Sec. II B), and M S is a "shape-control" map that makes the grid conform (approximately) to the shapes of the black holes. The third major technical development (which makes it possible now for us to perform robust merger simulations) consists of improvements in the choice of the shape-control map, M S , and improvements in the way this map is adapted to the dynamically changing shapes of the black holes. The kinematical map M K can itself be decomposed into more elementary maps: The maps M T , M E , and M R each move the centers of the black holes, but do not significantly distort their shapes: The map M T translates the grid to account for the motion of the center of the system due to linear momentum being exchanged with the near field [54] and being emitted in gravitational radiation. The map M E does a conformal rescaling that keeps the coordinate distance between the centers of the two black holes fixed in the grid frame, as they inspiral in the inertial frame. And M R rotates the frame so that the centers of the two black holes remain on thex axis in the grid frame, as they move along their orbits in the inertial frame. These maps are also used during inspiral simulations in the same way we use them for our merger simulations. So here we focus on the new shape-control map M S , which is one of the critical new developments that made the merger simulations reported in Sec. III possible. The interiors of the black holes are excised in our simulations at the spherical boundaries shown in Figs. 2 and 3. This is possible because, if this excision boundary is chosen wisely, the spacetime inside this boundary can not influence the spacetime region covered by our computational grid. For hyperbolic evolution systems, like the generalized harmonic form of the Einstein equations used in our code [45], boundary conditions must be placed on each incoming characteristic field at each boundary point. Apparent horizons are surfaces that are often used to study black holes, because they can be found numerically in a fairly straightforward way, and because (if they exist at all) they are always located within the true event horizons. If the excision bounaries for our computational domain were placed exactly on the apparent horizons, then all the characteristic fields of the Einstein system would be outflowing (with respect to the computational domain) since apparent horizons are spacelike (or null) hypersurfaces. Boundary conditions would not be needed on any field at such boundaries. Unfortunately we are not able to place excision boundaries precisely on the apparent horizons. So the best that can be done numerically is to place them slightly inside the apparent horizons (i.e., the apparent horizons must remain within the computational domain), if we are to avoid the need for boundary conditions on these excision surfaces.
However, if an excision boundary is placed inside an apparent horizon, the outflow condition is no longer auto-matic or simple: the condition depends on the shape and the location of the excision boundary, its motion with respect to the horizon, and the gauge. In our simulations, the excision boundaries are kept somewhat inside the apparent horizons, so the outflow condition can-and does-fail if we are not careful. One reason (and probably the principal reason) we need to control the shapes of the computational domains through the map M S is to keep pure outflow conditions on all the dynamical fields at the excision boundaries. We do this, as described in more detail below, by requiring that the excision boundaries closely track the shapes and sizes of the horizons.
Another (probably secondary) reason that shape and size control are needed in our numerical simulations is related to finite numerical resolution. If the excision boundaries had different shapes than the horizons, then some points on the boundaries would be located deeper into the black-hole inerior and hence closer to the spacetime singularity. Higher numerical resolution would be needed at these points, and for fixed resolution, the amount of constraint violation and errors in the solution would be largest there. In many situations we find that numerical instabilities at these points cause our simulations to fail. This mode of failure can be eliminated by keeping the shapes and sizes of the boundaries close to those of the horizons.
We have decomposed the map M, which connects the grid-frame coordinatesx i with inertial-frame coordinates x i , into kinematical and shape-control parts: M = M K • M S . It will be convenient to have a name for the intermediate coordinate system whose existence is implied by this split. The map M K connects the inertialframe coordinates x i with coordinatesx i in which the centers of the black holes are at rest (approximately). So it is natural to call these intermediate-frame coordinates, x i , the "rest frame" coordinate system. The shape control map M S connects the grid-frame coordinatesx i with these rest-frame coordinatesx i .
It is useful to express the shape-control map M S as the composition of maps acting on each black hole individually: M S = M S1 • M S2 . Such individual black-hole shape-control maps can be written quite generally in the form (16) where (r A ,θ A ,φ A ) and (r A ,θ A ,φ A ) are, respectively, the grid-frame and rest-frame spherical polar coordinates centered at the (fixed) grid-coordinate location of black hole A = {1, 2}. The functions f A (r A ,θ A ,φ A ) are fixed functions of space in the grid frame. We note that these maps become the identity whenever f A (r A ,θ A ,φ A ) = 0, so making f A (r A ,θ A ,φ A ) vanish asr A → ∞ ensures that the distortion is limited to the neighborhood of each black hole. The parameters λ ℓm A (t) specify the angular structure of the distortion map. They are determined by a feedback control system (discussed in Appendix ) that dynamically adjusts the shape and overall size of the grid frame coordinates relative to the rest-frame coordinates.
The Caltech/Cornell collaboration has used shapecontrol maps of this form, Eqs. (14)- (16), in all of our recent binary black-hole simulations [53,54,56]. In those cases the functions f A (r A ,θ A ,φ A ) were taken to be smooth functions ofr A alone: roughly constant near each black hole and falling off to zero rapidly enough that M S1 and M S2 approximately commute. These functions had large gradients which made it difficult to control the grid distortion near the horizons without introducing additional unwanted distortions elsewhere. As a result, we were never able to achieve very robust shape control using these maps.
One of our major breakthroughs, leading to the successful merger results reported here, is an improvement in our choice of the map functions f A (r A ,θ A ,φ A ). We now use simpler functions f A (r A ,θ A ,φ A ) that have much smaller gradients and that are exactly zero outside (nonintersecting) compact regions surrounding each black hole. This means the new maps M S1 and M S2 commute, exactly. These new maps can be defined everywhere between the black holes because of the new non-overlapping grid structure discussed in Sec. II B, but these improvements are achieved at the expense of smoothness: the new f A (r A ,θ A ,φ A ) are smooth except at subdomain boundaries where they may only be continuous, not differentiable 1 . Fortunately, this lack of smoothness at subdomain boundaries is not a problem for our basic evolution method. In our multidomain code the equations are solved independently on each subdomain, and subdomains communicate only by equating the appropriate characteristic fields at their mutual boundaries. These boundary conditions require no differentiation, so the grid coordinates themselves need not be smooth across these boundaries.
The new f A (r A ,θ A ,φ A ) are chosen to be unity inside a sphere of radiusb A centered on black hole A={1,2}, then decrease linearly with radius along a ray through the center of the black hole, and finally vanish on the surface of a centered cube of size 2ā A . A two-dimensional sketch of this grid-frame domain is shown in Fig. 4, which is merely an abstract version of the type of grid structures we use around each black hole as illustrated in Figs. 1  and 2. This function f A (r A ,θ A ,φ A ) can be expressed analytically as wherē is the value ofr A on the surface of the cube of size 2ā A centered on black hole A={1,2}. The different cases of the maxima which appear in the denominator of Eq. (18) correspond to the different faces of this cube. We note that the choice of f A (r A ,θ A ,φ A ) in Eq. (17) implies that the shape-control map M S , and hence the map between grid-frame and rest-frame coordinates is not smooth. Since the kinematical map M K is smooth, this also imples that the full map M relating grid-frame to inertial-frame coordinates is not smooth either. This non-smoothness does not cause a problem for our basic evolution method, but it implies that calculations that require smoothness must not be performed in the grid frame. For example, the representation of a smooth apparent horizon in grid coordinates will not be smooth. So apparent horizon finding and other calculations that require smoothness are done in the smooth rest-frame or inertial-frame coordinates.
The angular structures of the shape-control maps defined in Eq. (14)- (16) are determined by the parameters λ ℓm A (t). These are chosen dynamically to ensure that the shapes of the excision boundaries match the evolving shapes of the apparent horizons. These maps also control the sizes of the excision boundaries by adjusting λ 00 A , which are chosn to ensure that the excision boundaries remain close to but safely inside the apparent horizons. Figure 5 illustrates the effect of this distortion map on the structure of the grid surrounding a black hole in one of This grid structure has been distorted in relation to the inertial frame by the shape-control maps MS described here. the generic merger simulations described in Sec. III. The choices of suitable target shapes and sizes for these shapecontrol maps, and the feedback control system that we use to implement these choices dynamically in our simulations, are also among the critical technical advances that allow us now to perform robust merger and ringdown simulations. The details of exactly how this is done are somewhat complicated, however, so we defer their discussion to the Appendix .

III. MERGER AND RINGDOWN SIMULATIONS
In this section we present simulations of the merger and ringdown of binary black-hole systems using our Spectral Einstein Code (SpEC). In Section III A we describe the binary black-hole initial data sets that we evolve. In Section III B we discuss our evolution algorithm, pointing out at which stage we employ the various improvements described in Section II and what goes wrong when these improvements are not used. Finally, in Section III C we describe six different binary black-hole merger and ringdown simulations peformed using these new methods. We present snapshots of the shapes and locations of the horizons at various times during the merger and ringdown, and a demonstration of convergence of the constraint violations.

A. Initial Data
Our initial data are in quasi-equilibrium [71,72,73] (see also [74,75]), and are built using the conformal thin sandwich formalism [76,77] with the simplifying choices of conformal flatness and maximal slicing. Quasiequilibrium boundary conditions are imposed on spherical excision boundaries for each black hole, with the lapse boundary condition given by Eq. (33a) of Ref. [73]. The spins of the black holes are determined by boundary conditions on the shift vector at each excision surface [72].
This formalism for constructing initial data also requires the initial radial velocity v r of each black hole toward its partner and an initial orbital angular frequency Ω 0 (chosen to be about the z axis without loss of generality). These parameters determine the orbital eccentricity of the binary, and can be tuned by an iterative method [56,60] to produce data with very low eccentricity. This has been done for most of the initial data sets described here (cases A through D in Table I), but is unnecessary for the purposes of the present paper, since our goal here is simply to document our improved methods for evolving binaries through merger. For cases E and F, v r and Ω 0 are chosen roughly by using post-Newtonian formulae that ignore spins. Table I shows the mass ratio and initial spins for the black hole binary configurations we evolve here. The initial data and first 15 orbits of inspiral for case A are identical to those presented in [49,53]. The initial data and first 9 orbits of inspiral for case C are identical to those presented in [56]. Cases E and F are fully generic, with unequal masses and unaligned spins, and exhibit precession and radiation-reaction recoil.

B. Evolution procedure
In this section we summarize our procedure for evolving black hole binary systems through merger and ringdown, concentrating on the improvements discussed in Section II. This procedure can be divided into several stages, beginning with the early inspiral and extending through ringdown: 1. Evolve early inspiral using "quasi-equilibrium" gauge.
4. Turn on shape control.

5.
Replace remaining inner spherical shells with cubed spheres.
6. Turn on size control just before holes merge.
7. After merger, interpolate to a single-hole grid to run through ringdown.
The above ordering of these stages is not the only possible choice: we have exchanged the order of some of these stages without trouble. For example, for the runs described here we perform stage 3 at the instant we start stage 2. The transitions between these stages are relatively simple, require little or no fine tuning, and can be automated. We now discuss each of these stages in turn. For several of these stages, we will illustrate the effects of new improvements (Section II) for one particular evolution, case F of Table I.

Early inspiral
The gauge is chosen early in the evolution following the procedure of Ref. [49]: the initial data are constructed in quasi-equilibrium, so we choose the initial H a to make the time derivatives of the lapse and shift zero. We attempt to maintain this quasi-equilibrium condition by demanding that H a remains constant in time in the grid frame in the following sense: we require whereH a is a tensor (note that H a is not a tensor) defined so thatH a = H a in the inertial frame. The bars in Eq. (19) refer to the grid frame coordinates. We refer to this as "quasi-equilibrium" gauge, although strictly speaking this gauge maintains quasi-equilibrium only when the spin directions are constant in the grid frame.
For cases A-D, the simulations follow many orbits using this quasi-equilibrium gauge. For cases E and F, in which the spins change direction in the grid frame, this gauge is not appropriate and causes difficulties, so for these cases we transition to harmonic gauge H a = 0 very early in the inspiral.
During inspiral, the map M E uniformly contracts the grid so that the grid-coordinate distance between the centers of the two black holes remains constant. Because this contraction is uniform, the spherical-shell subdomains inside each apparent horizon shrink relative to the horizon. This means that the excision boundary (the inner boundary of the innermost spherical shell) inside each black hole moves further into the strong-field region in the interior, towards the singularity. If this motion continues unchecked, gradients on the innermost spherical shell increase until the solution is no longer resolved. However, the inner boundary of the next-to-innermost spherical shell is also shrinking, and it eventually shrinks enough that the characteristic fields on this boundary all become outflowing (into the black hole); when this happens, the innermost spherical shell no longer influences the exterior solution, so we simply drop that shell from the domain. This shell-dropping process occurs automatically as successive spherical shells shrink relative to the horizons, This process is illustrated in Figure 6, which shows the minimum of the (into the hole) characteristic speeds on spherical grid boundaries around one of the black holes. At any given time, the innermost of these boundaries is an excision boundary. All characteristic speeds at all points on the excision boundary (and hence the minimum) must be positive for the excision algorithm to be well-posed. The innermost shell is dropped whenever the characteristic speeds are positive on the inner boundary of the next shell. For example, surface A 3 in Figure 6 is the excision boundary at t = 24M , but after the characteristic speeds on A 4 become positive (around t = 26M ), the innermost shell is dropped and A 4 becomes the excision boundary. Similarly, A 5 becomes the excision boundary at t = 38, and so on. It is important that the characteristic speed on surface A n in Figure 6 does not become negative before the speed on surface A n+1 becomes positive; otherwise this shelldropping procedure will fail. Shell dropping works well during most of the evolution but it fails in the last stage before merger, which is discussed in Section III B 6. When describing mergers, it is useful to introduce a measure of distance between the black holes. One such measure of distance is indicated on the horizontal axis at the top of Figure 6. This is the spatial (in each time slice) proper separation between the apparent horizon surfaces, obtained by integrating along the rest-framex axis (recall that in the rest frame, the centers of the horizons always remain along this axis). This is not the true proper separation between the horizons because we do not minimize over all possible paths between all possible pairs of points on the two surfaces; however, we use this only as an approximate measure of how far the black holes are from merger. Note that the horizons touch each other when (or possibly before) this proper separation measure falls to zero.
It turns out that the improvements discussed in Section II are unnecessary during early inspiral for cases A-D and are not yet used at this stage. For instance, the domain decomposition uses overlapping cylinders and spherical shells, and not the cubed-sphere subdomains described in Sec. II B. We find no problems while the holes are far enough apart that they remain roughly in equilibrium. However, the binary eventually becomes extremely dynamical, and the black holes become significantly distorted. Unless the algorithm is modified, the simulation fails (typically because of large constraint violations) before the black holes merge.
Typically, it is possible to extend the quasi-equilibrium inspiral stage until the proper separation decreases to about 7M or 8M , but the next stages, Secs. III B 3 and III B 2, can be done sooner if desired. (For example, cases E and F use a non-overlapping cubed-sphere grid from t = 0.) For the simulations presented here, the length of this quasi-equilibrium inspiral stage is variable. For instance, for case A in Table I this stage lasts for 15 orbits, but for case F we transition to damped harmonic gauge and eliminate overlapping subdomains starting at t = 0 because the black holes are initially very close together.

Transition to damped harmonic gauge
When the inspirling black holes reach a proper separation of about 8M , we smoothly transition from quasiequilibrium gauge to damped harmonic gauge at t = t g (g stands for "gauge") by choosing whereH a (t) is the value of H a (t) obtained from the quasiequilibrium gauge condition (19), and σ g is a constant. The last line in Eq. (20) is the damped harmonic gauge condition, Eq. (8). The values of µ L and µ S are set according to Eq. (9), using where t d and σ d are constants. It is necessary to choose σ g and σ d sufficiently large or else the gauge becomes unnecessarily dynamical, possibly resulting in incoming characteristic fields on the excision boundary. For all runs shown here, t d corresponds approximately to a proper separation of 8M . We use the values σ g = 20M, σ d = 50M for the equal mass binaries (cases A, C, and D), σ g = 15M, σ d = 100M for case B, and σ g = 25M, σ d = 50M for case F. For these cases, t g = t d . For case E, the gauge is rolled off to harmonic early in the inspiral, so at t = t d we need only to roll on the damped harmonic gauge, and we use σ d = 40M . The constants in Eqs. (20) and (21) can be chosen quite flexibly: the only constraint is that σ g and σ d must be large enough to avoid significant spurious gauge dynamics, but small enough so that the simulation is using damped harmonic gauge before the black holes approach each other too closely. The σ g and σ d can always be made longer by choosing earlier transition times t g and t d .

Eliminate overlapping subdomains
Although the domain decomposition consisting of overlapping spherical shells, cylinders, and cylindrical shells is quite efficient for the early inspiral, it often suffers from numerical instabilities, particularly as the black holes approach merger. Therefore, we regrid onto the cubedsphere subdomains described in Sec. II B. Use of these non-overlapping subdomains eliminates these remaining instabilities, and effectively increases our resolution. For most of the runs described here, we happen to regrid at the same time t d that we begin the transition to damped harmonic gauge. However, for runs E and F, the entire simulation uses the non-overlapping cube-sphere subdomains. As far as we know, there is no reason that regridding cannot occur at times other than t d or t g .

Shape control
The next stage is to turn on shape control by introducing the map M S defined in Eqs. (14)(15)(16)(17)(18). The parameters λ ℓm A that appear in this map are determined (via a feedback control system) by Eqs. (A.3) and (A.6) for ℓ > 0, and λ 00 A is set to zero. This map deforms the grid so that the boundaries of the spherical shells inside the horizons, including the excision boundaries, are mapped to closely match the horizons' shapes.
We do not turn on shape control until time t g + 2.5σ g , which is when the coefficient multiplying the gauge term H a (t) in Eq. (20) becomes smaller than double-precision roundoff. This ensures that H a remains a smooth function of space in inertial coordinates. If shape control were turned on too early, H a would fail to be everywhere smooth because the quasi-equilibrium-gauge quantityHā is a smooth function of space in the grid-frame coordinates (and its numerical value is constant in time for each grid point), and the deformation is only C 0 across subdomain boundaries. Shape control is necessary for the shell-dropping procedure described in Section III B 1 to remain successful as the holes approach each other. For a shell to be dropped, the entire inner boundary of that shell must remain an outflow surface until the entire inner boundary of the next shell becomes an outflow surface. If these boundaries do not have the same shape as the horizon, then typically some portion of either boundary moves too close or too far from the horizon, and violates this condition. An example is shown in Figure 7. Like Figure 6, Figure 7 shows the minimum characteristic speeds on the boundaries of inner spherical shells for a portion of case F, but Figure 7 shows two separate evolutions: one with shape control turned on at t = 62M (solid curves) and one with no shape control (dashed curves). Without shape control, the minimum characteristic speed on surface A 7 stops growing around t = 100M and becomes negative soon thereafter. Because

Eliminate inner spherical shells
Spherical-shell subdomains are extremely efficient when the solution is nearly spherically symmetric, because these subdomains use Y lm basis functions, which are well-suited for nearly spherical functions. However, once the holes become sufficiently distorted, it is difficult to resolve the region near each horizon using spherical shells. At this point in the evolution, we replace the inner spherical shells with cubed-sphere shells. Figure 8 shows the constraints as a function of time for two evolutions of case F: one with inner spherical shells (dashed curves), and another in which the inner spherical shells were replaced with cubed-sphere shells at t = 80M (solid curves). The quantity plotted is ||C GH ||, the L 2 norm over all the constraint fields of our first-order generalized harmonic system, normalized by the L 2 norm over the spatial gradients of all dynamical fields (see Eq. (71) of Ref. [45]). The L 2 norms are taken over the portion of the computational volume that lies outside the apparent horizons. For the spherical-shell case in Figure 8, the constraints in the spherical shells grow large enough that the simulation terminates at t = 129.5M , when the proper separation has fallen to 0.5M . Note that these simulations typically proceed through merger even if we do not replace spherical shells with cubed-sphere shells. So while this replacement step is not strictly necessary, eliminating spherical shells improves the accuracy of the simulation (as measured by the constraint quantity ||C GH ||) by almost an order of magnitude.

Size control
Eventually, even using shape control, shell dropping typically fails before the horizons merge: the outflow condition on the inner boundary of the inner shell fails before the outflow condition on the inner boundary of the next shell becomes satisfied. This problem occurs because the map M E shrinks the entire grid-including regions inside each hole-faster and faster as the holes approach each other. As soon as the velocity of excision boundary towards the center of the hole becomes large enough, the boundary becomes timelike, and the excision algorithm fails. To remedy this, we turn off shell dropping and we turn on size control. Size control slows down the infall of the excision surface (by pulling the excision sphere towards the horizon), thus keeping the characteristic speeds from changing sign. Size control is implemented by changing the map parameters λ 00 A (t), which previously were set to zero, to values given by Eqs. (A.3) and (A.10). This causes the size of each excision boundary to be driven towards some fraction η of the size of the appropriate horizon. For the equal-mass runs A, C, and D we use η = 0.8 for both black holes. For the unequal-mass run B and the generic cases E and F we use η = 0.9 and η = 0.7 for the smaller and larger black holes, respectively. These values do not require finetuning; for instance, case E still runs through merger and ringdown if we change η to 0.9 for both black holes. The main criterion for choosing this parameter is that a grid sphere of radiusb A (see Figure 4) should not be mapped outside the cube of side 2ā A . If this occurs, the map M S becomes singular and the run fails. Figure 9 shows the minimum characteristic speed at particular domain boundaries for two evolutions of case F: one with size control and one without. Without size control, the outflow condition on the excision boundary fails at proper separation 1.67M and time t = 127.24M . With size control, the evolution proceeds to proper separation 0.015M , time t = 130.17, well past the formation of a common apparent horizon, which occurs at proper separation 1.4M (time t = 128.23M ).

Ringdown
At some time t = t m shortly after a common apparent horizon forms, we define a new grid, composed only of spherical shells. This grid has only a single excision boundary inside the common horizon. We also define a new map M ringdown :x i → x i that can be written as a composition of simpler maps: This is similar to the ringdown maps described in Refs. [53,56] (which had no translation) and Ref. [54] ( As described in Ref. [53], at t = t m the parameters λ ℓm A are chosen so that the common apparent horizon is stationary and spherical in the grid frame. Similarly, the map M Tf is chosen at t = t m so that the center of the horizon is located at the grid-frame origin and is stationary in the grid frame. Once we have defined the maps at t = t m , we then interpolate all variables from the old grid to the new grid. Note that the grid frame changes discontinuously in time at t = t m . Because of this, we take care to properly use both the new map M ringdown and the old map M in doing this interpolation, so that the inertial frame and quantities defined in that frame remain smooth. In particular, the gauge source function H a is still determined by Eq. (20) during ringdown, and remains smooth at t = t m .
For t > t m , the parameters λ ℓm A (t) are determined by a feedback control system that keeps the common apparent horizon stationary in the grid frame. Likewise, the map M Tf keeps the common apparent horizon centered at the origin in the grid frame.

C. Results
In this section we present some results from the simulations listed in Table I. We first show snapshots of the inertial-coordinate shapes of the apparent horizons during merger and after ringdown. Figures 10, 11, 12, and 13 show cross-sections of the horizons in the orbital plane for cases A-D; in these cases, the orbital plane is well-defined and is constant in time. For all cases, we show the apparent horizons of the individual black holes and the common apparent horizon at the time when the common horizon is first detected (solid curves), and at the time t m when we transition to a new grid that has a single excision boundary (dotted curves). We also show the apparent horizon of the final remnant black hole after it has reached equilibrium (dashed curves). Figures 14  and 15 show cross-sections of the horizons for cases E and F, in which the orbital plane precesses. For these cases we show cross-sections in the coordinate plane (defined  Figure 10, except for an anti-aligned-spin equal-mass binary black hole merger (case C).
using the flat metric) that is perpendicular to the instantaneous orbital angular velocity at the time the common horizon is first detected, and that passes through the coordinate center of the common horizon at this time. We show cross-sections at three different times: the time of horizon formation, the time t m , and a time 300M after merger, all with respect to the same plane. The remnant black hole in cases E and F have nonzero linear momentum both because of radiation reaction and because of some nonzero linear momentum present in the initial data; this is why the centers of the horizons shown in Figures 14 and 15 change with time.
For case F, we also show the constraint norm as a func-  Figure 10, except for an aligned-spin equal-mass binary black hole merger (case D).
FIG. 14: Three snapshots of the apparant horizons of a generic binary black-hole merger (case E). The plane of the cross-sections is a coordinate plane perpendicular to the instantenous orbital axis at the time the common horizon is first detected. The additional degree of freedom is fixed by having this plane go through the coordinate-center of the shape of the common horizon at its first detection. Solid curves are the cross section when a common horizon is first detected; dotted curves represent the time of transition between the binary merger and the single hole ringdown evolutions; dashed curve shows the horizon at a time well after merger.
tion of time in Figure 16. We plot the same quantity ||C GH || as shown in Figure 8. This quantity is shown for four numerical resolutions, and is convergent at all times (although the convergence rate is smaller near merger when the solution is most dynamical). The constraints are largest at t = t m , when we transition to a grid with  Figure 14, except for the generic binary black-hole merger case F in Table I. a single excision boundary. Just after t = t m the constraints decrease discontinuously by a small amount because part of the computational domain has been newly excised. The approximate number of grid points for the four resolutions is 79 3 , 87 3 , 95 3 , 103 3 just before merger and 44 3 , 51 3 , 57 3 , 63 3 during ringdown.

Acknowledgments
We thank Jan Hesthaven for helpful discussions. We are grateful to Luisa Buchman, Larry Kidder, and Harald Pfeiffer for the improved kinematical coordinate maps and associated control systems that handle inspirals of precessing and unequal-mass binaries, and we thank Harald Pfeiffer for help preparing the initial data. We also thank Luisa Buchman, Tony Chu, Geoffrey Lovelace, and Harald Pfeiffer for producing some of the BBH inspiral simulations that we here extend through merger and ringdown. We acknowledge use of the Spectral Einstein Code (SpEC). This work was supported in part by grants from the Sherman Fairchild Foundation, by NSF grants PHY-0601459, PHY-0652995, and by NASA grant NNX09AF97G.

APPENDIX: CONTROLLING THE MAPS
The purpose of the shape-control maps M SA , with A={1,2}, is to distort the grid structures around the black holes so the excision boundaries of the computational domain are mapped to surfaces lying just inside and having the same basic shapes as the apparent horizons. Choosing the right maps is equivalent to choosing the right target values for the parameters λ ℓm A (t) that define these maps, cf. Eq. (16). This Appendix describes in some detail how these parameters are chosen, and what target surfaces are used to fix these maps in the merger simulations described in Sec. III.
The grid structures used in our merger simulations, cf. Sec. II B, have excision boundaries that are grid-frame coordinate spheres. The target surfaces S T A to which we want to map these grid-frame spheres can be written in the formr where the expansion coefficients Λ ℓm A (t) define the target surface and (r A ,θ A ,φ A ) are rest-frame coordinates.
Letr target A denote the radii of the grid-frame coordinate spheres that are to be mapped onto the target surfaces S T A . The rest-frame representations of these coordinate spheres, using Eqs. (14)(15)(16), are: Note that the functions f A that appear in Eq. (17) are set equal to unity here because we always choose the target grid-spheresr target A to be smaller than the outer radii of the largest spherical subdomain layers:b a ≥r target A . It follows that the target spheres will be mapped to the shapes of S These conditions can not be imposed directly for a number of reasons, but they can be enforced approximately using a feedback control system. The control system used in the merger calculations described in Sec. III is the same one used to perform our earlier binary black-hole simulations [47], so we will not describe it in detail here.
To complete the specification of the shape-control maps we must choose the target grid-spheresr target A , and the target surface parameters Λ ℓm A . From Eq. (A.3) it follows that the choice ofr target A only effects the ℓ = 0 size control part of the distortion map. So until size control is activated late in our runs, there is no need to specify whatr target A actually is. We think of it as being the radius of the suitably scaled apparent horizon, but this fact does not influence the shape control part of the map in any way. Once size control is activated late in the run, then we requirer target Next consider the choice of target surfaces S T A , starting with the ℓ > 0 contributions that control the shape but not the overall scaling of the map. The idea is to choose the target surfaces S T A to be similar in shape to the apparent horizons H A . The H A can be represented as smooth surfaces in the rest-frame coordinate system: To keep the shapes of the target surfaces similar to the shapes of the apparent horizons, the target surface parameters Λ ℓm A should be made proportional to S ℓm A . We find it is appropriate to scale these coefficients by a factor, G(R A ), which depends on the average radius of the apparent horizon,R A = S 00 A Y 00 : The larger the apparent horizon radius in relation to the desired excision radius, the smaller this scaling factor must be to maintain an appropriate shape for the excision boundry. In practice, we find the scale factor works quite well for p = 2. This scaling factor is near unity whenR A <ā A and decreases like 1/R A for larger values ofR A . The tanh x function is introduced here because it is linear for small values of |x|, and approaches 1 for x ≫ 1. The p-dependence maintains these asymptotic forms, but allows the transition between them to occur more quickly than the p = 1 case. We have found no need to introduce the ℓ = 0 size control parts these maps, until the final plunge phase just before the black holes merge. So during most of our merger simulations we simply set λ 00 A (t) = 0. This allows the apparent horizons to grow relative to the grid coordinates as the map M E contracts the rest-frame relative to the inertial coordinates. This growth in the apparent horizons in the grid frame allows us to drop a number of spherical subdomain layers as the evolution progresses, and this helps remove unwanted constraint violations from the computational domain.
At some point however, we may not be able to drop additional spherical subdomain layers: either because we run out of layers, or because the inner boundary experiences incoming characteristic fields when the apparent horizon expands too quickly during the final plunge. When either of these conditions occurs we turn on the ℓ = 0 part of the shape-control maps, to keep the radius of the excision boundaryr ex A close to the size of the apparent horizon. It will be useful to define We choose η so that when ∆r A < 0 there is no need to impose size control, and use η in the range 0.7 ≤ η ≤ 0.9 for the merger simulations reported in Sec. III. When ∆r A becomes positive we want to turn on size control, but we want to do it in a fairly smooth and continuous way. This is done by defining the function: which vanishes for ∆r A ≤ 0, and asymptotically approaches P (∆r A ) → 1 for ∆r A ≫R A /40. We find that an appropriate value for the overall scale of the target surfaces S T A , i.e. their ℓ = 0 components, to be Once size control is activated, the target radius is given by Eq. (A.4):r target A =r ex A . This choice implies that the average rest-frame radius of the excision boundary, R ex A =r ex A − λ 00 A Y 00 , approaches ηR A when ∆r A R A /40.
To summarize, the appropriate shape and size control of the inner grid structures in our binary blackhole merger simulations are deterined by the target shape-control surfaces S T A , whose definitions are given in Eqs. (A.6) and (A.10). The shape-control map needed to distort the grid coordinates to the shape of S T A , is given in Eq. (A.3). This mapping is enforced with an appropriate feedback control system.