{'fraction_non_alphanumeric': 0.03211906141367323, 'fraction_numerical': 0.010175260718424102, 'mean_word_length': 4.91117294520548, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Status of Black-Hole Binary Merger Simulations with Numerical Relativity

The advent of long-term stability in numerical relativity has yielded a windfall of answers to long-standing questions regarding the dynamics of space-time, matter, and electromagnetic fields in the strong-field regime of black-hole binary mergers. In this review, we will briefly summarize the methodology currently applied to these problems, emphasizing the most recent advancements. We will discuss recent results of astrophysical relevance, and present some novel interpretation. Though we primarily present a review, we also present a simple analytical model for the time-dependent Poynting flux from two orbiting black holes immersed in a magnetic field, which compares favorably with recent numerical results. Finally, we will discuss recent advancements in our theoretical understanding of merger dynamics and gravitational waveforms that have resulted from interpreting the ever-growing body of numerical relativity results.


I. INTRODUCTION
In the last five years, the field of gravitational wave (GW) astronomy has witnessed a remarkable convergence of experimental and theoretical achievement. The current generation of GW interferometers, particularly the Laser Interferometer Gravitational Wave Observatory (LIGO) and VIRGO detectors, has positioned this community on the cusp of discovery, by achieving their design sensitivities and making the first direct detection of GWs a real possibility. Meanwhile, the GW signature from the merger of a black-hole binary (BHB), expected to be the most common source for these and future observatories, went from being a complete unknown to a reasonably well-understood, surprisingly smooth transition between the earlier inspiral and the final ringdown, thanks to tremendous advancements in the field of numerical relativity (NR). In this article, we will limit our scope to the advancements that have occurred in just the past year, with only a minimal amount of background information, in order to provide a snapshot of the most current research in the field. For a review of progress in the preceding year, we refer the reader to [1], and for a broader review of NR achievements, please see [2,3]. We consider it an indicator of the healthy progress of the field that, in addition to the general review articles typical of any field, in recent years the field of NR has progressed rapidly enough to warrant an annual review of the past year's progress.
In Sec. II, we give a brief review of the equations being solved by numerical relativists, and the methods that are most typically applied to solve them. In Sec. III, we discuss the most recent achievements in NR methodology, and the simulations of numerically challenging systems of astrophysical interest that have been facilitated. * Electronic address: sean@astro.columbia.edu In Sec. IV, we briefly mention novel theoretical developments resulting from the interpretation of numerical results. We summarize and conclude in Sec. V.

II. STANDARD METHODOLOGY
The community of numerical relativists attempts to solve Einstein's equations in the strong-field regime numerically (hence the name).
During the early, approximately-adiabatic inspiral, we can treat the system perturbatively through expansions in v/c or other small parameters related to it. After the merger, we can again treat the system perturbatively, as the quasinormal ringing of small deviations from a Kerr background. However, during the merger, unless the masses are disparate and q = m 1 /m 2 provides a tractably small expansion parameter, we can no longer treat the system perturbatively, and are forced to solve the full equations, given by where R µν , R, g µν , and T µν are the Ricci tensor, its trace, and the metric and stress-energy tensors, respectively. In the case of BHBs in vacuum, T µν vanishes. The resulting equation is generally evolved by splitting up space-time into a sequence of three-dimensional hypersurfaces. The choice of how to make this split is not unique, as different evolution equations with different variables can be used, as well as different coordinates. The coordinate freedom is fixed through the choice of dividing space-time into hypersurfaces. The lapse, α, determines the coordinate time interval between one hypersurface and the next. The shift, β i , determines the motion of the three spatial coordinates within each hypersurface. The splitting of space-time results in a set of evolution equations, plus constraints that need to be satisfied. To the extent that the constraints vanish, all evo-lution equations must be analytically equivalent. However, different choices will have different numerical properties. Indeed, the four decades that passed between the first attempts to perform NR simulations and the recent breakthroughs were spent, in large part, searching for formulations that would be numerically stable.
At present, there are primarily three analytic forms for Einstein's equations that are actively being applied to the BHB problem: the BSSN formulation [4][5][6] and two different Generalized Harmonic formulations [7][8][9]. Most typically, the equations are spatially discretized with finite differences (with the exception of SpEC [10], which uses pseudo-spectral methods). The state-of-theart codes employ eighth-order spatial differencing, but most codes have a mixture of orders for different parts of the code. All implementations of which we are aware discretize in time using Runge-Kutta methods, usually to fourth-order accuracy in the time step. Given the disparate scales involved in the problem (r > ∼ 100M for the wave zone, ∆ < ∼ M/100 for the near-field), a range of resolutions within each simulation is required, and this is achieved in a number of different ways by different groups. There are many codes being actively used and developed [11], [12], [8], [13], [14], [10], [15], [16], [17], and comparisons have been published that verify the consistency of results among different groups [18,19].
In the presence of matter and/or electromagnetic (EM) fields, the equations governing the evolution of spacetime must be coupled to an appropriate stress-energy tensor. It is noteworthy that matter is not needed for T µν to be nonvanishing. In the case of the Einstein-Maxwell system (see e. g. [20]), the sources for EM fields are not included in the stress-energy tensor. The inclusion of these sources, the fields themselves, and the equations of hydrodynamics, all coupled to Einstein's equations, constitute the equations of general relativistic magneto-hydrodynamics (GRMHD), while the coupling of the hydrodynamics and Einstein's equations, without EM fields, make up the general relativistic hydrodynamics (GRHD) equations. The recent advent of GR(M)HD codes has yielded some of the most interesting science in NR, and perhaps in all of theoretical astrophysics, in the past year, and will be discussed in Sec. III D.

III. NOVEL METHODOLOGY AND ASTROPHYSICS
In this section, we will discuss novel developments in NR methodology, and the novel astrophysical results that those developments have facilitated.

A. "Extreme" mass ratios
Whereas the typical definition of extreme mass ratio systems is q ∼ 10 −6 , such systems are of little interest to numerical relativists, as they can be treated well with perturbative methods, and they are, at present, completely impossible to simulate non-perturbatively. However, for more intermediate mass ratios, q > − 1/100, it is now possible to perform short simulations which include 1-2 orbits of the late inspiral, the merger, and the ringdown. This is an impressive achievement, because the resolution demands are set by the smaller black hole, so that a black hole 100 times smaller than its companion requires 100 times finer resolution. The trailblazers in the endeavor to simulate these systems have been the RIT group, using their LazEv code [21]. Their first effort involved combining a full numerical evolution of the punctures with a perturbative treatment of the radiation [22]. Later, they validated this approach by comparing with fully numerical simulations for a q = 1/10 system [23]. Most recently, they have achieved the aforementioned 1/100 mass ratio through a combination of improved gauge conditions, an improved mesh refinement scheme, and an improved allotment of cpu cycles [24].

B. Toward Maximal Kerr Initial Data
The Caltech-Cornell-CITA collaboration has overcome the obstacles of generically stable black-hole binary evolutions, and have demonstrated the ability to stably evolve systems with moderate mass ratios q < − 1/2 [25] and extremely high spins a = S/M 2 < − 0.95 [26], including merger and ringdown, using their SpEC code [10]. To facilitate the simulation of stable near-extremal spins, they had to employ novel methodology in order to exceed the bound of a = S/M 2 < ∼ 0.93 that is an intrinsic characteristic of Bowen-York initial data [27]. More specifically, while the spin in Bowen-York initial data can be set arbitrarily high, it will quickly relax to the aforementioned bound, so that larger spins cannot be evolved with that approach. To this end, the authors of [26] applied a novel method [28] wherein the extended-conformal-thin-sandwich equations [29] were combined with a conformally-curved metric, constructed through the weighted superposition of the metrics for two boosted, spinning black holes in Kerr-Schild coordinates. Using this approach facilitated the simulation of a pair of equal-mass black holes with spins of a = 0.95 anti-aligned with the orbital angular momentum (thereby minimizing the total system angular momentum) [26]. We note that the Illinois group also developed a novel approach for near-extremal initial data [30], but it has not yet been demonstrated in a binary evolution.

C. Nonlinear memory
An interesting consequence of the nonlinear nature of GW propagation is the emission of GWs generated by GWs. This effect has long been established analytically, and is often referred to as the "Christodoulou memory" [31]. However, this effect is exceedingly small, and is therefore difficult to simulate, as systematic error in the simulations must be smaller than the effect. The GW memory, though small, could conceivably be detectable by the space-based Laser Interferometer Space Antenna (LISA) [32]. The authors of [33] were able to accurately calculate the memory effect from the final merger and ringdown, by using novel numerical methodology implemented in their Llama code, which includes multipatch capabilities [15,34] and Cauchy characteristic extraction (CCE) [35,36]. The multipatch code facilitates multiple coordinate patches in different regions, allowing traditional Cartesian grids around the near zone, with a smooth transition to a six-patch polar grid in the wave zone. Since the waves do not gain angular structure in the wavezone, this allows one to increase the grid size by only adding radial elements, thereby greatly decreasing the cost of large grids, and permitting much more accuracy in the region where the waves propagate. The implementation of CCE represents the first fully gauge invariant GW calculation, wherein the GWs are formally measured at infinity, rather than extrapolating to infinity from a sequence of finite extraction radii (though we note that progress has also been made in the last year on compactification, including infinity explicitly in the computational domain [37]). Using Llama, the authors of [33] were able to accurately resolve the nonlinear memory for the first time. We note that, in [32,33], the h 20 component of the GW strain (which contains the memory) appears to be a significant fraction of the fundamental quadrupole modes, so that some readers might be confused why the signal-to-noise ratio (SNR) of the memory contribution is roughly two orders of magnitude below the SNR of the total signal (see Fig. 2 of [32]). In this context, one needs to bear in mind that constant offsets in strain are gauge-dependent and undetectable. Therefore, only the derivative of the memory contribution to the GW strain contributes to the SNR [38].

D. GR(M)HD
Perhaps the most exciting developments have been the coming-of-age of GR(M)HD codes. There are currently four codes capable of GRMHD with dynamical spacetime: WhiskyMHD [16], SACRA [17], the LSU-LIU-BYU-PI collaboration code [39], and the Illinois group code [40]. A number of other groups have also developed, or are currently developing, working GRHD codes. As many results of astrophysical interest do not require including electromagnetism (EM), these codes are also an exciting development. In addition to implementing several stable evolution schemes that satisfy the constraints of Einstein's equations along with other consistency requirements (such as preserving the divergencefree character of the magnetic field for GRMHD), most of the aforementioned codes have either incorporated, or are currently developing, more relevant physics for problems involving white dwarfs, neutron stars, and any gaseous environment, such as photon and neutrino transport, and stellar equations-of-state (EOS), ranging from simple polytropes to microphysical EOSs. While most of the work in this field has so far centered on mergers of binaries that include a neutron star or white dwarf, and is therefore beyond the scope of this review, there have been a few results from GR(M)HD codes applied to BHBs in gaseous environments.
The Illinois group estimated the accretion luminosity from a BHB immersed is a gaseous environment using GRHD, and compared the result with the luminosity of a single black hole of the same total mass, immersed in the same environment [41]. By estimating the contributions from bremsstrahlung and synchrotron emission, they predicted an enhancement of luminosity by three orders of magnitude at BHB merger, compared to the single-black-hole case. If true, this result would have a significant impact on the search for EM counterparts of GW signals from BHBs.
In addition to accretion luminosity from a surrounding disk, significant luminosity may accompany BHB mergers in the form of jets, particularly in the presence of an EM field. Much progress has been made in studying this possibility. In [20], which built upon earlier studies [39,42], the Einstein-Maxwell equations were evolved for the case of a BHB immersed in an EM field, with different BHB spin configurations. By calculating the Poynting flux, they showed that the direct EM emission was 13 orders of magnitude smaller than the GW emission. Furthermore, as the EM frequency evolution tracked that for GWs, the direct EM signature would occur in the mHz range for supermassive BHBs, and would therefore be at far too low a frequency to be detectable.
However, by combining the evolution of the Einstein-Maxwell equations with a tenuous plasma, which is evolved using the force-free approximation [43,44] wherein the inertia of the plasma is neglected, the authors of [45,46] were able to calculate the synchrotron radiation resulting from the acceleration of the plasma due to the EM Poynting flux previously discussed, upconverting that energy to GHz, and thereby making it detectable by X-ray observatories. The radiation transitions from an m = 2 multipolar structure to an m = 0 structure, which has implications for determining whether the EM emission is in fact of the same origin as a given GW observation, since there will likely be many variable EM sources within the error ellipse of any GW observation. However, the most interesting aspect of [45], from a theoretical standpoint, is the fact that the EM emission is generated by a BHB consisting of two nonspinning black holes. Since the emission from a single black hole is most typically associated with the Blandford-Znajek (BZ) mechanism [43], which requires that the black hole be spinning in order to operate, this new result may be surprising at first. The flux in the BZ process originates from the fact that a spinning black hole will twist magnetic field lines, which are assumed to be anchored to a distant accretion disk. This twisting gives rise to the EM Poynting flux, with the torquing of the magnetic field lines causing a decrease in the black hole's spin.
In the case of a BHB, it has now been shown that there is a significant contribution to EM emission that is a generic feature, depending only on the presence of an orbiting binary and an EM field. Indeed, given the observation of protostellar jets and their apparent similarities to their extragalactic counterparts [47], it is possible that the conditions for jet formation may be far more generic than was previously thought, with the presence of plasma and an appropriate magnetic field strength and geometry being the only clear prerequisites.
In the context of BHBs, the existence of jets, independent of spin, may not be as surprising when thought of within the framework of the membrane paradigm [48], wherein the black hole event horizon is reinterpreted as a two-dimensional viscous membrane, and the standard Maxwell equations in three dimensions govern the dynamics of EM fields evolving along a sequence of hypersurfaces. It is well understood that, from this viewpoint, the BZ effect can be calculated using simple circuit equations. In this picture, a magnetic field threading the black hole and frozen in to a distant accretion disk serves the role of electrical wires, with charges spiraling along them. The disk and the horizon membrane serve as resistors, and the potential results from the twisting of the magnetic field lines by the spin of the black hole.
The authors of [45,46] suggested that, in the case of a BHB, the Faraday induction results from the orbital motion of the two black holes through the fixed magnetic field, and could also be understood through application of the membrane paradigm. We carry out the calculation here, which does indeed bear out the scaling behavior suggested in [45,46]. Assuming that the plasma is lowdensity and the force-free approximation applies, that the orbit is in the {ρ, φ} plane with orbital velocity v (though any relative velocity between the black hole and the magnetic field will have the same effect), and that the background magnetic field is given by B = Bẑ (in cylindrical coordinates {ρ, φ, z}), the induced potential from one orbiting hole is then given by (2) where α ≡ 1 − 2M/r is the lapse between hypersurfaces of constant Schwarzschild time, M the combined rest mass of both black holes, and r H = M is the horizon radius of either hole. The contour integration for employing Faraday's law in Eq. 2 is somewhat subtle. We place one side along one hemisphere of the stretched horizon, and assume the opposing side resides far away, where the B-field is weak, so that its contribution to the integral vanishes. We note that this region is assumed to be far away from the computational domain in [45,46], where the B-field is uniform. The other sides of the contour are perpendicular to v × B where it remains uniform, so that those contributions also vanish, and only the integral along the horizon remains.
To leading order, the orbital velocity is given by v = M/r, so that Eq. 2 can be re-expressed as Using the remarkable result that the membrane has an effective resistance, R H , of 377 Ω (or 4π in geometrized units) [49], and bearing in mind that each hole emits in both hemispheres, we can solve for the total Poynting flux from both orbiting holes: We note that it is often assumed that the astrophysical load also contributes a resistance, R L ≈ R H [46,48], but as this is rather ad hoc and depends on poorly understood astrophysics, and there is no apparent way that a distant astrophysical load could manifest itself in the simulations of [45,46], we neglect this contribution. Eq. 4 predicts a peak luminosity at r = 4M , which we can express in cgs units with appropriate scaling relationships as where a is the dimensionless spin parameter of the Kerr black hole for the single black hole BZ process (0 < − a < − 1), and L BZ is understood to be the luminosity from the BZ process for a single black hole of spin a and mass M (to approximate the mass of the merged remnant). We emphasize that this estimate is only valid prior to merger, and indeed predicts a vanishing luminosity for r = 2M , so that the merger and its accompanying burst of luminosity seen in [45,46] are not included. However, if we combine Eq. 4 with the leading order relationship for r(t), we see excellent agreement between our simple model and the data from Fig. 4 in [45] (see Fig. 1).
We also note that this model differs significantly from that suggested in [50]. Using Eq. 7 in [50] to find the prediction for L max in that model, we find that the prediction is approximately four orders of magnitude larger than the numerical result found in [45]. While neither model should be expected to perform well in predicting the true peak luminosity, which occurs in a highly dynamical and strong-field regime, the suggestion in [50] that L ∝ M 3 B 3 yields vastly different predictions from Eq. 4, and disagrees with the scaling suggested in [45].

E. State of the art
It is useful to summarize the boundaries of parameter space that have currently been explored. While some-   what arbitrary, a summary of the extreme cases of different quantities that have been simulated is useful and serves as a snapshot of the current status of the field. We include a set of physical quantities from NR simulations, the superlative result for that quantity, and the reference wherein that result was found, in Table I.
The fact that a number of the extremes have not occurred within the past year is perhaps indicative that the tide of NR breakthroughs that began in 2005 has begun to ebb. However, as the pace of 2005-2009 could not be maintained indefinitely, it is our hope that the slowdown is simply the start of a period of more modest, but sustainable, growth in the field. We emphasize that the list of superlatives is restricted to results that have been achieved in full BHB simulations, as opposed to results from single-black-hole simulations, or extrapolations of a set of BHB simulations to some extreme. We also emphasize that we are summarizing the extremes in physical systems that have been simulated, and not including extremes in the numerics (such as smallest phase error, etc.), as any set of such results is arbitrary, and difficult to summarize in an unambiguous way.

IV. THEORETICAL ADVANCES
Apart from achievements in simulating new systems, a burgeoning subfield of NR is the ongoing effort to draw general theoretical insights from the ever-growing body of available NR simulations. Independent research projects at Caltech [56,57] and Chicago [58] have endeavored to gain analytical insight into the nonlinear dynamics of momentum transfer in merging BHBs, particularly with respect to the bobbing and ultimate recoil. Substantial work has been done to condense the body of simulations into a simple fitting formula, based on either the final state of the binary just prior to merger [59] or on the initial conditions at wider separations [60]. More general efforts to understand the complete final state of the merged remnant, including the final mass and spin, have been successful in predicting subsequent numerical results [61,62].
Progress has also been made in moving beyond purely phenomenological fits of merger waveforms. Previously, progress was made by fitting to natural extensions of the post-Newtonian [63] or effective-one-body formalisms [64,65]. More recently, attempts have been made to apply our improved theoretical understanding of the dominant dynamics during merger in order to formulate more physically-motivated merger waveforms. These attempts are an effort to move beyond simply bridging the gap between the well-understood inspiral and the wellunderstood ringdown, in order to encapsulate the key features of generic mergers, in the hopes that the results will prove to be predictive for as-yet-unsimulated (or currently unsimulatable) systems. In [66,67], the close-limit approximation was applied in an effort to predict the recoil of comparable-mass BHBs. This research program has a long history, but the novelty in these investigations is the quality and analyticity of the post-Newtonian initial data, and the existence of full NR results, which largely validate the effectiveness of this procedure. Similarly, in [68], a hybrid of post-Newtonian and black hole perturbation theory techniques was used to study the radiated energy and momentum and the gravitational waveforms for head-on collisions. The results compared favorably with full NR simulations.
In [69,70], the "implicit rotating source" model for merger waveforms was developed, and its results were rigorously compared with full numerical simulations for nonspinning systems with mass ratios q < − 1/6. This model was built on the observation that merger waveforms from numerical simulations behave as though they were generated by a rigid rotator. Specifically, the frequency evolution of all ℓ = m modes, when weighted by m, are equal across all modes and mass ratios, consistent with multipoles of a single object rotating with a welldefined frequency. The angular momentum at merger is also proportional to the frequency, as one would expect for a rotator with a well-defined moment of inertia. These and other indicators led the authors of [69] to formulate a physically-motivated merger waveform, which naturally ties to the final ringdown phase of the waveform, whereas most other methods tie to the late inspiral waveform, at a point where it may have already accumulated a large degree of error. Combined with the aforementioned models for the final mass and spin of the merged remnant across parameter space, this model may provide a novel method for predicting the merger waveform for generic systems, with a minimum of fine-tuning.
The commonality of all the approaches to gaining analytical insight into merger dynamics is the observation that the merger, though happening in the strong-field, nonlinear regime, seems to be well-described by tools developed primarily for linear perturbative analyses in the weak-field regime. This observation is one of the most fundamental results that has been made possible by numerical relativity.

V. CONCLUSIONS
We have presented a brief overview of the current state of black-hole binary simulations in numerical relativity. While progress has slowed somewhat since the "gold rush" period following the breakthroughs in 2005-2006, there have still been a number of novel results of significant astrophysical and theoretical interest. In the past year, progress has been made in simulating longer waveforms overall, as well as systems with more extreme mass ratios and larger spins. Advances in methodology have significantly enhanced the accuracy and efficiency achievable by state-of-the-art codes. Progress has also been made in using the body of available numerical simulations to inform a greater theoretical understanding of the strong-field behavior of black-hole binaries. The inclusion of matter and electromagnetic fields in blackhole binary simulations has led to unexpected discoveries, which will have significant astrophysical implications in the years to come. The field of numerical relativity has now entered a period of more modest, but hopefully sustainable growth, with a substantial amount of discovery space that remains to be explored.