Binary Neutron Star Mergers

Abstract

We review the current status of studies of the coalescence of binary neutron star systems. We begin with a discussion of the formation channels of merging binaries and we discuss the most recent theoretical predictions for merger rates. Next, we turn to the quasi-equilibrium formalisms that are used to study binaries prior to the merger phase and to generate initial data for fully dynamical simulations. The quasi-equilibrium approximation has played a key role in developing our understanding of the physics of binary coalescence and, in particular, of the orbital instability processes that can drive binaries to merger at the end of their lifetimes. We then turn to the numerical techniques used in dynamical simulations, including relativistic formalisms, (magneto-)hydrodynamics, gravitational-wave extraction techniques, and nuclear microphysics treatments. This is followed by a summary of the simulations performed across the field to date, including the most recent results from both fully relativistic and microphysically detailed simulations. Finally, we discuss the likely directions for the field as we transition from the first to the second generation of gravitational-wave interferometers and while supercomputers reach the petascale frontier.

1 Introduction

Binaries composed of neutron stars (NSs) and black holes (BHs) have long been of interest to astrophysicists. They provide many important constraints for models of massive star evolution and compact object formation, and are among the leading potential sources for detection by gravitational-wave (GW) observatories. While it remains uncertain whether mergers of compact binaries are an important contributor to the production of r-process elements, they are now thought to be the leading candidate to explain short-duration, hard-spectrum gamma-ray bursts (often abbreviated to “short-hard” GRBs, or merely SGRBs).

The first neutron star-neutron star (NS-NS) binary to be observed was PSR B1913+16, in which a radio pulsar was found to be in close orbit around another NS [135]. In the decades since its discovery, the decay of the orbit of PSR B1913+16 at exactly the rate predicted by Einstein’s general theory of relativity (see, e.g., [306, 325]) has provided strong indirect evidence that gravitational radiation exists and is indeed correctly described by general relativity (GR). This measurement led to the 1993 Nobel Prize in physics for Hulse and Taylor.

According to the lowest-order dissipative contribution from GR, which arises at the 2.5PN level (post-Newtonian; where the digit indicates the expansion order in [υ/c]² in the Taylor expansion term), and assuming that both NSs may be approximated as point masses, a circular binary orbit decays at a rate da/dt = −a/τ_GW where a is the binary separation and the gravitational radiation merger timescale τ_GW is given by

$$\begin{array}{*{20}c} {{\tau _{{\rm{GW}}}} = {5 \over {64}}{{{a^4}} \over {\mu {M^2}}} = {5 \over {64}}{{{a^4}} \over {q(1 + q)M_1^3}}}\quad\quad\quad\quad\quad\quad\quad\quad\\ {= 2.2 \times {{10}^8}{q^{- 1}}{{(1 + q)}^{- 1}}{{\left({{a \over {{R_ \odot}}}} \right)}^4}{{\left({{{{M_1}} \over {1.4{M_ \odot}}}} \right)}^{- 3}}{\rm{yr}},}\\ \end{array}$$

(1)

where M₁, M₂, and M ≡ M₁ + M₂ are the individual NS masses and the total mass of the binary, μ = M₁M₂/M is the reduced mass, q = M₂/M₁ is the binary mass ratio, and we assume geometrized units where G = c =1 (as we do throughout this paper, unless otherwise noted). The timescale for an elliptical orbit is shorter, and it can be shown that eccentricity is reduced over time by GW emission, leading to a circularization of orbits as they decay. A quick integration shows that the time until merger is given by τ_merge = τ_GW/4.

The luminosity of such systems in gravitational radiation is

$$\begin{array}{*{20}c} {{L_{{\rm{GW}}}} = - {{d{E_{{\rm{GW}}}}} \over {dt}} = {{32} \over 5}{{{\mu ^2}{M^3}} \over {{a^5}}} = {{32} \over 5}{{m_1^2m_2^2({m_1} + {m_2})} \over {{a^5}}}}\quad\quad\quad\quad\\ {= 5.34 \times {{10}^{32}}{q^2}(1 + q){{\left({{{{M_1}} \over {1.4{M_ \odot}}}} \right)}^5}{{\left({{a \over {{R_ \odot}}}} \right)}^{- 5}}{\rm{erg}}/{\rm{s}}}\quad\;\;\\ {= 8.73 \times {{10}^{51}}{q^2}(1 + q){{\left({{{{M_1}} \over {1.4{M_ \odot}}}} \right)}^5}{{\left({{a \over {100\;{\rm{km}}}}} \right)}^{- 5}}{\rm{erg}}/{\rm{s,}}}\\ \end{array}$$

(2)

which, at the end of a binary’s lifetime, when the components have approached to within a few NS radii of each other, is comparable to the luminosity of all the visible matter in the universe (∼ 10⁵³ erg/s). The resulting strain amplitude observed at a distance D from the source (assumed to be oriented face-on) is given approximately by

$$h = {{4{M_1}{M_2}} \over {aD}} = 5.53 \times {10^{- 23}}q{\left({{{{M_1}} \over {1.4{M_ \odot}}}} \right)^2}{\left({{a \over {100\;{\rm{km}}}}} \right)^{- 1}}{\left({{D \over {100\;{\rm{Mpc}}}}} \right)^{- 1}},$$

(3)

at a characteristic frequency

$${f_{GW}} = 2{f_{{\rm{orb}}}} = {1 \over \pi}\sqrt {{M \over {{a^3}}}} = 194{\left({{M \over {2.8{M_ \odot}}}} \right)^{1/2}}{\left({{a \over {100\;{\rm{km}}}}} \right)^{- 3/2}}{\rm{Hz}}.$$

(4)

The first measurement that will likely be made with direct GW observations is the orbital decay rate, with the period evolving (for the circular case) according to the relation

$${{dT} \over {dt}} = - {{192\pi} \over 5}{({{\mathcal M}_c}\omega)^{5/3}},$$

(5)

where T is the orbital period and ω the angular frequency, and thus the “chirp mass,”

$${{\mathcal M}_c} \equiv {\mu ^{3/5}}{M^{2/5}} = M_1^{3/5}M_2^{3/5}{({M_1} + {M_2})^{- 1/5}},$$

(6)

is likely to be the easiest parameter to determine from GW observations.

Several NS-NS systems are now known, including PSR J0737-3039 [59], a binary consisting of two observed pulsars, which allows for the prospect of even more stringent tests of GR [149]. Even with the handful of observed sources to date, one may use this sample to place empirical limits on the expected rate of NS-NS mergers [143] and to constrain the many parameters that enter into population synthesis calculations [217]. With regard to the former, the very short merger timescale for J0737, τ_merge = 85 Myr, makes it especially important for estimating the overall rate of NS-NS mergers since it is a priori very unlikely to detect a system with such a short lifetime.

Although black hole-neutron star (BH-NS) binaries are expected to form through the same processes as NS-NS binaries, none has been detected to date. This is generally thought to reflect their lower probability of detection in current surveys, in addition to intrinsically smaller numbers compared to NS-NS systems [39]. BH-NS systems are an expected byproduct of binary stellar evolution, and properties of the population may be inferred from population synthesis studies calibrated to the observed NS-NS sample (see, e.g., [37]).

In this review, we will summarize the current state of research on relativistic mergers, beginning in Section 2 with a description of the astrophysical processes that produce merging binaries and the expected parameters of these systems. The phases of the merger are briefly described in Section 3. In Section 4, we discuss the numerical techniques used to generate quasi-equilibrium (QE) sequences of NS-NS configurations, and we summarize the QE calculations that have been performed. These sequences yield a lot of information about NS physics, particularly with regard to the nuclear matter equation of state (EOS). They also serve as initial data for dynamical merger calculations, which we discuss next, focusing in turn on the numerical hydrodynamics techniques used to compute mergers and the large body of results that has been generated, in Sections 5 and 6, respectively. We pay particular attention to how numerical studies have taken steps toward answering a number of questions about the expected GW and electromagnetic (EM) emission from merging binaries, and we discuss briefly the possibility that they may be the progenitors of SGRBs and a source of r-process elements. We close with conclusions and a look to the future in Section 7.

While most of this review focuses on NS-NS mergers, many of the methods used to study NS-NS binaries are also used to evolve BH-NS binaries, and it has become clear that both merger types may produce similar observational signatures as well. For a review focusing on BH-NS merger calculations, we encourage the reader to consult the recent work by Shibata and Taniguchi [284].

2 Evolutionary Channels and Population Estimates

Merging NS-NS and BH-NS binaries, i.e., those for which the merger timescale is smaller than the Hubble time, are typically formed through similar evolutionary channels in stellar field populations of galaxies [37] (both may also be formed through dynamical processes in the high-density cores of some star clusters, but the overall populations are smaller and more poorly constrained; see [258] for a review). It is difficult to describe the evolutionary pathways that form NS-NS binaries without discussing BH-NS binaries as well, and it is important to note that the joint distribution of parameters such as merger rates and component masses that we could derive from simultaneous GW and EM observations will constrain the underlying physics of binary stellar evolution much more tightly than observing either source alone.

Population synthesis calculations for both merging NS-NS and BH-NS binaries typically favor the standard channel in which the first-born compact object goes through a common-envelope (CE) phase, although other models have been proposed, including recent ones where the progenitor binary is assumed to have very nearly-equal mass components that leave the main sequence and enter a CE phase prior to either undergoing a supernova [42, 54]. Simulations of this latter process have shown that close NS-NS systems could indeed be produced by twin giant stars with core masses ≳ 0.15 M_⊙, though twin main sequence stars typically merge during the contact phase [176].

In the standard channel (see, e.g., [44, 178], and Figure 1 for an illustration of the process), the progenitor system is a high-mass binary (with both stars of mass M ≳ 8–10 M_⊙ to ensure a pair of supernovae). The more massive primary evolves over just a few million years before it leaves the main sequence, passes through its giant phase, and undergoes a Type Ib, Ic, or II supernova, leaving behind what will become the heavier compact object (CO): the BH in a BH-NS binary or the more massive NS in an NS-NS binary. The secondary then evolves off the main sequence in turn, triggering a CE phase when it reaches the giant phase and overflows its Roche lobe. Dynamical friction shrinks the binary separation dramatically, until sufficient energy is released to expel the envelope. Without this step, binaries would remain too wide to merge through the emission of GWs within a Hubble time. Eventually, the exposed, Helium-rich core of the secondary undergoes a supernova, either unbinding the system or leaving behind a tight binary, depending on the magnitude and orientation of the supernova kick.

Figure 1

Cartoon showing standard formation channels for close NS-NS binaries through binary stellar evolution. Image reproduced from [178].

This evolutionary pathway has important effects on the physical parameters of NS-NS and BH-NS binaries, leading to preferred regions in phase space. The primary, which can accrete some matter during the CE phase, or during an episode of stable mass transfer from the companion Helium star, should be spun up to rapid rotation (see [178] for a review). In NS-NS binaries, we expect that this process will also reduce the magnetic field of the primary down to levels seen in “recycled” pulsars, typically up to four orders of magnitude lower than for young pulsars [180, 73]. The secondary NS, which never undergoes accretion, is likely to spin down rapidly from its nascent value, but is likelier to maintain a stronger magnetic field.

While this evolutionary scenario has been well studied for several decades, many aspects remain highly uncertain. In particular:

• The CE efficiency, which helps to determine the expected range of binary separations and the mass of the primary compact object after its accretion phase, remains very poorly constrained [229, 37, 78, 338]. If too much matter is accreted by the NS, it may undergo accretion-induced collapse to a BH [42], though the growing body of observed NS-NS systems does help place constraints on the allowed range of accretion-related parameters.

• The exact relation between a star’s initial mass and the eventual compact object mass is better understood, but significant theoretical uncertainties remain, and the relation is sensitive to the metallicity (often unknown), mainly through the effects of mass loss in stellar winds, and to the details of the explosion itself [336, 205].

• The maximum allowed NS mass will affect whether the primary remains a NS or undergoes accretion-induced collapse to a BH; its value is dependent upon the as-yet undetermined nuclear matter EOS. At present, the strongest limit is set by the binary millisecond pulsar PSR J1614-2230, for which a mass of M_NS = 1.97 ± 0.04 M_⊙ was determined by Shapiro time delay measurements [81]. Determining the NS EOS from GW observations may eventually provide stronger constraints [237, 221, 184], including a determination of whether supernova remnants are indeed classical hadronic NS or instead have cores consisting of some form of strange quark matter or other elementary particle condensate [223, 119, 230, 120, 23, 5].

• The supernova kick velocity distribution is only partially understood, leading to uncertainties as to which systems will become unbound after the second explosion [133, 324, 219, 152].

Given all these uncertainties, it is reassuring that most estimates of the NS-NS and BH-NS merger rate, expressed either as a rate of mergers per Myr per “Milky Way equivalent galaxy” or as a predicted detection rate for LIGO (the Laser Interferometer Gravitational-Wave Observatory) and Virgo (see Section 5.5 below), agree to within 1–2 orders of magnitude, which is comparable to the typical uncertainties that remain once all possible sources of error are folded into a population synthesis model. In Table 1, we show the predicted detection rates of NS-NS and BH-NS mergers for both the first generation LIGO detectors (“LIGO”), which ran at essentially their design specifications [2], and the Advanced LIGO (“AdLIGO”) configuration due to go online in 2015 [292]. We note that the methods used to generate these results varied widely. In [143], the authors used the observed parameters of close binary pulsar systems to estimate the Galactic NS-NS merger rate empirically (such results do not constrain the BH-NS merger rate). In [198, 128], the two groups independently estimated the binary merger rate from the observed statistics of SGRBs. In these cases, one does not get an independent prediction for the NS-NS and BH-NS merger rate, but rather some linear combination of the two. In both cases, the authors estimated that, if NS-NS and BH-NS mergers are roughly equal contributors to the observed SGRB sample, LIGO will detect about an order of magnitude more BH-NS mergers since their higher mass allows them to be seen over a much larger volume of the Universe. As they both noted, should either type of system dominate the SGRB sample, we would expect a doubling of LIGO detections for that class, and lose our ability to constrain the rate of the other using this method. Many population synthesis models have attempted to understand binary evolution within our galaxy by starting from a basic parameter survey of the various assumptions made about CE evolution, supernova kick distributions, and other free parameters. In [323, 79], population synthesis models are normalized by estimates of the star formation history of the Milky Way. In [140, 218], parameter choices are judged based on their ability to reproduce the observed Galactic binary pulsar sample, which allows posterior probabilities to be applied to each model in a Bayesian framework. A review by the LIGO collaboration of this issue may be found in [1].

Table 1 Estimated initial and advanced LIGO rates for BH-NS and NS-NS mergers from population synthesis calculations and other methods. The methods used are, in order, empirical constraints from the observed sample of binary pulsars (‘Empirical’), constraints on the combined NS-NS/BH-NS merger rate assuming that they are the progenitors of short-hard gamma-ray bursts (‘SGRBs’), population synthesis models calibrated to the star formation rate in the Milky Way (‘Pop. Synth. — SFR’), and population synthesis calibrated against the observed Galactic binary pulsar sample (‘Pop. Synth. — NS-NS’). We note that observations of binary pulsars do not yield constraints for BH-NS binaries. SGRB observations may produce constraints on NS-NS merger rates, BH-NS, merger rates, or both, depending on which sources are the true progenitors, but this remains unclear. Therefore, the table quotes results assuming a roughly equal split between the two. The official review of these results and their implications by the LIGO/Virgo Scientific Collaborations may be found in [1].

Should the next generation of GW interferometers begin to detect a statistically significant number of merger events including NSs, it should be possible to constrain several astrophysical parameters describing binary evolution much more accurately. These include

• The relative numbers of BH-NS and NS-NS mergers: Interferometric detections are sensitive to a binary’s “chirp mass” (see Eq. 6), and to the binary mass ratio as well [57, 11, 321, 320] if the signal-to-noise ratio is sufficiently high. Even if the merger signal takes place at frequencies too high to fall within the LIGO band, it should still be possible in most cases to determine whether the primary’s mass exceeds the maximum mass of a NS.

• The mass ratio probability distribution for merging binary systems: If both binary components’ masses are determined, we will be able to constrain both the BH mass distribution in merging binaries and the NS mass ratio distribution. Knowledge of the former would determine, e.g., whether the current low estimates for the mass accreted onto the primary CO core during the CE phase are correct [38], as this model predicts that BH masses in close BH-NS binaries should cluster narrowly around M_BH = 10 M_⊙. Previous calculations assuming larger accreted masses typically favored mass ratios closer to unity, since the primaries often began as NSs and underwent “accretion-induced collapse” to a BH during the CE phase. The NS-NS mass distribution will allow for tests of whether “twins”, i.e., systems whose zero-age main sequence (or “ZAMS”) masses are so close that they both leave the main sequence before either undergoes a supernova, play a significant role in the merging NS-NS population [42].

While Advanced LIGO or another interferometer will likely be required to make the first direct observations of NS-NS mergers and their immediate aftermath, it is possible that more than just the high-energy prompt emission from mergers may be observable using EM telescopes. Although the particular candidate source they identified resulted from a pointing error [110], Nakar and Piran suggest that the mass ejection from mergers should yield an observable radio afterglow [200], although the afterglows may be too faint to be seen by current telescopes at the observed distances of existing localized SGRBs [190]. While such outbursts could also result from a supernova, the luminosity required would be an order of magnitude larger than those previously observed. Given the length and timescales characterizing radio bursts, no NS-NS simulation has been able to address the model directly, but it certainly seems plausible that the time-variable magnetic fields within a stable hypermassive remnant could generate enough EM energy to power the resulting radio burst [280]. If mergers produce sufficiently large ejecta masses, M_ej ≳ 10⁻³ M_⊙, r-process nuclear reactions may produce a “kilonova” afterglow one day after a merger with a V-band optical luminosity νŁ_ν ≈ 10⁴¹ erg/s (roughly 1000 times brighter than a classical nova) [191]. These potential EM observations of mergers are likely to spur further research into the amount and velocity of merger ejecta, which could then be coupled to a larger-scale astrophysical simulation of the potential optical and radio afterglows.

3 Stages of a Binary Merger

The qualitative evolution of NS-NS mergers, or indeed any compact binary merger, has long been understood, and may be divided roughly into inspiral, merger, and ringdown phases, each of which presents a distinct set of challenges for numerical modeling and detection. As a visual aid, we include a cartoon summary in Figure 2, originally intended to describe black hole-black hole (BH-BH) mergers, and attributed to Kip Thorne. Drawn before the advent of the supercomputer simulations it envisions, we note that merger waveforms for all compact binary mergers have proven to be much smoother and simpler than shown here. To adapt it to NS-NS mergers, we note that NSs are generally assumed to be essentially non-spinning, and that the “ringdown” phase may describe either a newly formed BH or a NS that survives against gravitational collapse.

Figure 2

Cartoon picture of a compact binary coalescence, drawn for a BH-BH merger but applicable to NS-NS mergers as well (although NSs are generally assumed to be non-spinning). Image adapted from Kip Thorne.

Summarizing the evolution of the system through the three phases:

• After a pair of supernovae yields a relatively tight NS-NS binary, the orbital separation decays over long timescales through GW emission, a phase that takes up virtually all of the lifetime of the binary except the last few milliseconds. During the inspiral phase, binary systems may be accurately described by QE formalisms, up until the point where the gravitational radiation timescale becomes comparable to the dynamical timescale. The evolution in time is well-described by PN expansions, currently including all terms to 3PN [47], though small deviations can arise because of finite-size effects, especially at smaller separations (see Eqs. 1–6 for the lowest-order 2.5PN expressions for circular inspirals).

• Once the binary separation becomes no more than a few times the radii of the two NSs, binaries rapidly become unstable. The stars plunge together, following the onset of dynamical instability, and enter the merger phase, requiring full GR simulations to understand the complicated hydrodynamics that ensues. According to all simulations to date, if the NS masses are nearly equal, the merger resembles a slow collision, while if the primary is substantially more massive than the secondary the latter will be tidally disrupted during the plunge and will essentially accrete onto the primary. Since the NSs are most likely irrotational just prior to merging, there is a substantial velocity discontinuity at the surface of contact, leading to rapid production of vortices. Meanwhile, some fraction of the mass may be lost through the outer Lagrange points of the system to form a disk around the central remnant. This phase yields the maximum GW amplitude predicted by numerical simulations, but with a signal much simpler and more quasi-periodic than in the original cartoon version. GW emission during the merger encodes important information about the NS EOS, particularly the GW frequency f_cut at which the binary orbit becomes unstable (see Eq. 4) resulting in a characteristic cutoff in GW emission at those frequencies. Meanwhile, the merger itself may generate the thermal energy that eventually powers a SGRB, which occurs when the neutrinos and anti-neutrinos produced by shock-heated material annihilate around the remnant to produce high-energy photons.

• Finally, the system will eventually settle into a new, dynamically stable configuration through a phase of ringdown, with a particular form for the GW signal that depends on the remnant’s mass and rotational profile. If the remnant is massive enough, it will be gravitationally unstable and collapse promptly to form a spinning BH. Otherwise it must fall into one of three classes depending on its total mass. Should the remnant mass be less than the maximum mass M_iso supported by the nuclear matter EOS for an isolated, non-rotating configuration, it will clearly remain stable forever. Instead a remnant that is “supramassive”, i.e., with a mass above the isolated stationary mass limit but below that allowed for a uniformly rotating NS (typically ≲ 1.2 M_iso, with very weak dependence on the EOS; see, e.g., [147, 70, 71], and references therein) may become unstable. Supramassive remnants are stable against gravitational collapse unless angular momentum losses, either via pulsar-type emission or magnetic coupling to the outer disk, can drive the angular velocity below the critical value for stability. If the remnant has a mass above the supramassive limit, it may fall into the hypermassive regime, where it is supported against gravitational collapse by rapid differential rotation. Hypermassive NS (HMNS) remnants can have significantly larger masses, depending on the EOS (see, e.g., [31, 275, 86, 293, 114]), and will survive for timescales much longer than the dynamical time, undergoing a wide variety of oscillation modes. They can emit GWs should a triaxial configuration yield a significant quadrupole moment, and potentially eject mass into a disk around the remnant. Eventually, some combination of radiation reaction and magnetic and viscous dissipation will dampen the differential rotation and lead the HMNS to collapse to a spinning BH, again with the possibility that it may be surrounded by a massive disk that could eventually accrete. The energy release during HMNS collapse provides the possibility for a “delayed” SGRB, in which the peak GW emission occurs during the initial merger event, but the gamma-ray emission, powered by the collapse of the HMNS to a BH, occurs significantly later.

  Most calculations indicate that a geometrically thick, lower-density, gravitationally bound disk of material will surround whatever remnant is formed. Such disks, which are geometrically thick, are widely referred to as “tori” throughout the literature, though there is no clear distinction between the two terms, and we will use “disk” throughout this paper to describe generically the bound material outside a central merger remnant. Such disks are expected to heat up significantly, and much of the material will eventually accrete onto the central remnant, possibly yielding observable EM emission. Given the low densities and relatively axisymmetric configuration expected, disks are not significant GW emitters. There may be gravitationally unbound outflow from mergers as well, though dynamical simulations neither confirm nor deny this possibility yet. Such outflows, which can be the sites of exotic nuclear reactions, are frequently discussed in the context of r-process element formation, but their inherently low densities make them difficult phenomena to model numerically with high accuracy.

3.1 Comparison to BH-NS mergers

Since this three-phase picture is applicable to BH-NS mergers as well, it is worthwhile to compare the two merger processes at a qualitative level to understand the key similarities and differences. Inspiral for BH-NS mergers is also well-described by PN expansions up until shortly before the merger, but the parameter space is fundamentally different. First, since BHs are heavier than NSs, the dynamics can be quite different. Also, since BHs may be rapidly-spinning (i.e., have dimensionless spin angular momenta as large as J/M² ∼ 1), spin-orbit couplings can play a very important role in the orbital dynamics of the binary, imprinting a large number of oscillation modes into the GW signal (see, e.g., [126, 57]). From a practical standpoint, the onset of instability in BH-NS mergers should be easier to detect for LIGO, Virgo, and other interferometers, since the larger mass of BH-NS binaries implies that instability occurs at lower GW frequencies (see Eq. 4, noting that the separation a at which mass-transfer begins scales roughly proportionally to the BH mass).

The onset of instability of a BH-NS binary is determined by the interplay of the binary mass-ratio, NS compactness, and BH spin, with the first of these playing the largest role (see Figures 13–15 of [302] and the summary in [284]). In general, systems with high BH masses and/or more compact NSs tend to reach a minimum in the binding energy as the radius increases, leading to a dynamical orbital instability that occurs near the classical innermost stable circular orbit (ISCO). In these cases, the NS plunges toward the BH before the onset of tidal disruption, and is typically “swallowed whole”. leaving behind almost no material to form a disk. The GW emission from such systems is sharply curtailed after the merger event, yielding only a low-amplitude, high-frequency, rapidly-decaying “ringdown” signal (see, e.g., [154]). Numerical calculations have shown that even in borderline cases between dynamical instability and mass-shedding the NS is essentially swallowed whole, especially in cases where the BH in either non-spinning or spinning in the retrograde direction, which pushes the ISCO out to larger radii (see, e.g., [290, 289, 283, 91, 94]).

A richer set of phenomena occurs when the BH-NS mass ratio is closer to unity, the NS is less compact, the BH has a prograde spin direction, or more generally, some combination of those factors. In such cases, the NS will reach the mass-shedding limit prior to dynamical instability, and be tidally disrupted. Unlike what was described in semi-analytic Newtonian models (see, e.g., [68, 228, 77]) and seen in some early Newtonian and quasi-relativistic simulations (see, e.g., [165, 166, 138], stable mass transfer, in which angular momentum transfer via mass-shedding halts the inspiral, has never been seen in full GR calculations, nor even in approximate GR models (see the discussion in [96]). Even so, unstable mass transfer can produce a substantial disk around the BH, though in every GR simulation to data the substantial majority of the NS matter is accreted promptly by the BH (see [284] for a detailed summary of all current results). The exact structure of the disk and its projected lifetime depend on the binary system parameters, with the binary mass ratio and spin both important in determining the disk mass and the BH spin orientation critical for determining both the disk’s vertical structure and its thermodynamic state given the shock heating that occurs during the NS disruption. In general, the mass and temperature of the post-merger disks are comparable to those seen in some NS-NS mergers, and inasmuch as either is a plausible SGRB progenitor candidate then both need to be viewed as such. To date, no calculation performed in full GR has found any bound and self-gravitating NS remnant left over after the merger, including both NS cores that survive the initial onset of mass transfer by recoiling outward (predicted for cases in which stable mass transfer was thought possible, as noted above) or those in which bound objects form through fragmentation of the circum-BH disk. Motivated by observations of wide-ranging timescales for X-ray flares in both long and short GRBs [225], the latter channel has been suggested to occur in collapsars [227] and mentioned in the context of mergers [243], possibly on longer timescales than current simulations permit. Even so, there is no analogue to the HMNS state that may result from NS-NS mergers, nor any mechanism for a delayed SGRB as provided by HMNS collapse.

3.2 Qualitative numerical results

Constructing QE sequences for a given set of NS parameters requires sophisticated numerical schemes, but not supercomputer-scale resources, as we discuss in Section 4 below, focusing first on the numerical techniques used to construct QE binary data in GR, and the astrophysical information contained in the GW emission during the inspiral phase. Merger and ringdown, on the other hand, typically require large-scale numerical simulations, including some of the largest calculations performed at major supercomputer centers, as we discuss in detail Section 5 and 6 below.

To illustrate the various physical processes that occur during NS-NS mergers, we show the evolution of three different NS-NS merger simulations in Figures 3, 4, and 5, taken from Figures 4–6 of [144]. In Figure 3, we see the merger of two equal-mass NSs, each of mass M_NS = 1.4 M_⊙, described by the APR (Akmal-Pandharipande-Ravenhall) EOS [3]. In the second panel, clear evidence of “tidal lags” is visible shortly after first contact, leading to a slightly off-center collision pattern. By the third panel, an ellipsoidal HMNS has been formed, surrounded by a disk of material of lower density, which gradually relaxes to form a more equilibrated HMNS by the final panel. In Figure 4, we see a merger of two slightly heavier equal-mass NSs with M_NS = 1.5 M_⊙. In this case, the deeper gravitational potential limits the amount of mass that goes into the disk, and once a BH is formed (with a horizon indicated by the dashed blue circle in the final panel) it accretes virtually all of the rest mass initially present in the two NSs, with only 0.004% of the total remaining outside the horizon.

Figure 3

Isodensity contours and velocity profile in the equatorial plane for a merger of two equal-mass NSs with M_ns = 1.4 M_⊙ assumed to follow the APR model [3] for the NS EOS. The hypermassive merger remnant survives until the end of the numerical simulation. Image reproduced by permission from Figure 4 of [144], copyright by APS.

Figure 4

Isodensity contours and velocity profile in the equatorial plane for a merger of two equal-mass NSs with M_ns = 1.5 M_⊙ assumed to follow the APR model [3] for the NS EOS. With a higher mass than the remnant shown in Figure 3, the remnant depicted here collapses promptly to form a BH, its horizon shown by the dashed blue circle, absorbing all but 0.004% of the total rest mass from the original system. Image reproduced by permission from Figure 5 of [144], copyright by APS.

Figure 5

Isodensity contours and velocity profile in the equatorial plane for a merger of two unequal-mass NSs with M₁ = 1.3 M_⊙ and M₂ = 1.6 M_⊙, with both assumed to follow the APR model [3] for the NS EOS. In unequal-mass mergers, the lower mass NS is tidally disrupted during the merger, forming a disk-like structure around the heavier NS. In this case, the total mass of the remnant is sufficiently high for prompt collapse to a BH, but 0.85% of the total mass remains outside the BH horizon at the end of the simulation, which is substantially larger than for equal-mass mergers with prompt collapse (see Figure 4). Image reproduced by permission from Figure 6 of [144], copyright by APS.

Figure 6

Dimensionless binding energy E_b/M₀ vs. dimensionless orbital frequency M₀Ω, where M₀ is the total ADM (Arnowitt-Deser-Misner) mass of the two components at infinite separation, for two QE NS-NS sequences that assume a piecewise polytropic NS EOS. The equal-mass case assumes M_NS = 1.35 M_⊙ for both NSs, while the unequal-mass case assumes M₁ = 1.15 M_⊙ and M₂ = 1.55 M_⊙. The thick curves are the numerical results, while the thin curves show the results from the 3PN approximation. The lack of any minimum suggests that instability for these configurations occurs at the onset of mass shedding, and not through a secular orbital instability. Image reproduced by permission from Figure 16 of [305], copyright by AAS.

In Figure 5, we see the merger of an unequal-mass binary, with masses M₁ = 1.3 M_⊙ and M₂ = 1.6 M_⊙. In this case, the heavier NS partially disrupts the lighter NS prior to merger, leading to the secondary NS being accreted onto the primary. In this case, a much more massive disk is formed and, even after a BH forms in the center of the remnant, a substantial amount of matter, representing 0.85% of the total mass, remains outside the horizon. Later accretion of this material could potentially release the energy required to power a SGRB.

4 Initial Data and Quasi-Equilibrium Results

4.1 Overview

While dynamical calculations are required to understand the GW and EM emission from BH-NS and NS-NS mergers, some of the main qualitative features of the signals may be derived directly from QE sequences. From the variation of total system energy with binary angular velocity along a given sequence, it is possible to construct an approximate GW energy spectrum dE_GW/df immediately from QE results, essentially by performing a numerical derivative (see Figure 6). Doing so for a number of different sequences makes it possible to identify key frequencies where tidal effects may become measurable and to identify these with binary parameters such as the system mass ratio and NS radius. Similarly, since QE sequences should indicate whether a binary begins to shed mass prior to passage through the ISCO (see Figure 7), one may be able to classify observed signals into mass-shedding and non-mass-shedding events, and to use the critical point dividing those cases to help constrain the NS EOS. Single-parameter estimates have been derived for NS-NS binaries using QE sequences [98] (and for BH-NS binaries using QE [301] and dynamical calculations [283]). NS-NS binaries typically approach instability at frequencies f_GW ≳ 1 kHz, where laser shot noise is severely degrading the sensitivity of an interferometer detector. To observe ISCO-related effects with higher signal-to-noise, it may be necessary to operate GW observatories using narrow-band signal recycling modes, in which the sensitivity in a narrow range of frequencies is enhanced at the cost of lower sensitivity to broadband signals [56].

Figure 7

Mass-shedding indicator $\chi \equiv {\left({{{\partial (\ln \,h)} \over {\partial r}}} \right)_{{\rm{eq}}}}/{\left({{{\partial (\ln \,h)} \over {\partial r}}} \right)_{{\rm{pole}}}}$ vs. orbital frequency M₀Ω, where h is the fluid enthalpy and the derivative is measured at the NS surface in the equatorial plane toward the companion and toward the pole in the direction of the angular momentum vector, for a series of QE NS-NS sequences assuming equal-mass components. Here, χ = 1 corresponds to a spherical NS, while χ = 0 indicates the onset of mass shedding. More massive NSs are more compact, and thus able to reach smaller separations and higher angular frequencies before mass shedding gets underway. Image reproduced by permission from Figure 19 of [305], copyright by AAS.

It is important to note that, while the potential parameter space for NS EOS models is still very large, a much smaller set may serve to classify models for comparison with the first generation of GW detections. Indeed, by breaking up the EOS into piecewise polytropic segments, one may use as few as four parameters to roughly approximate all known EOS models, including standard nuclear models as well as models with kaon or other condensates [237]. To illustrate this, we show in Figure 8 four different QE models for NS-NS configurations with different EOS, taken from [305]; all have M₁ = 1.15 M_⊙ and M₂ = 1.55 M_⊙, and they correspond to the closest separation for which the QE code still finds a convergent result.

Figure 8

Isodensity contours for QE models of NS-NS binaries. In each case, the two NSs have masses M₁ = 1.15 M_⊙ (left) and M₂ = 1.55 M_⊙ (right), and the center-of-mass separation is as small as the QE numerical methods allow while able to find a convergent result. The models assume different EOS, resulting in different central concentrations and tidal deformations. Image reproduced by permission from Figures 9–12 of [305], copyright by AAS.

The inspiral of NS-NS binaries may yield complementary information about the NS structure beyond what can be gleaned from QE studies of tidal disruption. NSs have a wide variety of oscillation modes, including f-modes, g-modes, and r-modes, any of which may be excited by resonances with the orbital frequency as the latter sweeps upward. Should a particular oscillation mode be excited resonantly, it can then serve briefly as an energy sink for the system, potentially changing the phase evolution of the binary. For example, in a rapidly spinning NS, excitation of the m =1 r-mode can be significant, yielding a change of over 100 radians for the pre-merger GW signal phase in the case of a millisecond spin period [161]. For NS-NS mergers in the field, this would require one of the NSs to be a young pulsar that has not yet spun down significantly, which is unlikely because of the difficulty in obtaining such an extremely small binary separation after the second supernova. Other modes, such as the l = 2 f-mode, may be excited in less extreme circumstances, also yielding information about NS structure parameters [105].

4.2 Quasi-equilibrium formalisms

It has long been known that the GW emission from eccentric binaries is very efficient at radiating away angular momentum relative to the radiated energy [226]; as a result, the orbital eccentricity decreases as a binary inspirals, so that orbits should be very nearly circular long before they enter the detection range of ground-based interferometers. The only exception could be from a dynamical capture process that would create a binary with a significant eccentricity and very small orbital separation. Such eccentric binaries have been predicted to form in the nuclear cluster of our galaxy (see, e.g., [213]) and in core-collapsed globular clusters [127, 167]. However, at present, no formalism exists to construct initial data for such systems, besides superposing the individual components with sufficiently large initial separations to minimize constraint violations.

Using this circularity of primordial binaries as a starting point, one may use the constraint equations of GR, along with an assumption of quasi-circularity, to derive sets of elliptic equations describing compact binary configurations. For both QE and dynamical calculations, most groups typically make use of the Arnowitt-Deser-Misner (ADM) 3+1 splitting of the metric [9], which foliates the metric into a set of three-dimensional hypersurfaces by introducing a time coordinate. The resulting form of the metric, which is completely general, is written

$${g_{\mu \nu}} \equiv (- {\alpha ^2} + {\beta _i}{\beta ^i})d{t^2} + 2{\beta _i}dt\;d{x^i} + {\gamma _{ij}}d{x^i}d{x^j},$$

(7)

where α is known as the lapse function, β_i the shift vector, and γ_ij the spatial three-metric intrinsic to the hypersurface. We are following the standard relativistic notation here where Greek indices correspond to four-dimensional quantities and Latin indices to spatial three-dimensional quantities. Thus, the shift vector is a 3-vector, raised and lowered with the spatial 3-metric γ_ij rather than the spacetime 4-metric g_μν. To simplify matters, one typically defines a conformal factor ψ that factors out the determinant of the 3-metric, such that

$$\psi \equiv {[\det ({\gamma _{ij}})]^{1/12}},$$

(8)

introducing the conformal 3-metric ${\tilde \gamma _{ij}} \equiv {\psi ^{- 4}}{\gamma _{ij}}$ with unit determinant. While the 3-metric is a fundamental component of the geometric structure of the spacetime, the lapse function and shift vector are gauge quantities that simply reflect our choice of coordinates. Thus, while one often determines the lapse and shift in order to construct a appropriately “stationary” solution in the relevant coordinates between neighboring time slices, their values are often replaced to initialize dynamical runs with more convenient choices and thus different assumptions about coordinate evolution in time.

The field equations of general relativity take the deceptively simple form

$${G_{\mu \nu}} \equiv {R_{\mu \nu}} - {1 \over 2}{g_{\mu \nu}}R = 8\pi {T_{\mu \nu}},$$

(9)

where G_μν is the Einstein tensor, R_μν and R the Ricci curvature tensor and the curvature scalar, and T_μν the stress-energy tensor that accounts for the presence of matter, electromagnetic fields, and other physical effects that contribute to the mass-energy of the spacetime. Since GR is a second-order formulation, valid initial data must include not only the metric but also its first time derivative. It generally proves most convenient to introduce the time derivative of the metric after subtracting away the Lie derivative with respect to the shift, yielding a quantity known as the extrinsic curvature, K_ij:

$$({\partial _t} - {{\mathcal L}_\beta}){\gamma _{ij}} \equiv - 2\alpha {K_{ij}}.$$

(10)

Both the 3-metric and extrinsic curvature are symmetric tensors with six free parameters.

For systems containing NSs, one must consider the effects of nuclear matter through its presence in the stress-energy tensor T^μν. It is common to assume that the matter has the EOS describing a perfect, isotropic fluid, for which the stress energy tensor is given by

$${T^{\mu \nu}} \equiv ({\rho _0} + {\rho _0}\varepsilon + P){u^\mu}{u^\nu} + P{g^{\mu \nu}},$$

(11)

where ρ₀, ε, P and u^μ are the fluid’s rest-mass density, specific internal energy, pressure, and 4-velocity, respectively. Many calculations further assume that the NS EOS is described by an adiabatic polytrope, for which

$$P = (\Gamma - 1){\rho _0}\varepsilon = k\rho _0^\Gamma,$$

(12)

where Γ is the adiabatic index of the gas and k a constant, though a number of models designed to incorporate nuclear physics and/or strange matter condensates have also been constructed and studied (see Sections 4.4 and 6 below).

The problem in constructing initial data is not so much producing solutions that are self-consistent within GR, but rather to specify a sufficient number of assumptions to fully constrain a solution. Indeed, there are only four constraints imposed by the equations of GR, known as the Hamiltonian and momentum constraints. The Hamiltonian constraint is found by projecting Einstein’s equations twice along the direction defined by a normal observer, and describes the way stress-energy leads to curvature in the metric (see, e.g., [28] for a thorough review):

$$R + {K^2} - {K_{ij}}{K^{ij}} = 16\pi \rho,$$

(13)

where R is the scalar curvature of the 4-metric, $K = K_i^i$ is the trace of the extrinsic curvature, and

$$\rho \equiv {\bf{n}}\cdot{\bf{T}}\cdot{\bf{n}} = {\alpha ^2}{T^{00}} = \rho h{(\alpha {u^0})^2} - P$$

(14)

is the total energy density seen by a normal observer. The third term indicates that the total energy density is found by projecting the stress-energy tensor in the direction of the unit-length timelike normal vector n, whose components are given by

$${n_\mu} = (- \alpha,0,0,0){.}$$

(15)

In the final expression h ≡ 1 + ε + P/ρ₀ is the specific enthalpy of the fluid, and the combination Γ_n ≡ αu⁰ represents the Lorentz factor of the matter seen by an inertial observer. The notation here makes use of the standard summation convention, in which repeated indices are summed over.

Projecting Einstein’s equations in the space and time directions leads to the vectorial momentum constraint

$${D_i}K_j^i - {D_i}K = 8\pi {j_i},$$

(16)

where D_i represents a three-dimensional covariant derivative and j_i ≡ ρ₀hΓ_nu_i is the total momentum seen by a normal observer.

4.2.1 The Conformal Thin Sandwich formalism

In order to specify all the free variables that remain once the Hamiltonian and momentum constraints are satisfied, two different techniques have been widely employed throughout the numerical relativity community. One, known as the Conformal Transverse-Traceless (CTT) decomposition, underlies the Bowen-York [52] solution for black holes with known spin and/or linear momentum that is widely used in the “moving puncture” approach. To date, however, the CTT formalism has not been used to generate NS-NS initial data, and we refer readers to [284, 63] for descriptions of the CTT formalism applied to BH-NS and BH-BH initial data, respectively.

To date, most groups have used the Conformal Thin Sandwich (CTS) formalism to generate QE NS-NS data (see [28] for a review, [13, 137] for the initial steps in the formulation, and [326, 327, 333, 69] for derivations of the form in which it is typically used today). One first specifies that the conformal 3-metric is spatially flat, i.e., ${{\tilde \gamma}_{ij}} = {\delta _{ij}}$, where δ_ij is the Kronecker delta function. Under this assumption, the only remaining parameter defining the spatial metric is the conformal factor ψ, which serves the role of a gravitational potential. Indeed, in the limit of weak sources, it is linearly related to the standard Newtonian potential. Next, one specifies that there exists a helical Killing vector, so that, as the configuration advances forward in time, all quantities remain unchanged when properly rotated at constant angular velocity in the azimuthal direction. This is sufficient to fix all but the trace of the extrinsic curvature, with the other components forced to satisfy the relation

$${K^{ij}} = - {1 \over {2\alpha {\psi ^4}}}\left[ {{\nabla ^i}{\beta ^j} + {\nabla ^j}{\beta ^i} - {2 \over 3}{\gamma ^{ij}}{\nabla _k}{\beta ^k}} \right].$$

(17)

The trace of the extrinsic curvature K remains a free parameter in this approach. While one may choose arbitrary prescriptions to fix it, most implementations choose a maximal slicing of the spatial hypersurfaces by setting K = ∂_tK = 0. Under these assumptions the Hamiltonian and momentum constraints, along with the trace of the Einstein equations, yield five elliptic equations for the lapse, shift vector, and conformal factor:

$${\nabla ^2}\psi = - {\psi ^5}\left({{1 \over 8}{K^{ij}}{K_{ij}} - 2\pi \rho} \right),$$

(18)

$${\nabla ^2}(\alpha \psi) = \alpha {\psi ^5}({7 \over 8}{\psi ^4}{K_{ij}}{K^{ij}} + 2\pi {\psi ^4}(\rho + 2S),$$

(19)

$${\nabla ^2}{\beta ^i} + {1 \over 3}{\nabla ^j}{\nabla _j}{\beta ^i} = 2\alpha {\psi ^4}{K^{ij}}{\nabla _j}(\alpha {\psi ^{- 6}}) + 16\alpha {\psi ^4}{j^i},$$

(20)

where

$$S = ({g_{\mu \nu}} + {n_\mu}{n_\nu}){T^{\mu \nu}} = S_j^j = 3P + (\rho + P)[1 - \Gamma _n^{- 2}]$$

(21)

is the trace of the stress-energy tensor projected twice in the spatial direction. While these five equations are linked and the right-hand sides are nonlinear, they are amenable to solution using iterative methods. Boundary conditions are set by assuming asymptotic flatness: at large radii, the metric takes on the Minkowski form so α → 1, ψ → 1, and $\beta _{{\rm{rot}}}^i \rightarrow \Omega \times \vec r$. We note that a purely corotating shift term yields zero when we apply the left-hand side of Eq. 20, so we may subtract it away and solve the equation with a boundary condition of zero instead.

The breakdown in Eqs. 18, 19, and 20 is not unique. The Meudon group [125, 124], to pick one example, has often chosen to define ν ≡ ln α and β ≡ ln(αψ²), and replace Eqs. 18 and 19 with the equivalent pair

$$\begin{array}{*{20}c} {{\nabla ^2}\nu = {\psi ^4}{K^{ij}}{K_{ij}} - {\nabla _i}\nu {\nabla ^i}\beta + 4\pi {\psi ^4}(\rho + S),}\quad\quad\quad\quad\\ {{\nabla ^2}\beta = {3 \over 4}{\psi ^4}{K^{ij}}{K_{ij}} - {1 \over 2}({\nabla _i}\nu {\nabla ^i}\nu + {\nabla _i}\beta {\nabla ^i}\beta) + 4\pi {\psi ^4}S.}\\ \end{array}$$

This approach is sufficient to define the field component of the configuration, but one still needs to solve for the matter quantities as well. One starts by assuming that there is a known prescription for reconstructing the density, internal energy, and pressure from the enthalpy h. Next, one has to assume some model for the spin of the NS. While corotation is often a simpler choice, since the velocity field of the matter is zero in the corotating frame, the more physically reasonable condition is irrotational flow. Indeed a realistic NS viscosity is never sufficiently large to tidally lock the NS to its companion during inspiral [45, 146]. If we define the co-momentum vector w_i = hu_i, irrotational flow implies the vanishing of its curl:

$$\nabla \times w = {\nabla _\mu}{w_\nu} - {\nabla _\nu}{w_\mu} = {\partial _\mu}{w_\nu} - {\partial _\nu}{w_\mu} = 0,$$

(22)

which allows us to define a velocity potential Ψ such that w ≡ ∇Ψ. Using these quantities, one may write down the integrated Euler equation

$$h\alpha {\Gamma _n}/(1 - {\gamma _{ij}}{U^i}U_0^j) = {\rm{const}}{.},$$

(23)

where the 3-velocity Uⁱ of the fluid with respect to an Eulerian observer is given by

$${U^i} \equiv {{{u^i} + {\beta ^i}{u^0}} \over {{\Gamma _n}}},$$

(24)

and the orbital 3-velocity $U_0^i$ with respect to the same observer by $U_0^i = {\beta ^i}/\alpha$. For details on the ways in which one may construct an elliptic equation for the velocity potential, we refer to the derivation in [125].

To date, all QE sequences and dynamical runs in the literature have assumed that NSs are either irrotational or synchronized, but it is possible to construct the equations for arbitrary NS spins so long as they are aligned [29, 309]. While suggestions are also given there on how to construct QE sequences with intermediate spins using the new formalism, none have yet appeared in the literature. Similarly, a formalism to add magnetic fields self-consistently to QE sequences has been constructed [314], as current dynamical simulations typically begin from data assuming either zero magnetic fields or those that only contribute via magnetic pressure.

4.2.2 Other formalisms

The primary drawback of the CTS system is the lack of generality in assuming the spatial metric to be conformally flat, which introduces several problems. The Kerr metric, for example, is known not to be conformally flat, and conformally flat attempts to model Kerr BHs inevitably include spurious GW content. The same problem affects binary initial data: in order to achieve a configuration that is instantaneously time-symmetric, one actually introduces spurious gravitational radiation into the system, which can affect both the measured parameters of the initial system as well as any resulting evolution.

Other numerical formalisms to specify initial data configurations in GR have been derived using different assumptions. Usui and collaborators derived an elliptic set of equations by allowing the azimuthal component of the 3-metric to independently vary from the radial and longitudinal components [319, 318], finding good agreement with the other methods discussed above. A number of techniques have been developed to construct helically symmetric spacetimes in which one actually solves Einstein’s equations to evaluate the non-conformally-flat component of the metric, which are typically referred to as “waveless” or “WL” formalisms [260, 50, 291]. In terms of the fundamental variables, rather than specifying the components of the conformal spatial metric by ansatz, one specifies instead the time derivative of the extrinsic curvature using a physically motivated prescription. These methods are designed to match the proper asymptotic behavior of the metric at large distances, and may be combined with techniques designed to enforce helical symmetry of the metric and gauge in the near zone the near zone helical symmetry, or “NHS” formalism) to produce a global solution [315, 334, 316]. QE sequences generated using this formalism [315] have shown that the resulting conformal metric is indeed non-flat, with deviations of approximately 1% for the metric components, and similar differences in the system’s binding energy when compared to equivalent CTS results. They suggest [316] that underestimates in the quadrupole deformations of NS prior to merger may result in total phase accumulation errors of a full cycle, especially for more compact NS models.

QE formalisms reflect the assumption that binaries will be very nearly circular, since GW emission acting over very long timescales damps orbital eccentricity to negligible values for primordial NS-NS binaries between their formation and final merger. Binaries formed by tidal capture and other dynamical processes, which may be created with much smaller initial separations, are more likely to maintain significant eccentricities all the way to merger (see, e.g., [213] for a discussion of such processes for BH-BH binaries) and it has been suggested based on simple analytical models that such mergers, likely occurring in or near dense star clusters, may account for a significant fraction of the observed SGRB sample [167]. However, more detailed modeling is required to work out accurate estimates of merger rates given the complex interplay between dynamics and binary star evolution that determines the evolution of dense star clusters, and given the large uncertainties in the distributions of star cluster properties in galaxies throughout the universe. No initial data have ever been constructed in full GR for merging NS-NS binaries with eccentric orbits since the systems are then highly time-dependent, while the calculations performed to evolve them generally use a superposition of two stationary NS configurations [122].

4.3 Numerical implementations

There are a number of numerical techniques that have been used to solve these elliptic systems. The first calculations of NS-NS QE sequences, in both cases for synchronized binaries, were performed by Wilson, Mathews, and Marronetti [327, 328, 189] and Baumgarte et al. [25, 26]. The former used a finite differencing scheme, and centered different quantities at cells, vertices, and faces in order to construct a system of equations that was solved using fast matrix inversion techniques, while the latter used a Cartesian multigrid scheme, restricted to an octant to increase computational efficiency.

After a formalism for evaluating QE irrotational NS-NS sequences was developed [51, 313], some of the first results were obtained by Uryū and Eriguchi, who developed a finite-differencing code in spherical coordinates allowing for the solution of relativistic NS-NS binaries using Green’s functions[313, 317]. Their method extended the self-consistent field (SCF) work of [147], which had previously been applied to axisymmetric configurations. Irrotational configurations were also generated by Marronetti et al. [185], using the same finite difference scheme as found in the work on synchronized binaries.

The most widely used direct grid-based solver in numerical relativity is the Bam_Elliptic solver [55], which solves elliptic equations on single rectangular grids or multigrid configurations. It is included within the Cactus code, which is widely used in 3-D numerical relativity [307]. In particular it has been used to initiate a number of single and binary BH simulations, including one of the original breakthrough binary puncture works [61].

Lagrangian methods, typically based on smoothed particle hydrodynamics (SPH) [181, 118, 194]) have been used to generate both synchronized and irrotational configurations for PN [10, 99, 101, 100] and conformally flat (CF) [211, 210, 97, 212, 207, 209, 208, 34, 35, 33] calculations of NS-NS mergers, but they have not yet been extended to fully GR calculations, in part because of the difficulties in evolving the global spacetime metric.

The most widely used data for numerical calculations are those generated by the Meudon group (see Section 4.4 below for details on their calculations and [125] for a detailed description of their methods). The code they developed, Lorene [124], uses multidomain spectral methods to solve elliptic equations (while the code has been used primarily for relativistic stellar and binary configurations, it can be used as a more general solver). Around each star, one creates a set of nested, contiguous grids, with points arrayed in the radial and angular directions. The innermost grid has spheroidal geometry, and the surrounding grids are annular. The outermost grid may be allowed to extend to spatial infinity through a compactification transformation of the radial coordinate. To solve elliptic equations for various field quantities, one breaks each into a sum of two components, each of whose source terms are concentrated in one NS or the other. Similarly, the source terms themselves are split into two pieces, ideally, so each component is well-described by spheroidal spectral coefficients centered around each star. Using the spectral expansion, one may pass values from one star to the other and then recalculate spectral coefficients for the other grid configuration. This scheme has several efficiency advantages over direct grid-based methods, which helps to explain its popularity. First, the domain geometry may be chosen to fit to a NS surface, which eliminates Gibbs phenomenon-related errors and allows for exponential convergence with respect to the number of grid points, rather than the geometric convergence that characterizes finite difference-based grid codes. Second, the use of spectral methods requires much less computer memory than grid-based codes, and, as a result, Lorene is a serial code that can run easily on any off-the-shelf PC, rather than requiring a supercomputer platform.

4.4 Quasi-equilibrium and pre-merger simulations

NS-NS binaries may be well approximated by QE configurations up until they reach separations comparable to the sizes of the binary components themselves, that latter phase lasting a fraction of a second after an inspiral of millions of years or more. The eventual merger will occur after the binary undergoes one of two possible orbital instabilities. If the total binary energy and angular momenta reach a minimum at some separation, which defines the ISCO, the binary becomes dynamically unstable and plunges toward merger. Alternately, if the NS fills its Roche lobe (typically the lower density NS) mass will transfer onto the primary and the secondary will be tidally disrupted. The parameters of some NS-NS systems could technically allow for stable mass transfer, in which mass loss from a lighter object to a heavier one leads to a widening of the binary separation. This does occur for some binaries containing white dwarfs, but every dynamical calculation to date using full GR or even approximate GR has found that the rapid inspiral rate leads to inevitably unstable mass transfer and the prompt merger of a binary.

Many of the results later confirmed using relativistic QE sequences were originally derived in Newtonian and PN calculations, particularly as explicit extensions of Chandrasekhar’s body of work (see [65]). Chandrasekhar’s studies of incompressible fluids were first extended to compressible binaries by Lai, Rasio, and Shapiro [156, 155, 158, 157, 159], who used an energy variational method with an ellipsoidal treatment for polytropic NSs. They established, among other results, the magnitude of the rapid inspiral velocity near the dynamical stability limit [156], the existence of a critical polytropic index (n ≈ 2) separating binary sequences undergoing the two different terminal instabilities [155], the role played by the NS spin and viscosity and magnitude of finite-size effects in relation to 1PN terms [158, 157], and the development of tidal lag angles as the binary approaches merger [159]. They also determined that for most reasonable EOS models and nonextreme mass ratios, as would pertain to NS-NS mergers, an energy minimum is inevitably reached before the onset of mass transfer through Roche lobe overflow. The general results found in those works were later confirmed by [201], who used a SCF technique [131, 132], finding similar locations for instability points as a function of the adiabatic index of polytropes, but a small positive offset in the radius at which instability occurred. Similar results were also found by [311, 312], but with a slight modification in the total system energy and decrease in the orbital frequency at the onset of instability.

The first PN ellipsoidal treatments were developed by Shibata and collaborators using self-consistent fields [270, 269, 279, 281, 299] and by Lombardi, Rasio, and Shapiro [177]. Both groups found that the nonlinear gravitational effects imply a decrease in the orbital separation (increase in the orbital frequency) at the instability point for more compact NS. This result reflects a fairly universal principle in relativistic binary simulations: as gravitational formalisms incorporate more relativistic effects, moving from Newtonian gravity to 1PN and on to CF approximations and finally full GR, the strength of the gravitational interaction inevitably becomes stronger. The effects seen in fully dynamical calculations will be discussed in Section 6, below.

The first fully relativistic CTS QE data for synchronized NS-NS binaries were constructed by Baumgarte et al. [26, 25], using a grid-based elliptic solver. Their results demonstrated that the maximum allowed mass of NSs in close binaries was larger than that of isolated NSs with the same (polytropic) EOS, clearly disfavoring the “star-crushing” scenario that had been suggested by [327, 187] using a similar CTS formalism (but see also the error in these latter works addressed in [104], discussed in Section 6.3 below). Baumgarte et al. also identified how varying the NS radius affects the ISCO frequency, and thus might be constrained by GW observations. Using a multigrid method, Miller et al. [192] showed that while conformal flatness remained valid until relatively near the ISCO, the assumption of syncronized rotation broke down much earlier. Usui et al. [319] used the Green’s function approach with a slightly different formalism to compute relativistic sequences and determined that the CTS conditions were valid up until extremely relativistic binaries were considered.

The first relativistic models of physically realistic irrotational NS-NS binaries were constructed by the Meudon group [51] using the Lorene multi-domain pseudo-spectral method code. Since then, the Meudon group and collaborators have constructed a wide array of NS-NS initial data, including polytropic NS models [125, 303, 304], as well as physically motivated NS EOS models [36] or quark matter condensates [170]. Irrotational models have also been constructed by Uryū and collaborators [313, 317] for use in dynamical calculations, and nuclear/quark matter configurations have been generated by Oechslin and collaborators [212, 209]. A large compilation of QE CTS sequences constructed using physically motivated EOS models including FPS (Friedman-Pandharipande) [222], SLy (Skyrme Lyon) [83], and APR [3] models, along with piecewise polytropes designed to model more general potential cases (see [237]), was published in [305].

The most extensive set of results calculated using the waveless/near-zone helical symmetry condition appear in [316], with equal-mass NS-NS binary models constructed for the FPS, Sly, and APR EOS in addition to Γ = 3 polytropes. Results spanning all of these QE techniques are summarized in Table 2.

Table 2 A summary of various studies focusing on QE sequences of NS-NS binaries. Please refer to Section 6 for a discussion of papers that focus on dynamical simulations instead. Gravitational schemes include Newtonian gravity (‘Newt.’), lowest-order post-Newtonian theory (‘PN’), conformal thin sandwich (‘CTS’) including modified forms of the spatial metric (‘Mod. CTS’), and waveless/near-zone helical symmetry techniques. Numerical methods include ellipsoidal formalisms (‘Ellips.’), self-consistent fields (‘SCF’), numerical grids (‘Grid’), multigrids, and multipatch, Green’s function techniques (‘Green’s’), spectral methods (‘Spectral’), or SPH relaxation (‘SPH’). With regard to EOS models, ‘WD’ refers to the exact white dwarf EOS assuming a cold degenerate electron gas [64]. The ‘Physical’ EOS models include the FPS [222], SLy [83], and APR [3] nuclear EOS models, along with their parameterized approximations and other physically motivated models. The compactness ${\mathcal C} = M/R$ refers to the value for a NS in isolation before it is placed in a binary, and plays no role in Newtonian physics. The mass ratio q = M₂/M₁ is defined to be less than unity, and ‘spin’ refers to either synchronized or irrotational configurations.

5 Dynamical Calculations: Numerical Techniques

5.1 Overview: General relativistic (magneto-)hydrodynamics and microphysical treatments

NS-NS binaries are highly relativistic systems, numerous groups now run codes that evolve both GR metric fields and fluids self-consistently, with some groups also incorporating an ideal magnetohydrodynamic evolution scheme that assumes infinite conductivity. The codes that evolve the GR hydrodynamics or magnetohydrodynamics (GRHD and GRMHD, respectively) equations are many and varied, incorporating different spatial meshes, relativistic formalisms, and numerical techniques, and we will summarize the leading variants here. All full GR codes now make use of the significant insight gained from BH-BH merger calculations, but much work on these systems predates the three 2005 “breakthrough” papers by Pretorius [231], Goddard [22], and the group then at UT Brownsville (now at RIT) [61], with the first successful NS-NS merger calculations announced already in 1999 [287]. A list of the groups that have performed NS-NS merger calculations using full GR is presented below; note that many of these groups have also performed BH-NS simulations, as discussed in the review by Shibata and Taniguchi [284].

Of the full GR codes used to evolve NS-NS binaries, almost all are grid-based and make use of some form of adaptive mesh refinement. The one exception is the SpEC code developed by the SXS collaboration, formed originally by Caltech and Cornell, which has used a hybrid spectral-method field solver with grid-based hydrodynamics. Most make use of the BSSN formalism for evolving Einstein’s equations (see Section 5.2.1 below), while the HAD code uses the alternate Generalized Harmonic Gauge (GHG) approach. This technique is also used by the SXS collaboration and the Princeton group, who have both performed simulations of merging BH-NS binaries (see Section 6.6) but have yet to report any results on NS-NS mergers. Three groups have reported results for NS-NS mergers including MHD (HAD, Whisky, and UIUC), while the KT (Kyoto/Tokyo) group has reported magnetized evolutions of HMNS remnants (see [280] and references therein for a discussion of their work and that of other numerical relativity groups), but have yet to use that code for a NS-NS merger calculation.

While full GR codes were being developed to study NS-NS binaries, a parallel and rather independent track developed to study detailed microphysical effects in binary mergers using approximate relativistic schemes. This includes codes like that developed by the MPG group that accurately track the production of neutrinos and antineutrinos and their annihilation during a merger, as well as post-processing routines that use extensive nuclear chains to track the production of rare high-atomic number r-process elements in merger ejecta [123]. Meanwhile, the Bremen group’s SPH code includes variable-temperature physically motivated equations of state [247] and magnetohydrodynamics [233], and has been used with a multi-group flux-limited diffusion neutrino code to generate expected neutrino signatures from merger calculations [82]. A summary of groups performing NS-NS merger calculations is presented in Table 3.

Table 3 A summary of groups reporting NS-NS merger calculation results. The asterisk for the KT collaboration’s MHD column indicates that they have used an MHD-based code for other projects, but not yet for NS-NS merger simulations. Gravitational formalisms include full GR, assumed to be implemented using the BSSN decomposition except for the HAD collaborations’s GHG approach, the CF approximation, or Newtonian gravity. Microphysical treatments include physically motivated EOS models or quark-matter EOS and neutrino leakage schemes.

5.2 GR numerical techniques

5.2.1 GR formalisms and gauge choice

There are two distinct schemes used in all binary merger calculations performed to date, the BSSN (Baumgarte-Shapiro-Shibata-Nakamura) [277, 27] and Generalized Harmonic formalisms. For general reviews of these formalisms, as well as other developments in numerical relativity, we refer readers to two recent books on numerical relativity [4, 30]. Here we only briefly summarize these two schemes.

The BSSN formalism was adapted from the 3+1 ADM approach, with quantities defined as in Eqs. 7 and 8. While the original ADM scheme proved to be numerically unstable, defining the auxiliary quantities ${{\tilde \Gamma}^i} = - \tilde \gamma _{,j}^{ij}$ and treating these expressions as independent variables stabilized the system and allowed for long-term evolutions. While slight variants exist, the 19 evolved variables are typically either the conformal factor ψ or its logarithm ϕ, the conformal 3-metric, the trace K of the extrinsic curvature, the trace free extrinsic curvature A_ij and the conformal connection functions ${{\tilde \Gamma}^i}$ The evolution equations themselves are given in Appendix A.

The BSSN scheme was used in the binary merger calculations of the KT collaboration [287, 288, 285, 286], the first completely successful NS-NS calculations ever performed in full GR. Ever since the UTB/RIT [61] and Goddard [22] groups showed simultaneously that a careful choice of gauge allows BHs to be evolved in BSSN schemes without the need to excise the singularity, these “puncture gauges” have gained widespread hold, and have been used to evolve NS-NS binaries (and in some cases, BH-NS binaries) by the KT collaboration [332], UIUC [172], and Whisky [17].

The generalized harmonic formalism, developed over about two decades from initial theoretical suggestions up to its current numerical implementation [112, 115, 232, 130, 231] was used to perform the first calculations of merging BH-BH binaries by Pretorius [231], and has since been applied to NS-NS binaries by the HAD collaboration [6] and to BH-NS mergers by HAD [66], SXS [85, 84, 108], and the Princeton group [294, 88]. It involves introducing a set of auxiliary quantities denoted H^μ representing the action of the wave operator on the spacetime coordinates themselves

$${H^\mu} \equiv {\Gamma ^\mu} \equiv {g^{\alpha \beta}}\Gamma _{\alpha \beta}^\mu = - {1 \over {\sqrt {- g}}}{\partial _\nu}(\sqrt {- g} {g^{\mu \nu}}) = {g^{\alpha \beta}}{\nabla _\alpha}{\nabla _\beta}{x^\mu} = \square{x^\mu},$$

(25)

which are treated as independent gauge variables whose evolution equation must be specified. Current first-order formulations [171, 6] evolve equations for the spacetime metric g_μν along with its spatial derivative, Φ_iμν = ∂_ig_μν and projected time derivative Π_μν = −n^α∂_αg_μν, subject to consistency constraints on the spatial derivatives. The first BH-NS merger calculations in the GH formalism used a first-order reduction [171] of the Einstein equations and specified the source functions to damp to zero exponentially in time, while the first binary NS merger work [6] used a similar first-order reduction and chose harmonic coordinates with H^μ = 0.

In both formalisms, most groups employ grid-based finite differencing to evaluate spatial derivatives. While finite differencing operators may be easily written down to arbitrary orders of accuracy, there is a trade-off between the computational efficiency achievable by using smaller, second-order stencils and the higher accuracy that can be attained using larger, higher-order ones. For the moment, many groups are now moving to at least fourth-order accurate differencing techniques, with a high likelihood that at least the field sector of NS-NS merger calculations will soon be performed at comparable order to BH-BH calculations, at either sixth [136] or eighth-order [179] accuracy, if not higher. The main limitation to date involves the complexity of shock capturing using higher-order schemes, as we discuss in Section 5.2.3 below.

Imposing physically realistic and accurate boundary conditions remains a difficult task for numerical codes. Ideally, one wishes to impose boundary conditions at large distances that preserve the GR constraints and yield a well-posed initial-boundary value problem. On physical grounds, the boundary should only permit outgoing waves, preventing the reflection of spurious waves back into the numerical grid. Otherwise, reflections may be avoided by placing the outer boundaries so far away from the matter that they remain causally disconnected from the merging objects for the full duration of the simulation. Building upon previous work (see, e.g., [113, 151, 12, 150, 256, 242], Winicour [329] derived boundary conditions that satisfy all of the desired conditions for the generalized harmonic formalism. No such conditions have been derived for the full BSSN formalism, though progress has been made (see, e.g., [43, 129]) so that we may now construct well-posed boundary conditions in the weak-gravity linearized limit of BSSN [204] and for related first-order gravitational formalisms [49].

5.2.2 Grid decompositions

While unigrid schemes are extremely convenient, they tend to be extremely inefficient, since one must resolve small-scale features in the central regions of a merger but also extend grids out to the point where the GW signal has taken on its roughly asymptotic form. Thus, nearly every code incorporates some means of focusing resolution on the high-density regions via some form of mesh refinement. A simple approach, for instance, is to use “fisheye” coordinates that represent a continuous radial deformation of a grid [172], a technique that had previously been used successfully, e.g., for BH-BH mergers [21, 62].

While fixed mesh refinement offers the chance for greater computational efficiency and accuracy, much current work focuses on adaptive mesh refinement, in which the level of refinement of the grid is allowed to evolve dynamically to react to the evolving binary configuration. Several codes, both public and private, now implement this functionality. The publicly available Carpet computational toolkit [263, 262], which is distributed to the community as part of the open source Einstein [90, 175] uses a “moving boxes” approach, and has been designed to be compatible with the widely used and publicly available Cactus framework. It has been successfully implemented by the Whisky group [17] to perform NS-NS mergers, by UIUC for their BH-NS mergers [94], and a host of other groups for BH-BH mergers and additional problems. The KT code SACRA also implements an adaptive mesh refinement (AMR) scheme for NS-NS and BH-NS mergers [332], as does the most recent version of the HAD collaboration’s code [6], which is based on the publicly available infrastructure of the same name [169, 8], and the BAM code employed by the Jena group [308, 41, 122]. The Princeton group also has an AMR code, which has been used to perform BH-NS mergers to date [294, 88, 89]

One drawback of employing rectangular grids is that memory costs scale like N³, where N is the number of grid cells across a side, and total computational effort like N⁴ once one accounts for the reduction in the timestep due to the Courant-Friedrich-Levy criterion. Since one does not necessarily need high angular resolution at large radii, there is great current interest in developing schemes that use some form of spheroidal grid, for which the memory scaling is merely ∝ N. A group at LSU has implemented a multi-patch method [337], in which space is broken up into a number of non-overlapping domains in such a way that the six outermost regions (projections of the faces of a cube onto spheres of constant radius), maintain constant angular resolution and thus produce linear dependence of the total number of grid points on the number of radial points. To date, it has been used primarily for vacuum spacetimes [224] and hydrodynamics on a fixed background. The SXS collaboration, begun at Caltech and Cornell and now including members at CITA and Washington State, has used a spectral evolution code with multiple domains to evolve BH-NS binaries, which achieves the same scaling by expanding the metric fields in modes rather than in position space [85]. Their first published results on NS-NS binaries are currently in preparation (see [141] for a brief summary of work to date).

While all of these grid techniques produce tremendous advantages in computational efficiency, each required careful study since deformations of a grid or the introduction of multiple domains can introduce inaccuracies and non-conservative effects. As an example, in AMR schemes, one must deal with the same reflection issues at refinement boundaries that are present at the physical boundaries of the grid, as discussed above, though the interior nature of the boundaries allows for additional techniques, such as tapered grid boundaries [168], to be used to minimize reflections there. The study of how to minimize spurious effects in these schemes continues, and will represent an important topic for years to come, especially as evolution schemes become more complicated by including more physical effects.

5.2.3 Hydrodynamics, MHD, and high-resolution shock capturing

Fluids cannot be treated in the same way as the spacetime metric because finite differencing operators do not return meaningful results when placed across discontinuities induced by shocks. Instead, GR(M)HD calculations must include some form of shock-capturing that accounts for these jumps. These are typically implemented in conservative schemes, in which the fluid variables are transformed from the standard “primitive” set ${\vec P}$, which includes the fluid density, pressure, and velocity (and magnetic field in MHD evolutions), into a new set ${\vec U}$ so that the evolution equations may be written in the form

$${\partial _t}(\vec U) = \nabla \cdot \vec F + \vec S,$$

(26)

where the flux functions $\vec F(\vec P)$ and source terms $\vec S(\vec P)$ can be expressed in terms of the primitive variables but not their derivatives. These schemes allow one to evolve the resulting MHD set of equations so that numerical fluxes are conserved to numerical precision across cell walls as the fluid evolves in time. One widely used scheme, often referred to as the Valencia formulation [186], is described in Appendix B..

There are important mathematical reasons for casting the GRHD/GRMHD system in conservative form, primarily since the mathematical techniques describing Godunov methods may be called into play [121]. In such methods, we assume that the evolution of the quantities ${\vec U}$ may be expressed in the form

$$\vec U(x = {x_i},t = {t_{n + 1}}) = \vec U(x = {x_i},t = {t_n}) + {{\Delta t} \over {\Delta x}}\left({\vec F(x = {x_{i - 1/2}}) - \vec F(x = {x_{i + 1/2}})} \right),$$

(27)

where the points have spatial coordinates x_i ≡ x₀ + iΔx and the timesteps satisfy t_n = t₀ + nΔt. The fluxes must be determined by solving the Riemann problem at each cell face (thus the half-integer indices), either exactly or approximately. It can be shown that solutions constructed this way, when convergent, must converge to a solution of the original problem, even when shocks are present [163].

First one reconstructs the primitives from the conserved variables on both sides of an interface, using interpolation schemes designed to respect the presence of shocks. Simple schemes involving limiters yield second-order accuracy in general but first-order accuracy at shocks, while higher-order methods such as PPM (piecewise parabolic method) and essentially non-oscillatory (ENO) schemes such as CENO (central ENO) and WENO (weighted ENO) yield higher accuracy but at much higher computational cost. Once primitives are interpolated to the cell interfaces, flux terms are evaluated there and one solves the Riemann problem describing the evolution of two distinct fluid configurations placed on either side of a membrane (see [106] for a description). While complete solutions of the Riemann problem are painstaking to evolve, a number of approximation techniques exist and do not reduce the order of accuracy of the code. Finally, one must take the conservative solution, now advanced forward in time, and recover the underlying primitive variables, a process that requires solving as many as eight simultaneous equations in the case of GRMHD or five for GRHD systems. A number of methods to do this have been widely studied [203], and simplifying techniques are known for specific cases (for the case of polytropic EOSs in GRHD evolutions, one need only invert a single non-analytic expression and the remaining variables can then be derived analytically).

The inclusion of magnetic fields in hydrodynamic calculations adds another layer of complexity beyond shock capturing. Magnetic fields must be evolved in such a way that they remain divergence-free, much in the same way that relativistic evolutions must satisfy the Hamiltonian and momentum constraints. Brute force attempts to subtract away any spurious divergence often lead to instabilities, so more intricate schemes have been developed. “Divergence cleaning” schemes typically introduce a new field representing the magnetic field divergence and use parabolic/hyperbolic equations to damp the divergence away while moving it off the computational domain; the approach is relatively simple to implement but prone to small-scale numerical errors [24]. “Constrained transport” schemes stagger the grids on which different physical terms are calculated to enforce the constraints (see, e.g., [310] for a particular implementation), and have been applied widely to many different physical configurations. Recently, a new scheme in which the vector potential is used rather than the magnetic field was introduced by Etienne, Liu, and Shapiro [93, 95], and found to yield successful results for a variety of physical configurations including NSs and BHs.

5.3 Microphysical numerical techniques

5.3.1 Neutron star physics and equations of state

One of the largest uncertainties in the input physics of NS-NS merger simulations is the true behavior of the nuclear matter EOS. To date, EM observations have yielded relatively weak constraints on the NS mass-radius relation, with the most precise simultaneous measurement of both as of now resulting from observations of Type 1 X-ray bursts from accreting NSs in three different sources [220]. In each case, the NS mass was found to lie in the range 1.3 M_⊙ ≲ M_NS ≲ 2 M_⊙ and the radius 8 km ≲ R_NS ≲ 12 km, implying a NS compactness

$${\mathcal C} \equiv {{G{M_{{\rm{NS}}}}} \over {{R_{{\rm{NS}}}}{c^2}}} = 0.1476\left({{{{M_{{\rm{NS}}}}} \over {{N_ \odot}}}} \right){\left({{{{R_{{\rm{NS}}}}} \over {10\,{\rm{km}}}}} \right)^{- 1}} \approx 0{.}16 - 0{.}37.$$

(28)

A more stringent constraint on the NS EOS is provided by observations of the Shapiro time delay in the binary millisecond pulsar PSR J1614-2230, which was found to have a mass M_NS = 1.97 ± 0.04 M_⊙ [81], which would rule out extremely soft EOS models incapable of supporting such a massive NS against collapse. As we discuss in more detail below, GW observations are likely to eventually yield tighter constraints than our current EM-based ones, though BH-NS mergers, which can undergo stronger tidal disruptions than NS-NS mergers at frequencies closer to LIGO and other GW observatories’ maximum frequency sensitivity band, may prove to be more useful for the task than NS-NS mergers.

Given the large theoretical uncertainties in describing the proper physical NS EOS, many groups have chosen the simplest possible parameterization: a polytrope (see Eq. 12). Under this choice, the enthalpy h takes the particularly simple form

$$h \equiv 1 + \varepsilon + P/\rho = 1 + \Gamma \varepsilon.$$

(29)

Initial data are generally assumed to follow the relation

$$P = K{\rho ^\Gamma},$$

(30)

where K is constant across the fluid. In the presence of shocks, the value of K for a particular fluid element will increase with time. We note that the Whisky group [17, 18] uses the term “polytropic” to refer to simulations in which Eq. 30 is enforced throughout, which implies adiabatic evolution without shock heating, and use the term “ideal fluid” to describe an EOS that includes the effects of shock heating and enforces Eq. 29.

Since the temperatures of NSs typically yield thermal energies per baryon substantially below the Fermi energy, one may treat nearly all NSs as effectively cold, except for the most recently born ones. During the merger process for NS-NS binaries, the matter will remain cold until the two NSs are tidally disrupted and a disk forms, at which point the thermal energy input and substantially reduced fluid densities require a temperature evolution model to properly model the underlying physics. In light of these results, some groups adopt a two-phase model for the NS EOS (see, e.g., [286]), where a cold, zero-temperature EOS, evaluated as a function of the density only, encodes as much information about as we possess about the NS EOS, and the hot phase depends on both the density and internal energy, typically in a polytropic way,

$$P(\rho, \varepsilon) = {P_{{\rm{cold}}}}(\rho) + {P_{{\rm{hot}}}}(\rho, \varepsilon)\quad [{P_{{\rm{hot}}}}(\rho, \varepsilon) = (\Gamma - 1)\rho (\varepsilon - {\varepsilon _{{\rm{cold}}}})].$$

(31)

There are a number of physically motivated EOS models that have been implemented for merger simulations, whose exact properties vary depending on the assumptions of the underlying model. These include models for which the pressure is tabulated as a function of the density only: FPS [222], SLy [83], and APR [3]; as well as models including a temperature dependence: Shen [268, 267] and Lattimer-Swesty [162]. A variety of models have been used to study the effects of quarks, kaons, and other condensates, which typically serve to soften the EOS, leading to reduced maximum masses and more compact NSs [223, 119, 230, 120, 23, 5].

Given the variance among even the physically motivated EOS models, it has proven useful to parameterize known EOS models with a much more restricted set of parameters. In a series of works, a Milwaukee/Tokyo collaboration determined that essentially all current EOS models could be fit using four parameters, so that their imprint on GW signal properties could be easily analyzed [237, 238, 184]. Their method assumes that the SLy EOS describes NS matter at low densities, and that the EOS at higher densities can be described by a piecewise polytropic fit with breaks at ρ = 10^14.7 and 10¹⁵ g/cm³. The four resulting parameters are P₁ = P (ρ = 10^14.7), the pressure at the first breakpoint density, which normalizes the overall density scale, as well as Γ₁, Γ₂, Γ₃, the adiabatic exponents in the three regions. Their results indicate that advanced LIGO should be able to determine the NS radius to approximately 1 km at an effective distance of 100 Mpc, which would place tight constraints on the value of P₁ in particular.

5.4 Electromagnetic and neutrino signature modeling

Motivated by the evidence that SGRBs frequently appear in galaxies with very low star formation rates [40, 109], astronomers have suggested that their progenitors are likely to be mergers of either NS-NS and/or BH-NS binaries. While soft-gamma repeaters (SGRs) have been confirmed as an SGRB source from observations of the system SGR 1806-20, they make up no more than approximately 15% of the total observed SGRB fraction according to the leading population estimates [164, 199]. There has been much interest in predicting the EM signatures of NS-NS and BH-NS mergers, along with the associated neutrino emission. The simplest models estimate a local radiation cooling rate for the matter but do not attempt to follow the paths of the photons and/or neutrinos after they are emitted, instead calculating the time-dependent luminosity assuming free streaming. Such models have been used in non-GR simulations of binary mergers going back more than a decade [253, 246], and recently such schemes have been used to perform full GR NS-NS mergers [265], including a self-consistent evolution of the electron fraction of the material Y_e, rather than a passive advection approach.

More complicated flux-limited diffusion schemes, in which the neutrino fluxes for given species and energies are given by explicit formulae that limit to the correct values for zero optical depth (free-streaming) and very large optical depth (diffusion), have been used as a post-processing tool to investigate the merger remnants in Newtonian NS-NS mergers [82], but have yet to be applied to full GR simulations. Finally, radiation transport schemes to evolve EM and neutrino fluxes passing through fluid configurations have been implemented in numerical GR codes [80, 103], but have yet to be used in binary merger simulations.

5.5 GW signal modeling

Measuring the GW signal from a dynamical merger calculation is a rather difficult task. One must determine, using a method unaffected by gauge effects, the perturbations at asymptotically large distances from a source by extrapolating various quantities measured at large but finite distances from the merger itself.

In the early days of numerical merger simulations, most groups typically assumed Newtonian and/or quasi-Newtonian gravitation, for which there is no well-defined dynamical spacetime metric. GW signals were typically calculated using the quadrupole formalism, which technically only applies for slow-moving, non-relativistic sources (see [193] for a thorough review of the theory). Temporarily reintroducing physical constants, the strains of the two polarizations for signals emitted in the z-direction are

$$\begin{array}{*{20}c} h_{+} = {G \over rc^{4}}(\ddot{\rlap{\!-}\bf{I}}_{xx} - \ddot {\rlap{\!-}\bf{I}}_{yy}), \\ h_ {\times} = {2G \over rc^{4}} \ddot{\rlap{\!-}\bf{I}}_{xy}, \qquad\quad \end{array}$$

where r is the distance from the source to the observer and ${\ddot{{\rlap{\!-}{\bf I}}}_{ij}}$ it the traceless quadrupole moment of the system. The energy and angular momentum loss rates of the system due to GW emission are given, respectively, by

$$\begin{array}{*{20}c} {- {{\left({{{dE} \over {dt}}} \right)}_{{\rm{GW}}}} = {G \over {5{c^5}}}{{\overset{\ldots}{{\rlap{\!-}{\mathbf I}}}}_{ij}}{{\overset{\ldots}{{\rlap{\!-}{\mathbf I}}}}_{ij}},}\quad\\ {- {{\left({{{d{L_k}} \over {dt}}} \right)}_{{\rm{GW}}}} = {{2G} \over {5{c^5}}}{\epsilon _{ijk}}{\ddot{{\rlap{\!-}\mathbf{I}}}_{il}}{{\overset{\ldots}{{\rlap{\!-}{\mathbf I}}}}_{lj}}.}\\ \end{array}$$

while only approximate, the quadrupole formulae do yield equations that are extremely straightforward to implement in both grid and particle-based codes using standard integration techniques.

Quadrupole methods were adopted for later PN and CF simulations, again because the metric was assumed either to be static or artificially constrained in such a way that made self-consistent determination of the GW signal impossible. One important development from this period was the introduction of a simple method to calculate the GW energy spectrum dE/df from the GW time-series through Fourier transforming into the frequency domain [330]. GW signals analyzed in the frequency domain allowed for direct comparison with the LIGO noise curve, making it much easier to determine approximate distances at which various GW sources would be detectable and the potential signal-to-noise ratio that would result from a template search. To constrain the nuclear matter EOS, one can examine where a GW merger spectrum deviates in a measurable way from the quadrupole point-mass form,

$${\left({{{dE} \over {df}}} \right)_{{\rm{GW}}}} = {{\pi G{m_1}{m_2}} \over 3}{(\pi G({m_1} + {m_2}){f_{{\rm{GW}}}})^{- 1/3}};\quad {f_{{\rm{GW}}}} \equiv 2{f_{{\rm{orb}}}} = {{{\omega _{{\rm{orb}}}}} \over \pi},$$

(32)

because of finite-size effects, and then link the deviation to the properties of the NS [98], as we show in Figure 9.

Figure 9

Approximate energy spectrum dE_GW/df derived from QE sequences of equal-mass NS-NS binaries with isolated ADM masses M_NS = 1.35 M_⊙ and a Γ = 2 EOS, but varying compactnesses (denoted M/R here), originally described in [303]. The diagonal lines show the energy spectrum corresponding to a point-mass binary, as well as values with 90%, 75%, and 50% of the power at a given frequency. Asterisks indicate the onset of mass-shedding, beyond which QE results are no longer valid. Image reproduced by permission from Figure 2 of [98], copyright by APS.

Full GR dynamical calculations, in which the metric is evolved according to the Einstein equations, generally use one of two approaches to calculate the GW signal from the merger, if not both. The first method, developed first by by Regge and Wheeler [239] and Zerilli [335] and written down in a gauge-invariant way by Moncrief [195] involves analyzing perturbations of the metric away from a Schwarzschild background. The second uses the Newman-Penrose formalism [202] to calculate the Weyl scalar ψ₄, a contraction of the Weyl curvature tensor, to represent the outgoing wave content on a specially constructed null tetrad that may be calculated approximately [60]. The two methods are complementary since they incorporate different metric information and require different numerical integrations to produce a GW time series. Regardless of the method used to calculate the GW signal, results are often presented by calculating the dominant s = −2 spin-weighted spherical harmonic mode. For circular binaries, the l = 2, m = 2 mode generally carries the most energy, followed by other harmonics; in cases where the components of the binary have nearly equal masses and the orbit is circular, the falloff is typically quite rapid, while extreme mass ratios can pump a significant amount of the total energy into other harmonics. For elliptical orbits, other modes can dominate the signal, e.g., a 3:1 ratio in power for the l = 2, m = 0 mode to the l = 2, m = 2 mode observed for high-ellipticity close orbits in [122]. A thorough summary of both methods and their implementation may be found in [257].

6 Dynamical Calculations

NS-NS merger simulations address a broad set of questions, which can be roughly summarized as follows (note the same questions apply to BH-NS mergers as well):

• 1. What is the final fate of the system, assuming a given set of initial parameters? Do we get a prompt collapse to a BH or the formation of a HMNS supported against collapse by differential rotation? Other outcomes are disfavored, at least for pre-merger NSs with masses M_NS ≳ 1.4 M_⊙ since the supramassive limit is at most 20% larger than that of a non-rotating NS [70, 71], and even for the stiffest EOS these values are typically less than 2.8 M_⊙.

• 2. What is the GW signal from the merger, and how does it inform us about the initial pre-merger parameters of the system?

• 3. What fraction of the system mass is left in a disk around the central BH or HMNS? While deriving exact EM emission profiles from a hydrodynamical configuration remains a challenge for the future, minimum conditions that would allow for the energy release observed in SGRBs have been established based on scaling arguments.

• 4. What is the neutrino and EM emission from the system, in both the time and energy domains? Obviously, the answer to this question and those that follow depend critically on the answers above.

• 5. What role do B-fields play in the GW, EM, and neutrino emission, and how does that tie in with other models suspected of having the same disk/jet geometry and gamma-ray emission like active galactic nuclei, pre-main sequence stars, etc.?

• 6. Do mergers produce a cosmologically significant quantity of r-process elements, or do those likely get produced by other astrophysical events instead?

The influence of the gravitational formalism used in a numerical simulation on the answer one finds for the questions above differs item by item. Determining the final fate of a merging system is highly dependent on the gravitational formalism; NS-NS merger remnants only undergo collapse in quasi-relativistic and fully GR schemes. Moreover, orbital dynamics at separations comparable to the ISCO and even somewhat larger depend strongly on the gravitational scheme. In particular, mass loss rates into a disk are often suppressed by orders of magnitude in GR calculations when compared to CF simulations, and even more so in comparison to PN and Newtonian calculations. EM emission profiles from a disk are difficult to calculate accurately without the use of full GR for this reason. On the other hand, while GR is required to calculate the exact GW signal from a merger, even early Newtonian simulations predicted many of the qualitative GW emission features correctly, and PN and CF schemes yielded results with some degree of quantitative accuracy about the full wavetrain. B-fields have only begun to be explored, but it already seems clear that they will affect the hydrodynamical evolution primarily after the merger in cases where differential rotation in a HMNS or disk winds up magnetic field strengths up to energy equipartition levels, vastly stronger than those found in pre-merger NSs. For such configurations, non-relativistic calculations can often reproduce the basic physical scenario but full GR is required to properly understand the underlying dynamics. Finally, the production of r-process elements, which depends sensitively on the thermodynamic evolution of the merger, seems to generally disfavor binary mergers as a significant source of the observed stellar abundances since the temperature and thus the electron fraction of the fluid remains too small [252], regardless of the nature of the gravitational treatment used in the calculations. This picture may need to be revised if significant mass loss occurs from the hot accretion disk that forms around the central post-merger object, possibly due to energy release from the r-process itself, but numerical calculations do not currently predict sufficient mass loss to match observations [297]. We will address each of these topics in greater detail in the sections below.

Since the first NS-NS merger calculations, there have been two main directions for improvements: more accurate relativistic gravitation, resulting in the current codes that operate using a self-consistent fully GR approach, and the addition of microphysical effects, which now include treatments of magnetic fields and neutrino/EM radiation. Noting that several of the following developments overlapped in time, e.g., the first full GR simulations by Shibata and Uryū [287] are coincident with the first PN SPH calculations, and predate the first CF SPH calculations, we consider in turn the original Newtonian calculations, those performed using approximate relativistic schemes, the calculations performed using full GR, and finally those that have included more advanced microphysical treatments.

6.1 Quasi-equilibrium and semi-analytic methods vs fully dynamical results

Before reviewing fully dynamical calculations of NS-NS mergers, it is worthwhile to ask how much information can already be deduced from QE calculations, which may be performed at much smaller computational cost, as well as from semi-analytic PN treatments and related approximate techniques. Clearly, the details of the merger and ringdown phases fall outside the QE regime, so only dynamical calculations can yield reliable information about the stability of remnants, properties of ejecta, or other processes that arise during the merger itself or in its aftermath. Thus, the primary point of comparison is the GW signal just prior to merger, which is also easier to detect (for first and second generation interferometers).

The strength of QE calculations lies in their ability to model self-consistently finite-size effects not captured in PN treatments (which always assume two orbiting point masses). The increased tidal interaction between the objects typically results in a more rapid phase advance of the binary orbit, which is important for constructing template waveforms that cover the entire NS-NS inspiral, merger, and ringdown. While QE sequences potentially offer a wealth of information about well-separated binaries and can help fix the phase evolution of the inspiraling binary, they do have two weaknesses arising as the binary approaches the stability limit. First, most QE methods, including the CTS formalism described in Section 4.2.1, are time-symmetric, and assume that the NS possess a symmetry plane perpendicular to the direction of motion (i.e., a front-back symmetry whose axis is perpendicular to the orbital angular momentum and the binary separation vector). In reality, tidal lags develop prior to final plunge, with the innermost edge of each NS rotating forward and the outer edge backwards. This effect has been captured in analytic and semi-analytic approaches (see, e.g., [159] for an early example), and is clearly seen in dynamical calculations (see Figure 3), but is not captured in CTS-based schemes (tidal lags also develop in BH-NS merger calculations when the BH has a non-zero spin, since this breaks the front-back symmetry; see [300] for an example).

A second weakness of QE methods is the treatment of the ISCO, particularly its importance as a characteristic point along an evolutionary sequence that, in theory, could encode information about the NS EOS. Simple estimates of the infall trajectory derived solely from QE sequences predict a very sudden and rapid infall near the ISCO, i.e., the point where the binding energy reaches a minimum along the sequence (see, e.g., the argument in [98]). However, this is clearly an oversimplification. In reality, binaries transition more gradually to the merger phase, and the inward plunge may occur significantly before reaching the formal ISCO; this in turns leads to more rapidly growing deviations from the QE approximation. Looking at the GW energy spectrum, one typically sees minor deviations from the point-mass predictions at frequencies below those characterizing the ISCO, but substantially more power at frequencies above it. Equivalently, the cutoff frequency for GW emission f_cut, where the spectrum starts deviating strongly from the point-mass prediction, is usually higher than the QE frequency near the ISCO, f_ISCO, while simple QE estimates assume these two frequencies to coincide.

To date, most attempts to generate waveforms extended back to arbitrarily early starting points involve numerically matching PN signals, typically generated using the Taylor T4 approach [53], onto the early stages of numerically generated waveforms, with some form of maximum overlap method used to provide the most physical transition from one to the other. These approaches may be improved by adding tidal effects to the evolution, typically parameterized by the tidal Love numbers that describe how tidal gravity fields induce quadrupole deformations [105]. Tidal effects can be placed into a relativistic framework [46, 74], which may be included within the effective one-body (EOB) formalism to produce high-accuracy waveforms [75]. In the EOB approach [58], resummation methods are used to include higher-order PN effects, though some otherwise unfixed parameters need to be set by comparing to numerical simulations.

Work is in its early stages to compare directly the GW spectra inferred from QE sequences of NS-NS binaries with those generated in numerical relativity simulations, but this comparison has been discussed at some length with regard to BH-NS mergers. Noting that NS-NS mergers generally correspond more closely to the BH-NS cases in which an ISCO is reached prior to the onset of tidal disruption, the KT collaboration [283, 276] concluded that the cutoff frequency marking significant deviations from PN point-mass behavior is roughly 30% higher than that marking emission near the classical ISCO for BH-NS systems (f_cut ≃ 1.3 f_ISCO).

A more detailed study has now been performed comparing EOB methods to numerical evolutions. By comparing to long-term simulations of NS-NS mergers, Baiotti et al. [15] find that EOB models may be tuned, via careful choices of their unfixed parameters, to reproduce the GW phases and amplitudes seen in NR evolutions up until the onset of merger. They further suggest that the EOB approach seems to cover a wider range of phase space than the Taylor T4 approach, presumable because of a more consistent representation of tidal effects, and offers the best route forward for construction of more accurate NS-NS inspiral templates.

6.2 Early dynamical calculations

The earliest NS-NS merger calculations were performed in Newtonian gravity, sometimes with the addition of lowest-order 2.5PN radiation reaction forces, and typically assumed that the NS EOS was polytropic. Both Eulerian grid codes [214, 196, 215, 197, 278, 255, 201, 298] and Lagrangian SPH codes [234, 235, 236, 76, 330, 331] were employed, and GW signals were derived under the assumptions of the quadrupole formalism. Configurations in Newtonian gravity cannot collapse, so a stable (possibly hypermassive) remnant was always formed. For polytropic EOS models with adiabatic indices larger than the classical minimum for production of a Jacobi ellipsoid, Γ ≳ 2.6, remnants were typically triaxial and maintained a significant-amplitude GW signal until the end of the simulation. For simulations using smaller values of Γ, remnants rapidly relaxed to spheroidal configurations, quickly damping away the resulting GW signal. Mass loss from the central remnant was often quite significant, with thick accretion disks or completely unbound material comprising up to 10–20% of the total system mass

Mass loss was suppressed in numerical simulations by constructing irrotational, rather than synchronized, initial data. Irrotational flow is widely thought to be the more physically realistic case, since viscous forces are much too weak to synchronize a NS prior to merger [45, 146]. When irrotational NSs (which are counter-rotating in the corotating frame of the binary) first make contact, a vortex sheet forms. Since the low-density fluid layers at the contact surface are surrounded at first contact by the denser fluid layers located originally within each NS, the configuration is well understood to be Kelvin-Helmholtz unstable, resulting in rapid mixing through vortex production. Meanwhile, mass loss through the outer Lagrange points is hampered by the reduced rotational velocity along the outer halves of each NS.

The GW emission from these mergers is composed of a “chirp,” increasing in frequency and amplitude as the NSs spiral inward, followed by a ringdown signal once the stars collide and merge. In [330, 331], a procedure to calculate the energy spectrum in the frequency band was laid out, with the resulting signal following the quadrupole, point-mass power-law form up to GW frequencies characterizing the beginning of the plunge. Above the plunge frequency, a sharp drop in the GW energy was seen, followed in some cases by spikes at kHz frequencies representing coherent emission during the ringdown phase.

6.3 Approximate relativistic schemes

The first steps toward approximating the effects of GR included the use of 1PN dynamics or the CF approximation. Using a formalism derived by Blanchet, Damour, and Schäfer [48], the 1PN equations of motion require the solution of eight Poisson-like equations in the form

$${\nabla ^2}\psi = f(\vec x),$$

(33)

where the source terms $f(\vec x)$ are compactly supported, and thus the fields ψ may be determined using the same techniques already in place to find the Newtonian potential. Adding in the lowest-order dissipative radiation reaction effects requires solution of a ninth Poisson equation for a reaction potential. The 1PN formalism was implemented in both grid-based [216] and SPH codes [10, 99, 101, 100]. Unfortunately, physically realistic NSs are difficult to model using a PN expansion, since the characteristic NS compactness ${\mathcal C}\gtrsim 0.15$, leads to first order “corrections” that often rival Newtonian terms in magnitude. To deal with this problem Ayal et al. [10] considered large (R ≈ 30 km), low-mass (< 1 M_⊙) NSs, allowing them to study relativistic effects but making results more difficult to interpret for physically realistic mergers. In [99, 101, 100], a dual speed of light approach was used, in which all 1PN effects were scaled down by a constant factor to yield smaller quantities while Newtonian and radiation reaction terms were included at full-strength. Both SPH groups found that the GW signal in PN mergers is strongly modulated, whereas Newtonian merger calculations typically yielded smooth, either monotonically decreasing or nearly constant-amplitude ringdown signals. Even reduced 1PN effects were shown to suppress mass loss by a factor of 2–5 for initially synchronized cases, and disk formation was seen to be virtually non-existent for initially irrotational, equal-mass NSs with a stiff (Γ = 3 polytropic) EOS [101].

Moving beyond the linearized regime, several groups explored the CF approximation, which incorporates many of the nonlinear effects of GR into an elliptic, rather than hyperbolic, evolution scheme. While nonlinear elliptic solvers are expensive computationally, they typically yield stable evolution schemes since field solutions are always calculated instantaneously from the given matter configuration. Summarized quickly, the CF approach involves solving the CTS field equations, Eqs. 18, 19, and 20, at every timestep, and evolving the matter configuration forward in time. The metric fields act like potentials, with various gradients appearing in the Euler and energy equations. While the CTS formalism remains the most widely used method to construct NS-NS (and BH-NS) initial data, it does not provide a completely consistent dynamical solution to the GR field equations. In particular, while it reproduces spherically symmetric configurations like the Schwarzschild solution exactly, it cannot describe more complicated configurations, including Kerr BHs. Moreover, because the CF approximation is time-symmetric, it also does not allow one to consistently predict the GW signal from a merging configuration. As a result, most dynamical calculations are performed by adding the lowest-order dissipative radiation reaction terms, either in the quadrupole limit or via the radiation reaction potential introduced in [48].

The CTS equations themselves were originally written down in essentially complete form by Isenberg in the 1970s, but his paper was rejected and only published after a delay of nearly 30 years [137]. In the intervening years, Wilson, Mathews, and Marronetti [327, 328, 188, 187] independently re-derived the entire formalism and used it to perform the first nonlinear calculations of NS-NS mergers (as a result, the formalism is often referred to as the “Wilson-Mathews” or “Isenberg-Wilson-Mathews” formalism). The key result in [327, 328, 188, 187] was the existence of a “collapse instability,” in which the deeper gravitational wells experienced by the NSs as they approached each other prior to merger could force one or both to collapse to BHs prior to the orbit itself becoming unstable. Unfortunately, their results were affected by an error, pointed out in [104], which meant that much of the observed compression was spurious. While their later calculations still found some increase in the central density as the NSs approached each other [189], these results have been contradicted by other QE sequence calculations (see, e.g., [313]). Furthermore, using a “CF-like” formalism in which the nonlinear source terms for the field equations are ignored, dynamical calculations demonstrated the maximum allowed mass for a NS actually increases in response to the growing tidal stress [273].

The CF approach was adapted into a Lagrangian scheme for SPH calculations by the same groups that had investigated PN NS-NS mergers, with Oechslin, Rosswog, and Thielemann [211] using a multigrid scheme and Faber, Grandclément, and Rasio [97] a spectral solver based on the Lorene libraries [124]. The effects of nonlinear gravity were immediately evident in both sets of calculations. In [211], NS-NS binaries consisting of initially synchronized NSs merged without appreciable mass loss, with no more that ∼ 10⁻⁴ of the total system mass ejected, strikingly different from previous Newtonian and PN simulations. When evolving initially irrotational systems, [97] found no appreciable developments of “spiral arms” whatsoever, indicating a complete lack of mass loss through the outer Lagrange points. Both groups also found strong emission from remnants for a stiff EOS, as the triaxial merger remnant produced an extended period of strong ringdown emission. Neither set of calculations indicated that the remnant should collapse promptly to form a BH, but given the high spin of the remnant it was noted that conformal flatness would have already broken down for those systems.

6.4 Full GR calculations

A summary of full GR calculations of NS-NS mergers is presented in Table 4. The KT collaboration was responsible for the only full GR calculations of NS-NS mergers that predate the breakthrough calculations of numerically stable binary BH evolutions [231, 22, 61], which have since transformed the field of GR hydrodynamics and MHD in addition to vacuum relativistic evolutions (Miller et al. [192] performed NS-NS inspiral calculations in full GR, but were not able to follow binaries through to merger). The first calculations of NS-NS mergers using a completely self-consistent treatment of GR were performed by Shibata and collaborators in the KT collaboration using a grid based code and the BSSN formalism [287]. CTS initial data consisting of equal-mass NSs described by a Γ = 2 polytropic EOS were constructed via SCF techniques [313], for both synchronized and irrotational configurations. The hyperbolic system was evolved on a grid, with an approximate maximal slicing condition that results in a parabolic equation for the lapse [271] and an approximate minimal distortion condition for the shift vector requiring the solution of an elliptic equation at every time step [272]. The shift vector gauge condition was found to fail when BHs were produced in the merger remnant, a well-known problem that had long bedeviled simulations involving binary and even single BH evolutions, so modifications were introduced to extend the stability of the algorithm as far as possible. Among the key results from this early work was a clear differentiation between mergers of moderately low-compactness NSs $({\mathcal C} \gtrsim 0.11)$, where the remnant collapsed promptly to a BH, and very low-compactness models, which yielded hypermassive remnants stabilized against gravitational collapse by differential rotation. Virtually all the NS matter was contained within the remnant for initially irrotational models, which served as evidence against equal-mass NSs mergers being a leading source of r-process elements in the universe through ejection. The lack of significant mass loss in equal-mass mergers, together with insignificant shock-heating of the material, also argued against the likelihood of such mergers as progenitors for SGRBs if the gamma-ray emission was assumed to be coincident with the GW burst; instead a delayed burst following the collapse of a HMNS to a BH appeared more likely.

Table 4 A summary of Full GR NS-NS merger calculations. EOS models include polytropes, piecewise polytropes (PP), as well as physically motivated models including cold SLy [83], FPS [222], and APR [3] models to which one adds an ideal-gas hot component to reflect shock heating, as well as the Shen [268, 267] finite temperature model and EOS that include Hyperonic contributions [264]. “Co/Ir” indicates that both corotating and irrotational models were considered; “BHB” indicates that BH binary mergers were also presented, including both BH-BH and BH-NS types, “ν-leak” indicates a neutrino leakage scheme was included in the calculation, “GH” indicates calculations were performed using the GHG formalism rather than BSSN, “non-QE” indicates superposition initial data were used, including cases where eccentric configurations were studied (“Eccen.”); “MHD” indicates MHD was used to evolve the system.

Later works, in particular a paper by Shibata, Taniguchi, and Uryū [285], introduced several new techniques to perform dynamical calculations that most codes at present still include in nearly the same or lightly modified form. These included the use of a high-resolution shock-capturing scheme for the hydrodynamics, as well as a Gamma-driver shift condition closely resembling the moving puncture gauge conditions that later proved instrumental in allowing for long-term BH evolution calculations. In the series of papers that followed their original calculations, the KT group established a number of results about NS-NS mergers that form the basis for much of our thinking about their hydrodynamic evolution:

• By varying the EOS model for the NS as well as the mass ratio, it was possible to constrain the binary parameters separating cases that form a HMNS rather than producing prompt collapse to a BH, and it was quickly determined that the total system mass as a proportion of the maximum allowed mass for an isolated NS is the key parameter, with only weak dependence on the binary mass ratio.

  For polytropic EOS models, the critical compactness values leading to prompt collapse for equal-mass binary mergers were found to be ${\mathcal C} = 0.14$ for Γ = 2 and ${\mathcal C} = 0.16$ for Γ = 2.25 [288]. As a rough rule, collapse occurred for polytropic EOS when the total system rest-mass was at least 1.7 M_max, where M_max is the maximum mass of an isolated non-rotating NS for the given EOS. For physically motivated EOS models [286], the critical mass was significantly smaller; indeed, the critical NS mass was found to be ∼ 1.35 M_max for the SLy EOS [83] (i.e., collapse for M_tot ≥ 2.7 M_⊙ with M_max = 2.04 M_⊙) and ∼ 1.39 M_max for the FPS EOS [222] (collapse for M_tot ≥ 2.5 M_⊙ with M_max = 1.8 M_⊙). This was not a complete surprise, since for the physically motivated EOS the NS radius is nearly independent of the mass across much of the parameter space, limiting the ability of the HMNS to expand in response to the extra mass absorbed during the merger.

• The mass ratio was found to play a critical role in the evolution of the remnant/disk configuration, since unequal-mass cases are better characterized as disruptions of the smaller secondary followed by its accretion onto the primary, rather than a true merger between the two NSs. Disk masses from full GR calculations [285] are generically smaller than those predicted from non-GR calculations. For polytropic EOS, disks contain approximately 4% of the total system mass for mass ratios q ≃ 0.8, varying roughly ∝ (1 − q) for a fixed total mass, with the disk mass decreasing for heavier binaries (and thus larger compactnesses) given the stronger gravity of the central remnant. Using the stiffer APR EOS [3], the dependence on the mass ratio was seen to be much steeper for a physical EOS than for polytropes, scaling like (1 − qp)^p, where p ≃ 3–4 [282].

• With respect to GW emission, it was determined in [288] that in low-mass cases in which a HMNS was formed, stiffer polytropic EOS models were able to support the development of a bar-mode instability, leading to transient spiral arm formation from the remnant and an extended period of strong GW emission, in the characteristic modulated form that results from differentially rotating ellipsoids (see, e.g., [160, 274, 259, 16, 183, 72]). Unequal-mass cases typically yielded one high-frequency peak at roughly GW ≈ 2 kHz corresponding to non-axisymmetric oscillations, and equal-mass cases yielded multiple peaks including those associated with quasi-radial oscillations as well [285]. For physical EOS models [286], mass loss into a disk is reduced relative to the polytropic case given the higher compactness of the central region, and GW oscillation peaks, while very strong, occur at correspondingly higher frequencies. By contrast, prompt formation of a BH led to a ringdown signal with rapidly decreasing amplitude becoming negligible within a few dynamical times.

  The GW signals were evaluated under the gauge-dependent assumption of transverse tracelessness, and energy and angular momentum loss rates into each spherical harmonic mode were computed using the gauge-invariant Zerilli-Moncrief formalism [239, 335, 195] in much the same way that is used by some groups in numerical relativity today (many BH-BH and hydrodynamics simulations report GW signals derived from the alternate ψ₄ Weyl scalar formulation [202, 60], or use both methods).

• In [286], it was concluded that mergers of NSs with comparable masses made poor SGRB progenitor candidates, assuming prompt emission (because of the lack of energy available for neutrino annihilation), but that the energy budget in the HMNS case is orders of magnitude larger. Remarkably, this discussion from 2005 predates the first identifications of SGRBs with older populations, which greatly improved our theoretical understanding of compact object mergers as their likely progenitors. In the first work that followed the initial localizations of SGRBs, mergers of binaries with relatively small mass ratios, q ≈ 0.7, were seen to form sufficiently hot and massive disk to power a SGRB, albeit a relatively brief, low-luminosity one. It was suggested that the more likely SGRB progenitor is indeed a HMNS, since dissipative effects within the remnant can boost temperatures up to ∼ 10¹¹ K.

  Further approximate relativistic investigations of NS-NS mergers, along with BH-NS mergers, as potential SGRB sources quickly swept through the community after the initial localizations of SGRBs, with several groups using a wide variety of methods all concluding that mergers were plausible progenitors, but finding it extremely difficult to constrain the scenario in quantitative ways given the extremely complicated microphysics ultimately responsible for powering the burst (see, e.g., [207, 206] who investigated potential disk energies; [243], who modeled the fallback accretion phase onto a BH; and [266], who considered the thermodynamic and nuclear evolution of disks around newly-formed BHs produced by mergers). We will return to this topic below in light of recent GRMHD Simulations.

In the past few years, five groups have reported results from NS-NS mergers in full GR; KT, HAD, Whisky, UIUC, and Jena. Much of the work of the HAD and Whisky groups, developers respectively of the code of those names, began at Louisiana State University (HAD) and the Albert Einstein Institute in Potsdam (Whisky), though both efforts now include several other collaborating institutions. Two other groups, the SXS collaboration that originated at Caltech and Cornell, and the Princeton group, have reported BH-NS merger results and are actively studying NS-NS mergers as well, but have yet to publish their initial papers about the latter. All of the current groups use AMR-based Eulerian grid codes, with four evolving Einstein’s equations using the BSSN formalism and the HAD collaboration making use of the GHG method instead. HAD, Whisky, and UIUC have all reported results about magnetized NS-NS mergers (the KT collaboration has used a GRMHD code to study the evolution of magnetized HMNS, but not complete NS-NS mergers). The KT collaboration has considered a wide range of EOS models, including finite-temperature physical models such as the Shen EOS, and have also implemented a neutrino leakage scheme, while all other results reported to date have assumed a Γ = 2 polytropic EOS model.

Given the similarities of the various codes used to study NS-NS mergers, it is worthwhile to ask whether they do produce consistent results. A comparison paper between the Whisky code and the KT collaboration’s SACRA codes [20] found that both codes performed well for conservative global quantities, with global extrema such as the maximum rest-mass density in agreement to within 1% and waveform amplitudes and frequencies differing by no more than 10% throughout a full simulation, and typically much less.

Several of the the groups listed above have also been leaders in the field of BH-NS simulations: the KT, HAD, and UIUC groups have all presented BH-NS merger results, as have the SXS collaboration [85, 84, 108], and Princeton group [294, 88] (see [284] for a thorough review).

We discuss the current understanding of NS-NS mergers in light of all these calculations below.

6.4.1 HMNS and BH remnant properties

Using their newly developed SACRA code [332], the KT group [144], found that when a hybrid EOS is used to model the NS, in which the cold part is described by the APR EOS and the thermal component as a Γ = 2 ideal gas, the critical total binary mass for prompt collapse to a BH is M_tot = 2.8–2.9 M_⊙, independent of the initial binary mass ratio, a result consistent with previous explorations of other polytropic and physically motivated NS EOS models (see above). In all cases, the BH was formed with a spin parameter a ≈ 0.78 depending very weakly on the total system mass and mass ratio.

They further classified the critical masses for a number of other physical EOS in [134], finding that binaries with total masses M_tot ≲ 2.7 M_⊙ should yield long-lived HMNSs (> 10 ms) and substantial disk masses with M_disk > 0.04 M_⊙ assuming that the current limit on the heaviest observed NS, M = 1.97 M_⊙ [81] is correct. In Figure 10, we show the final fate of the merger remnant as a function of the total pre-merger mass of the binary. “Type I” indicates a prompt collapse of the merger remnant to a BH, “Type II” a short-lived HMNS, which lasts for less than 5 ms after the merger until its collapse, and “Type III” a long-lived HMNS which survives for at least 5 ms. See [134] for an explanation of the EOS used in each simulation.

Figure 10

Type of final remnant corresponding to different EOS models. The vertical axis shows the total mass of two NSs. The horizontal axis shows the EOSs together with the corresponding NS radii for M_NS = 1.4 M_⊙. Image reproduced by permission from Figure 3 of [134], copyright by APS.

While all of the above results incorporated shock heating, the addition of both finite-temperature effects in the EOS and neutrino emission modifies the numerically determined critical masses separating HMNS formation from prompt collapse. Adding in a neutrino leakage scheme for a NS-NS merger performed using the relatively stiff finite-temperature Shen EOS, the KT collaboration reports in [265] that HMNSs will form generically for binary masses ≲ 3.2 M_⊙, not because they are centrifugally-supported but rather because they are pressure-supported, with a remnant temperature in the range 30–70 MeV. Since they are not supported by differential rotation, these HMNSs were predicted to be stable until neutrino cooling, with luminosities of ∼ 3–10 × 10⁵³ erg/s, can remove the pressure support. Even for cases where the physical effects of hyperons were included, which effectively soften the EOS and reduce the maximum allowed mass for an isolated NS to 1.8 M_⊙, the KT collaboration [264] still finds that thermal support can stabilize HMNS with masses up to 2.7 M_⊙.

Using a Carpet/Cactus-based hydrodynamics code called Whisky [19] that works within the BSSN formalism (a version of which has been publicly released as GRHydro within the Einstein [90]), the Whisky collaboration has analyzed the dependence of disk masses on binary parameters in some detail. For mass ratios q = 0.7–1.0 [240], they found that bound disks with masses of up to 0.2 M_⊙ can be formed, with the disk mass following the approximate form

$${M_{{\rm{disk}}}} = 0{.}039({M_{\max}} - {M_{{\rm{tot}}}}) + 1{.}115(1 - q)({M_{\max}} - {M_{{\rm{tot}}}});\quad {M_{\max}} = 1{.}139(1 + q){M_\ast},$$

(34)

where M_max the maximum mass of a binary system for a given EOS (Γ = 2 ideal gas for these calculations), M_* is the maximum mass of an isolated non-rotating NS for the EOS, and M_tot the mass of the binary, with all masses here defined as baryonic. The evolution of the total rest mass present in the computational domain for a number of simulations is shown in Figure 11.

Figure 11

Evolution of the total rest mass M_tot of the remnant disk (outside the BH horizon) normalized to the initial value for NS-NS mergers using a Γ = 2 polytropic EOS with differing mass ratios and total masses. The order of magnitude of the mass fraction in the disk can be read off the logarithmic mass scale on the vertical axis. The curves referring to different models have been shifted in time to coincide at t_coll. Image reproduced by permisison from Figure 5 of [240], copyright by IOP.

6.4.2 Magnetized NS-NS mergers

Using the HAD code described in [8] that evolves the GHG system on an AMR-based grid with CENO reconstruction techniques, Anderson et al. [7] performed the first study of magnetic effects in full GR NS-NS mergers [6]. Beginning from spherical NSs with extremely strong poloidal magnetic fields (9.6 × 10¹⁵ G, as is found in magnetars), their merger simulations showed that magnetic repulsion can delay merger by 1–2 orbits and lead to the formation of magnetically buoyant cavities at the trailing end of each NS as contact is made (see Figure 12), although the latter may be affected by the non-equilibrium initial data. Both effects would have been greatly reduced if more realistic magnetic fields strengths had been considered. Magnetic fields in the HMNS remnant, which can be amplified through dynamo effects regardless of their initial strengths, helped to distribute angular momentum outward via the magnetorotational instability (MRI), leading to a less differentially rotating velocity profile and a more axisymmetric remnant. The GW emission in the magnetized case was seen to occur at lower characteristic frequencies and amplitudes as a result.

Figure 12

Fluid density isocontours and magnetic field distribution (in a plane slightly above the equator) immediately after first contact for a magnetized merger simulation. The cavities at both trailing edges are attributed to magnetic pressure inducing buoyancy. Image reproduced by permission from Figure 1 of [6], copyright by APS.

The UIUC group was among the first to produce fully self-consistent GRMHD results [87]. Using a newly developed Cactus-based code, they performed the first studies of unequal-mass magnetized NS-NS mergers [172]. Using poloidal, magnetar-level initial magnetic fields, Liu et al. found that magnetic effects are essentially negligible prior to merger, but can increase the mass in a disk around a newly formed BH moderately, from 1.3% to 1.8% of the total system mass for mass ratios of q = 0.85 and Γ = 2. They point out that MHD effects can efficiently channel outflows away from the system’s center after collapse [295], and may be important for the late-stage evolution of the system.

In [116], the Whisky group performed simulations of magnetized mergers with field strengths ranging from 10¹² to 10¹⁷ G. Agreeing with the UIUC work that magnetic field strengths would have essentially no effect on the GW emission during inspiral, they note that magnetic effects become significant for the HMNS, since differential rotation can amplify B-fields, with marked deviations in the GW spectrum appearing at frequencies of _GW ≳ 2 kHz. They also point out that high-order MHD reconstruction schemes, such as third-order PPM, can produce significantly more accurate results that second-order limiter-based schemes. A follow-up paper [117] showed that a plausible way to detect the effect of physically realistic magnetic fields on the GW signal from a merger was through a significant shortening of the timescale for a HMNS to collapse, though a third-generation GW detector could perhaps observe differences in the kHz emission of the HMNS as well.

More recently, they have used very long-term simulations to focus attention on the magnetic field strength and geometry found after the remnant collapses to a BH [241]. They find that the large, turbulent magnetic fields (B ∼ 10¹² G) present in the initial binary configuration are boosted exponentially in time up to a poloidal field of strength 10¹⁵ G in the remnant disk, with the field lines maintaining a half-opening angle of 30° along the BH spin axis, a configuration thought to be extremely promising for producing a SGRB. The resulting evolution, shown in Figure 13, is perhaps the most definitive result indicating that NS-NS mergers should produce SGRBs for some plausible range of initial parameters.

Figure 13

Evolution of the density in a NS-NS merger, with magnetic field lines superposed. The first panel shows the binary shortly after contact, while the second shows the short-lived HMNS remnant shortly before it collapses. In the latter two panels, a BH has already formed, and the disk around it winds up the magnetic field to a poloidal geometry of extremely large strength, ∼ 10¹⁵ G, with an half-opening angle of 30°, consistent with theoretical SGRB models. Image reproduced by permission from Figure 1 of [241], copyright by AAS.

It is worth noting that all magnetized NS-NS merger calculations that have been attempted to date have made use of unphysically large magnetic fields. This is not merely a convenience designed to enhance the role of magnetic effects during the merger, though it does have that effect. Rather, magnetic fields are boosted in HMNS remnants by the MRI, whose fastest growing unstable mode depends roughly linearly on the Alfvén speed, and thus the magnetic field strength. In order to move to physically reasonable magnetic field values, one would have to resolve the HMNS at least a factor of 100 times better in each of three dimensions, which is beyond the capability of even the largest supercomputers at present, and likely will be for some time to come.

6.4.3 GW emission

In [145], the KT collaboration found a nearly linear relationship between the GW spectrum cutoff frequency f_cut and the NS compactness, independent of the EOS, as well as a relationship between the disk mass and the width of the kHz hump seen in the GW energy spectrum. While f_cut is a somewhat crude measure of the NS compactness, it occurs at substantially lower frequencies than any emission process associated with merger remnants, and thus is the parameter most likely to be accessible to GW observations with a second generation detector.

The qualitative form of the high-frequency components of the GW spectrum is primarily determined by the type of remnant formed. In Figures 14 and 15, we show h(t) and $\tilde h(f)$, respectively, for four of the runs calculated by the KT collaboration and described in [134]. Type I collapses are characterized by a rapid decrease in the GW amplitude immediately after the merger, yielding relatively low power at frequencies above the cutoff frequency. Type II and III mergers yield longer periods of GW emission after the merger, especially the latter, with the remnant oscillation modes leading to clear peaks at GW frequencies f_GW = 2–4 khz that should someday be detectable by third generation detectors like the Einstein Telescope, or possibly even by advanced LIGO should the source be sufficiently close (D ≲ 20 Mpc) and the high-frequency peak of sufficiently high quality [265].

Figure 14

Dimensionless GW strain Dh/m₀, where D is the distance to the source and m₀ the total mass of the binary, versus time for four different NS-NS merger calculations. The different merger types become apparent in the post-merger GW signal, clearly indicating how BH formation rapidly drives the GW signal down to negligible amplitudes. Image reproduced by permission from Figures 5 and 6 of [134], copyright by APS.

Figure 15

Effective strain at a distance of 100 Mpc shown as a function of the GW frequency (solid red curve) for the same four merger calculations depicted in Figure 14. Post-merger quasi-periodic oscillations are seen as broad peaks in the GW spectrum at frequencies f_GW = 2–4 kHz. The blue curve shows the Taylor T4 result, which represents a particular method of deducing the signal from a 3PN evolution. The thick green dashed curve and orange dot-dashed curves depict the sensitivities of the second-generation Advanced LIGO and LCGT (Large Scale Cryogenic Gravitational Wave Telescope) detectors, respectively, while the maroon dashed curve shows the sensitivity of a hypothetical third-generation Einstein Telescope. Image reproduced by permission from Figures 5 and 6 of [134], copyright by APS.

Using new multi-orbit simulations of NS-NS mergers, Baiotti et al. [14, 15] showed that the semi-analytic effective one-body (EOB) formalism severely underestimates high-order relativistic corrections even when lowest-order finite-size tidal effects were included. As a result, phase errors of almost a quarter of a radian can develop, although these may be virtually eliminated by introducing a second-order “next-to-next-to-leading order” (NNLO) correction term and fixing the coefficient to match numerical results. The excellent agreement between pre-merger numerical waveforms and the revised semi-analytic EOB approximant is shown in Figure 16.

Figure 16

Comparison between numerical waveforms, shown as a solid black line, and semi-analytic NNLO EOB waveforms, shown as a red dashed line (top panel). The top panels show the real parts of the EOB and numerical relativity waveforms, and the middle panels display the corresponding phase differences between waveforms generated with the two methods. There is excellent agreement between with the numerical waveform almost up to the time of the merger as shown by the match of the orbital frequencies (bottom panel). Image reproduced by permission from Figure 14 of [15], copyright by APS.

6.4.4 Binary eccentricity

The effects of binary eccentricity on NS-NS mergers was recently studied by the Jena group [122]. Such systems, which would indicate dynamical formation processes rather than the long-term evolution of primordial binaries, evolve differently in several fundamental ways from binaries that merge from circular orbits. For nearly head-on collisions, they found prompt BH formation and negligible disk mass production, with only a single GW burst at frequencies comparable to the quasi-normal mode of the newly formed BH. For a collision in which mass transfer occurred at the first passage but two orbits were required to complete the merger and form a BH, a massive disk was formed, containing 8% of the total system mass even at time Δt = 100 M_tot ≈ 280 M_⊙ after the formation of the BH. During that time, the black hole accreted an even larger amount of mass, representing over twice the mass of the remaining disk. Between the first close passage and the second, during which the two NS merged, the GW signal was seen to be quasi-periodic, and a a frequency comparable to the fundamental oscillation mode of the two NS, a result that was duplicated in a calculation for which the periastron fell outside the Roche limit and the eccentric binary survived for the full duration of the run, comprising several orbits.

6.5 Simulations including microphysics

In parallel to efforts in full GR, there has also been great progress in numerical simulations that include approximate relativistic treatments but a more detailed approach to microphysical issues. The first simulations to use a realistic EOS for NS-NS mergers were performed by Ruffert, Janka, and collaborators [253, 139, 254], who assumed the Lattimer-Swesty EOS for their Newtonian PPM-based Eulerian calculations. They were able to determine a physically meaningful temperature for NS-NS merger remnants of 30–50 MeV, an overall neutrino luminosity of roughly 10⁵³ erg/s for tens of milliseconds, and a corresponding annihilation rate of 2–5 × 10⁵⁰ erg/s given the computed annihilation efficiencies of a few parts in a thousand. This resulted in an energy loss of 2–4 × 10⁴⁹ erg over the lifetime of the remnant [139], a value later confirmed in multigrid simulations that replaced the newly formed HMNS by a Newtonian or quasi-relativistic BH surrounded by the bound material making up a disk [251]. The temperatures in the resulting neutron-rich (Y_e ≈ 0.05–0.2) remnant were thought to be encouraging for the production of r-process elements [254], although numerical resolution of the low-density ejecta limited the ability to make quantitatively accurate estimates of its exact chemical distribution. Further calculations, some of which involved unequal-mass binaries, indicated that the temperatures and electron fractions in the ejecta were likely not sufficient to produce solar abundances of r-process elements [252], with electron fractions in particular smaller than those set by hand in the r-process production model that appears in [111, 245]. More recently, it was suggested [123] that the decompression of matter originally located in the inner crust of a NS and ejected during a merger has a nearly solar elemental distribution for heavy r-process elements (A > 140). This indicates that NS-NS mergers may be the source of the observed cosmic r-process elements should there be sufficient mass loss per merger event, M_ej ∼ 3–5 × 10⁻⁵ M_⊙, although these amounts have yet to be observed in full GR simulations which have often admittedly been performed using cruder microphysical treatments.

In [244], Rosswog and Davies included a detailed neutrino leakage scheme in their calculations and also adopted the Shen EOS for several calculations, finding in a later paper [248] that the gamma-ray energy release is roughly 10⁴⁸ erg, in line with previous results from other groups, but noting that the values would be significantly higher if temperatures in the remnant were higher, since the neutrino luminosity scales like a very high power of the temperature. These calculations also identified NS-NS mergers as likely SGRB candidates given the favorable geometry [249], and the possibility that the MRI in a HMNS remnant could dramatically boost magnetic fields on the sub-second timescales characterizing SGRBs [250]. Rosswog and Liebendörfer [246] found that electron antineutrinos ${{\bar \nu}_e}$ dominate the emission, as had Ruffert and Janka [252], though the exact thermodynamic and nuclear profiles were found to be somewhat sensitive to the properties of the EOS model. More recently, using the VULCAN 2-dimensional multi-group flux-limited-diffusion radiation hydrodynamics code [173] to evaluate slices taken from SPH calculations, Dessart et al. [82] found that neutrino heating of the remnant material will eject roug hly 10⁻⁴ M_⊙ from the system.

Price and Rosswog [233, 247] performed the first MHD simulation of merging NS-NS binaries using an SPH code that included magnetic field effects, finding that the Kelvin-Helmholtz unstable vortices formed at the contact surface between the two NSs could boost magnetic fields rapidly up to ∼ 10¹⁷ G. These results were not seen in GRMHD simulations, where gains in the magnetic field strength generated by dynamos were limited by the swamping of the vortex sheet at the surface of contact by rapidly infalling NS material that went on to form the eventual HMNS or BH [117]. Longer-term simulations did note that shearing instabilities were able to support power-law, or perhaps even exponential, growth of the magnetic fields on longer timescales (∼ 10 s of ms), which augurs well for NS-NS mergers as the central engines of SGRBs [241].

An effort to identify potential observational differences between NSs and COs with quark-matter interiors has been led by Oechslin and collaborators. Using an SPH code with CF gravity, Oechslin et al. [210, 212] considered mergers of NSs with quark cores described by the MIT bag model [67, 102, 142], which have significantly smaller maximum masses than traditional NSs. They found the hybrid nuclear-quark EOS yielded higher ISCO frequencies for NSs with masses ≳ 1.5 M_⊙ and slightly larger GW oscillation frequencies for any resulting merger remnant compared to purely hadronic EOS. Bauswein et al. [34] followed up this work by investigating whether “strangelets”, or small lumps of strange quark matter, would be ejected in sufficient amounts throughout the interstellar medium to begin the phase transition that would convert traditional hadronic NSs into strange stars. They determined that the total rate of strange matter ejection in NS-NS mergers could be as much as 10⁻⁸ M_⊙ per year per galaxy or essentially zero depending on the parameters input into the MIT bag model, with the upper values clearly detectable by orbiting magnetic spectrometers such as the AMS-02 detector that was recently installed on the International Space Station [182, 148]. Further calculations concluded that the mergers of strange stars produce a much more tenuous halo than traditional NS mergers, more rapid formation of a BH, and higher frequency ringdown emission [35], as we show in Figure 17.

Figure 17

Evolution of a binary strange star merger performed using a CF SPH evolution. The “spiral arms” representing mass loss through the outer Lagrange points of the system are substantially narrower than those typically seen in CF calculations of NS-NS mergers with typical nuclear EOS models. Image reproduced by permission from Figure 4 of [35], copyright by APS.

Oechslin, Janka, and Marek also analyzed a wide range of EOS models using their CF SPH code, finding that matter in spiral arms was typically cold and that the dynamics of the disk formed around a post-merger BH depends on the initial temperature assumed for the pre-merger NS [209]. They also determined that the kHz GW emission peaks produced by HMNSs may help to constrain various parameters of the original NS EOS, especially its high-density behavior [208], with further updates to the prediction provided by Bauswein and Janka [32]. Most recently, Stergioulas et al. [296] studied the effect of nonlinear mode couplings in HMNS oscillations, leading to the prediction of a triplet peak of frequencies being present or low mass (M_NS = 1.2–1.35 M_⊙) systems in the kHz range.

6.6 Comparison to BH-NS merger results

While there is a history of Newtonian, quasi-relativistic, post-Newtonian, and CF gravitational formalisms being used to perform BH-NS merger simulations, their results are nearly always quantitatively, if not qualitatively, different than full GR simulations, and we focus here on the latter (see [284] for a more thorough historical review). Most of the groups that have performed full GR NS-NS merger calculations have also published results on BH-NS mergers, including Whisky (for head-on collisions) [174], KT [290, 289, 283, 332, 276, 154, 153], UIUC [91, 94, 92], HAD [66], as well as the SXS collaboration [85, 84, 108, 107] and Princeton (for elliptical mergers) [294, 88]. Summarizing the results of these works, we get a rather coherent picture, which we describe below.

The GW signal from BH-NS mergers is somewhat “cleaner” than that from NS-NS mergers, since the disruption of the NS and its accretion by the BH rapidly terminate the GW emission. In general, 3PN estimates model the signal well until tidal effects become important. The more compact the NS, the higher the dimensionless “cutoff frequency” M_totf_cut at which the GW energy spectrum plummets, with direct plunges in which the NS is swallowed whole typically yielding excess power near the frequency maximum from the final pre-merger burst. For increasingly prograde BH spins, there is more excess power over the 3PN prediction at lower frequencies, but also a lower cutoff frequency marking the plunge (see the discussion in [284]). From an observational standpoint, the deviations from point-mass form become more visible for a higher mass BH-NS system, because frequencies scale characteristically like the inverse of the total mass. The distinction is particularly important for Advanced LIGO, as systems with BH ≳ 3 M_⊙ typically yield cutoff frequencies within the advanced LIGO band at source distances of D ∼ 100 Mpc, while for lower-mass systems the cutoff occurs at or just above the upper end of the frequency band. This is significantly different than the situation for NS-NS mergers, in which the characteristic frequencies corresponding to the merger itself typically fall at frequencies above the advanced LIGO high-frequency sensitivity limit, and those corresponding to remnant oscillations in the range 2–4 kHz, which will prove a challenge even for third-generation GW detectors.

Disk masses for BH-NS mergers were found to be extremely small in the first calculations, all performed using non-spinning BHs [290, 289, 91], but have since been corrected to larger values once more sophisticated grid-based schemes and atmosphere treatments were added to those codes. More recent results indicate disk masses for reasonable physical parameters can be as large as 0. 4 M_⊙, for highly-spinning (a_BH/M = 0.9) mergers [108], with values of 0.035–0.05 M_⊙ characterizing non-spinning models with mass ratios q ≈ 1/5 [153]. Mass loss into a disk is suppressed by misaligned spins, especially for highly-inclined BHs, so the aligned cases should currently be interpreted as upper limits for the disk mass when alignment is varied [108]. Overall, disk masses for BH-NS merger remnants are comparable to those from NS-NS merger remnants, and may not be clearly distinguishable from them based solely on the emission properties of the disk. For BH-NS mergers with mass ratios q =1/3 and prograde spins of dimensionless magnitude a_BH/M = 0.5, the disk parameters found after a run performed with the inclusion of a finite-temperature NS EOS [84] indicated that the neutrino luminosity from the disk might be as high as 10⁵³ erg/s. While NS-NS merger simulations have led to predictions of neutrino luminosities a few times larger than this, the result does indicate that BH-NS mergers are also plausible SGRB progenitor candidates, possibly with lower characteristic luminosities than bursts resulting from NS-NS mergers.

The role of magnetic fields in BH-NS mergers has only been investigated recently [66, 92], in simulations that apply an initially poloidal magnetic field to the NSs in the binary. Magnetic fields were found to have very little effect on the resulting GW signal and the mass accretion rate for the BH for physically reasonable magnetic field strengths, with visible divergences appearing only for B ∼ 10¹⁷ G [92], which is not particularly surprising. Just as in NS-NS mergers, magnetic fields play very little role during inspiral, and unlike the case of NS-NS mergers there is no opportunity to boost fields at a vortex sheet that forms when the binary makes contact, nor in a HMNS via differential rotation. While the MRI may be important in determining the thermal evolution and mass accretion rate in a post-merger disk, such effects will likely be observable primarily on longer timescales.

Just as full-GR NS-NS simulations do not indicate that such mergers are likely sources of the r-process elements we observe in the universe, BH-NS simulations in full GR make the same prediction: no detectable mass loss from the system whatsoever, at least in the calculations performed to date. The picture may change when even larger prograde spins are modeled, since this should lead to maximal disk production, or if more detailed microphysical treatments indicate that a significant wind can be generated from either a HMNS or BH disk and unbind astrophysically interesting amounts of material, but neither has been seen in the numerical results to date.

As is seen in NS-NS mergers, the pericenter distance plays a critical role in the evolution of eccentric BH-NS mergers as well. Large disk masses containing up to 0.3 M_⊙, with an unbound fraction of roughly 0.15 M_⊙, can occur when the periastron separation is located just outside the classical ISCO, with GW signals taking on the characteristic zoom-whirl form predicted for elliptical orbits [294]. In between pericenter passages, radial oscillations of the neutron star produce GW emission at frequencies corresponding to the f-mode for the NS as well [88].

7 Summary and Likely Future Directions

Returning to the questions posed in Section 6, we can now provide the current state of the field’s best answers, though this remains a very active area of research and new results will certainly continue to modify this picture.

• 1. With regard to the final fate of the merger remnant, calculations using full GR are required, but the details of the microphysics do not seem to play a very strong role. It is now possible to determine whether or not a pair of NSs with given parameters and specified EOS will form a BH or HMNS promptly after merger, and to estimate whether a HMNS will collapse on a dynamical timescale or one of the longer dissipative timescales (see, e.g., [134]). For NS-NS binaries with sufficiently small masses, it is also possible to determine quickly whether the remnant mass is below the supramassive limit for which a NS is stabilized against collapse by uniform rotation alone, and thus would be unlikely to collapse, barring a significant amount of fallback accretion, unless pulsar emission or magnetic field coupling to the outer disk reduced the rotation rate below the critical value. This scenario likely applies only for mergers where the total system mass is relatively small M_tot ≲ 2.5–2.6 M_⊙ [70], even taking into account the current maximum observed NS mass of M_NS = 1.97 M_⊙ [81]. Based on the wide arrays of EOS models already considered, it is entirely possible to infer the likely fate for sets of parameters and/or EOS models that have not yet been simulated, although no one has yet published a “master equation” that summarizes all of the current work into a single global form. While magnetic fields with realistic magnitudes are unlikely to affect the BH versus HMNS question [172, 117], finite-temperature effects might play a nontrivial role should NSs be sufficiently hot prior to merger [265] (and see also [209]). In the end, by the time the second generation of GW detectors make the first observations of mergers, the high-frequency shot-noise cutoff will prove to be a bigger obstacle to determining the fate of the remnant than any numerical uncertainty. A schematic diagram showing the possible final fates for a NS-NS merger along with the potential EM emission (see Figure 21 of [282]) is shown in Figure 18.

• 2. GW emission during merger is also well-understood, though there are a few gaps that need to be filled, with full GR again a vital requirement. While the PN inspiral signal prior to merger is very well understood, finite-size tidal effects introduce complications beyond those seen in BH-BH mergers, yet the longest calculations performed to date [15] encompass fewer orbits prior to merger than the longest BH-BH runs [261]. As noted in [15] and elsewhere, longer calculations will likely appear over time, helping to refine the prediction for the NS-NS merger GW signal as the binary transitions from a PN phase into one that can only be simulated using full GR, and teasing out the NS physics encoded in the GW signal. It seems clear from the published work that the emission during the onset of the merger is well-understood, as is the very rapid decay that occurs once the remnant collapses to a BH, either promptly or following some delay. GW emission from HMNSs has been investigated widely, and there have been correlations established between properties of the initial binary and the late-stage high-frequency emission (see, e.g., [145, 117]), but given that magnetic fields, neutrino cooling, and other microphysical effects seem to be important, a great deal of work remains to be done. Perhaps more importantly, since HMNSs emit radiation at frequencies well beyond the shot-noise limit of even second-generation GW detectors, while the final inspiral occurs near peak sensitivity, it is likely that the first observations of NSs will constrain the nuclear EOS (or perhaps the quark matter EOS [35]) primarily via the detection of small finite-size effects during inspiral. Since QE calculations are computationally inexpensive compared to numerical merger simulations, there should be much more numerical data available about the inspiral stage than other phases of NS-NS mergers, which should help optimize the inferences to be drawn from future observations.

• 3. Determining the mass of the thick disk that forms around a NS-NS merger remnant remains a very difficult challenge, since its density is much lower and harder to resolve using either grid-based or particle-based simulations. The parameterization given by Eq. 34 [240] is generally consistent with the GR calculations of other groups (see, e.g., [134]), and seems to reflect a current consensus. It is also clear that disk masses around HMNSs (up to 0.2 M_⊙) are significantly larger than those forming around prompt collapses, which are limited to about 0.05 M_⊙.

• 4. It is likely that several orders of magnitude more mass-energy are present in the remnant disk than is observed in EM radiation from a SGRB. Modeling the emission from the disk (and possibly a HMNS) remains extremely challenging. Neutrino leakage schemes have been applied in both approximate relativistic calculations [252, 246] and full GR [265], and a more complex flux-limited diffusion scheme has been applied to the former as a post-processing step [82], but there are no calculations that follow in detail the neutrinos as they flow outward, annihilate, and produce observable EM emission. At present, nuclear reactions are typically not followed in detail; rather, the electron fraction of the nuclear material, Y_e, is evolved, and used to calculate neutrino emission and absorption rates.

• 5. Magnetic fields, on the other hand, are starting to be much better understood. B-fields do seem to grow quite large through winding effects, even during the limited amount of physical time that can currently be modeled numerically [6, 172, 332, 241], with some calculations indicating exponential growth rates. The resulting geometries seem likely to produce the disk/jet structure observed throughout astrophysics when magnetized objects accrete material, which span scales from stellar BHs or pre-main sequence stars all the way up to active galactic nuclei [241].

• 6. While recent numerical simulations have strengthened the case for NS-NS mergers as SGRB progenitors, full GR calculations have not generated much support for the same events yielding significant amounts of r-process elements. Noting the standard caveat that low-density ejecta are difficult to model, and that nuclear reactions are rarely treated self-consistently, there is still tension between CF calculations producing ejecta with the proper temperatures and masses to reproduce the observed cosmic r-process abundances (see, e.g., [123]), and full GR calculations that produce almost no measurable unbound material whatsoever.

Figure 18

Summary of potential outcomes from NS-NS mergers. Here, M_thr is the threshold mass (given the EOS) for collapse of a HMNS to a BH, and Q_m is the binary mass ratio. ‘Small’, ‘massive’, and ‘heavy’ disks imply total disk masses M_disk ≪ 0.01 M_⊙, 0.01 M_⊙ ≲ M_disk ≲ 0.03 M_⊙, and M_disk ≳ 0.05 M_⊙, respectively. ‘B-field’ and ‘J-transport’ indicate potential mechanisms for the HMNS to eventually lose its differential rotation support and collapse: magnetic damping and angular momentum transport outward into the disk. Spheroids are likely formed only for the APR and other stiff EOS models that can support remnants with relatively low rotational kinetic energies against collapse. Image reproduced by permission from [282], copyright by APS.

While NS-NS merger calculations have seen tremendous progress in the past decade, the future remains extremely exciting. Between the addition of more accurate and realistic physical treatments, the exploration of the full phase space of models, and the linking of numerical relativity to astrophysical observations and GW detection, there remain many unsolved problems that will be attacked over the course of the next decade and beyond.

Appendices

Field evolution equations

In the BSSN evolution system, we define the following variables in terms of the standard ADM 4-metric g_ij, 3-metric γ_ij, and extrinsic curvature K_ij:

$$\phi \equiv {1 \over {12}}\log [\det {\gamma _{ij}}] = \log \psi,$$

(35)

$${\tilde \gamma _{ij}} \equiv {e^{- 4\phi}}{\gamma _{ij}},$$

(36)

$$K \equiv {g^{ij}}{K_{ij}},$$

(37)

$${\tilde A_{ij}} \equiv {e^{- 4\phi}}\left[ {{K_{ij}} - {1 \over 3}{g_{ij}}K,} \right]$$

(38)

$${\tilde \Gamma ^i} \equiv {\tilde \gamma ^{jk}}\tilde \Gamma _{jk}^i.$$

(39)

The evolution system consists of 15 equations for the various field terms,

$${\partial _0}K = - {\gamma ^{ij}}{\tilde D_i}{\tilde D_j}\alpha + \alpha \left({{{\tilde A}^{ij}}{{\tilde A}_{ij}} + {1 \over 3}{K^2}} \right) + 4\pi (E + {\gamma ^{ij}}{S_{ij}}),$$

(40)

$${\partial _0}\phi = - {1 \over 6}(\alpha K - {\partial _k}{\beta ^k}),$$

(41)

$${{\partial _0}{{\tilde \gamma}_{ij}} = - 2\alpha {{\tilde A}_{ij}} + 2{{\tilde \gamma}_{k(i}}{\partial _{j)}}{\beta ^k} - {2 \over 3}{{\tilde \gamma}_{ij}}{\partial _k}{\beta ^k},}$$

(42)

$$\begin{array}{*{20}c}{{\partial _0}{{\tilde A}_{ij}} = {e^{- 4\phi}}{{\left[ {\alpha {{\tilde R}_{ij}} + \alpha R_{ij}^\phi - {{\tilde D}_i}{{\tilde D}_j}\alpha} \right]}^{TF}}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {+ \alpha K{{\tilde A}_{ij}} - 2\alpha {{\tilde A}_{ik}}\tilde A_j^k + 2{{\tilde A}_{k(i}}{\partial _{j)}}{\beta ^k} - {2 \over 3}{{\tilde A}_{ij}}{\partial _k}{\beta ^k} - 8\pi \alpha {e^{- 4\phi}}S_{ij}^{TF},}\end{array}$$

(43)

$$\begin{array}{*{20}c}{{\partial _0}{{\tilde \Gamma}^i} = - 2{{\tilde A}^{ij}}{\partial _j}\alpha + 2\alpha \left[ {\tilde \Gamma _{kl}^i{{\tilde A}^{kl}} + 6{{\tilde A}^{ij}}{\partial _j}\phi - {2 \over 3}{{\tilde \gamma}^{ij}}{K_{,j}}} \right]\quad \quad \quad \quad \quad \quad \quad} \\ {- {{\tilde \Gamma}^j}{\partial _j}{\beta ^i} + {2 \over 3}{{\tilde \Gamma}^i}{\partial _j}{\beta ^j} + {1 \over 3}{{\tilde \gamma}^{ik}}\beta _{,jk}^j + {{\tilde \gamma}^{jk}}\beta _{,jk}^i - 16\pi \alpha {{\tilde \gamma}^{ik}}{S_k},}\end{array}$$

(44)

where the matter source terms contain various projections of the stress-energy tensor, defined through the relations

$$E \equiv {n^\mu}{n^\nu}{T_{\mu \nu}} = {1 \over {{\alpha ^2}}}({T_{00}} - 2{\beta ^i}{T_{0i}} + {\beta ^i}{\beta ^j}{T_{ij}}),$$

(45)

$${S_i} \equiv - {n^\mu}\gamma _i^\nu {T_{\mu \nu}} =- {1 \over \alpha}({T_{0i}} - {\beta ^j}{T_{ij}}),$$

(46)

$${S_{ij}} \equiv \gamma _i^\mu \gamma _j^\nu {T_{\mu \nu}}.$$

(47)

We have introduced the notation ∂₀ = ∂_t − β^j∂_j. All quantities with a tilde involve the conformal 3-metric ${{\tilde \gamma}_{ij}}$, which is used to raise and lower indices. In particular, ${{\tilde D}_i}$ and $\tilde \Gamma _{ij}^k$ refer to the covariant derivative and the Christoffel symbols with respect to ${{\tilde \gamma}_{ij}}$. Parentheses indicate symmetrization of indices, and the expression […]^TF denotes the trace-free part of the expression inside the brackets. In the BSSN approach, the Ricci tensor is typically split into two pieces, whose respective contributions are given by

$${\tilde R_{ij}} = - {1 \over 2}{\tilde \gamma ^{kl}}{\partial _k}{\partial _l}{\tilde \gamma _{ij}} + {\tilde \gamma _{k(i}}{\partial _{j)}}{\tilde \Gamma ^k} - {\partial _k}{\tilde \gamma _{l(i}}{\partial _{j)}}{\tilde \gamma ^{kl}} + {1 \over 2}{\tilde \Gamma ^l}{\tilde \gamma _{ij,l}} - \tilde \Gamma _{ik}^l\tilde \Gamma _{jl}^k,$$

(48)

$$R_{ij}^\phi = - 2{\tilde D_i}{\tilde D_j}\phi - 2{\tilde \gamma _{ij}}{\tilde D^k}{\tilde D_k}\phi + 4{\tilde D_i}\phi {\tilde D_j}\phi - 4{\tilde \gamma _{ij}}{\tilde D^k}\phi {\tilde D_k}\phi.$$

(49)

These equations must be supplemented with gauge conditions that determine the evolution of the lapse function α and shift vector βⁱ. Noting that some groups introduce slight variants of these, the moving puncture gauge conditions that have become popular for all GR merger calculations involving BHs and NSs typically take the form

$${\partial _0}\alpha = - 2\alpha K,$$

(50)

$${\partial _t}{\beta ^i} = {3 \over 4}{B^i},$$

(51)

$${\partial _t}{B^i} = {\partial _t}{\tilde \Gamma ^i} - \eta {B^i},$$

(52)

where Bⁱ is an intermediate quantity used to convert the second-order “Gamma-driver” shift condition into a pair of first-order equations, and η is a user-prescribed term used to control dissipation in the simulation. Note that it is possible to replace the three instances of ∂_t in the shift evolution equations 51 and 52 by the shift-advected time derivative (∂_t − β^j∂_j) without changing the stability or hyperbolicity properties of the evolution scheme; in both cases moving punctures translate smoothly across a grid over long time periods [322] and both systems are strongly hyperbolic so long as the shift vector does not grow too large within the simulation domain [129].

The generalized harmonic formulation involves recasting the Einstein field equations, Eq. 9 in the form

$${R_{\mu \nu}} = 8\pi \left({{T_{\mu \nu}} - {1 \over 2}{g_{\mu \nu}}T} \right)$$

(53)

and, after some tensor algebra, rewriting the Ricci tensor in the form

$$- 2{R_{\mu \nu}} = {g^{\gamma \delta}}{g_{\alpha \beta, \gamma \delta}} + {g^{\gamma \delta}}_{,\beta}{g_{\alpha \delta, \gamma}} + g_{,\alpha}^{\gamma \delta}{g_{\beta \delta, \gamma}} + 2{H_{(\alpha, \beta)}} - 2{H_\delta}\Gamma _{\alpha \beta}^\delta + 2\Gamma _{\delta \beta}^\gamma \Gamma _{\gamma \alpha}^\delta.$$

(54)

The Christoffel coefficients are calculated from the full 4-metric,

$$\Gamma _{\beta \gamma}^\alpha \equiv {1 \over 2}{g^{\alpha \delta}}[{g_{\beta \delta, \gamma}} + {g_{\gamma \delta, \beta}} - {g_{\beta \gamma, \delta}}],$$

(55)

and the gauge source terms H^μ are defined in Eq. 25.

Given well-posed initial data for the metric and its first time derivative (since the system is second-order in time according to Eq. 54), the evolution of the system may be treated by a first-order reduction that specifies the evolution of the four functions H^μ along with the spacetime metric g_μν, its projected time derivatives Π_μν = −n>^α∂_αg_μν, and its spatial derivatives

$${\Phi _{i\mu \nu}} = {\partial _i}{g_{\mu \nu}},$$

(56)

subject to a constraint specifying that the derivative terms Φ_iμν remain consistent with the metric g_μν in time. In practice, one typically introduces a constraint for the source functions, defining

$${C^\mu} = {H^\mu} - \square{x^\mu}$$

(57)

and then modifies the evolution equation by appending a constraint damping term to the RHS of the stress energy-term (following [130, 231]

$$- 8\pi (2{T_{\mu \nu}} - {g_{\mu \nu}}T)\quad \Rightarrow \quad - 8\pi (2{T_{\mu \nu}} - {g_{\mu \nu}}T) - \kappa ({n_\mu}{C_\nu} + {n_\nu}{C_\mu} - {g_{\mu \nu}}{n^\delta}{C_\delta}),$$

(58)

where n_μ is the unit normal vector to the hypersurface (see Eq. 15). The gauge conditions used in the first successful simulations of merging BH binaries [231] consisted of the set

$$\square{H_t} = - {\xi _1}{{\alpha - 1} \over {{\alpha ^\eta}}} + {\xi _2}{H_{t,\nu}}{n^\nu};\quad {H_i} = 0,$$

(59)

with ξ₁, ξ₂, η constant parameters used to tune the evolution. The first one drives the coordinates toward the ADM form and the latter provides dissipation. The binary NS-NS merger work of [6] chose harmonic coordinates with H^μ = 0.

GR Hydrodynamical and MHD equations

In what follows, we will adopt the stress energy tensor of an ideal relativistic fluid,

$${T^{\mu \nu}} = \rho h{u^\mu}{u^\nu} + {g^{\mu \nu}}P,$$

(60)

where ρ, P, and u^μ are the rest mass density, pressure, and fluid 4-velocity, respectively, and

$$h = 1 + \epsilon + P/\rho$$

(61)

is the relativistic specific enthalpy, with ϵ the specific internal energy of the fluid.

The equations of ideal GR hydrodynamics [186] may be derived from the local GR conservation laws of mass and energy-momentum:

$${\nabla _\mu}{J^\mu} = 0,\quad {\nabla _\mu}{T^{\mu \nu}} = 0,$$

(62)

where ∇_μ denotes the covariant derivative with respect to the 4-metric, and J^μ = ρu^μ is the mass current.

The 3-velocity υⁱ can be calculated in the form

$${v^i} = {{{u^i}} \over W} + {{{\beta ^i}} \over \alpha},$$

(63)

where

$$W \equiv \alpha {u^0} = {(1 - {v^i}{v_i})^{- 1/2}}$$

(64)

is the Lorentz factor. The contravariant 4-velocity is then given by:

$${u^0} = {W \over \alpha},\quad {u^i} = W\left({{v^i} - {{{\beta ^i}} \over \alpha}} \right),$$

(65)

and the covariant 4-velocity is:

$${u_0} = W({v^i}{\beta _i} - \alpha),\quad {u_i} = W{v_i}.$$

(66)

To cast the equations of GR hydrodynamics as a first-order hyperbolic flux-conservative system for the conserved variables D, Sⁱ, and τ, defined in terms of the primitive variables ρ, ϵ, υⁱ, we define

$$D = \sqrt \gamma \rho W,$$

(67)

$${S^i} = \sqrt \gamma \rho h{W^2}{v^i},$$

(68)

$$\tau = \sqrt \gamma (\rho h{W^2} - P) - D,$$

(69)

where γ is the determinant of γ_ij. The evolution system then becomes

$${{\partial {\bf{U}}} \over {\partial t}} + {{\partial {{\bf{F}}^i}} \over {\partial {x^i}}} = {\bf{S}},$$

(70)

with

$$\begin{array}{*{20}c} {{\bf{U}} = [D,{S_j},\tau ],}\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\\ {{{\bf{F}}^i} = \alpha [D{{\tilde v}^i},{S_j}{{\tilde v}^i} + \delta _j^iP,\tau {{\tilde v}^i} + P{v^i}],}\quad\quad\quad\\ {{\bf{S}} = \alpha \left[ {0,{T^{\mu \nu}}\left({{{\partial {g_{\nu j}}} \over {\partial {x^\mu}}} - \Gamma _{\mu \nu}^\lambda {g_{\lambda j}}} \right),} \right.}\quad\quad\quad\quad\\ {\left. {\alpha \left({{T^{\mu 0}}{{\partial \ln \alpha} \over {\partial {x^\mu}}} - {T^{\mu \nu}}\Gamma _{\mu \nu}^0} \right)} \right].}\\ \end{array}$$

(71)

Here, ${{\tilde v}^i} = {v^i} - {\beta ^i}/\alpha$ and $\Gamma _{\mu \nu}^\lambda$ are the 4-Christoffel symbols.

Magnetic fields may be included in the formalism, in the ideal MHD limit under which we assume infinite conductivity, by adding three new evolution equations and modifying those above to include magnetic stress-energy contributions of the form

$$T_{\mu \nu}^{EM} = {b^2}\left({{u^\mu}{u^\nu} + {1 \over 2}{g^\mu}\nu} \right) - {b^\mu}{b^\nu},$$

(72)

where the magnetic field seen by a comoving observer, b^μ is defined in terms of the dual Faraday tensor *F^νμ by the condition

$${b^\mu}{= ^{\ast}}{F^{\nu \mu}}{u_\nu},$$

(73)

where b² = b^μb_μ represents twice the magnetic pressure. With magnetic terms included, we may rewrite the stress-energy tensor in a familiar form by introducing magnetically modified pressure and enthalpy contributions:

$${T^{\mu \nu}} = \rho {h^{\ast}}{u^\mu}{u^\nu} + {P^{\ast}}{g^{\mu \nu}};\;\;{h^{\ast}} \equiv 1 + \epsilon + {{P + {b^2}} \over \rho},\;\;{P^{\ast}} = P + {{{b^2}} \over 2}$$

(74)

and redefine the conserved momentum and energy variables Sⁱ and τ accordingly:

$${S^i} = \sqrt \gamma \rho {h^{\ast}}{W^2}{v^i} - \alpha {b^0}{b^i},$$

(75)

$$\tau = \sqrt \gamma (\rho {h^{\ast}}{W^2} - {P^{\ast}} - {(\alpha {b^0})^2}) - D{.}$$

(76)

Defining the (primitive) magnetic field 3-vector as

$${B^i} = {- ^{\ast}}{F^{\mu i}}{n_\mu} = \alpha {\;^{\ast}}{F^{0i}}$$

(77)

and the conserved variable ${{\mathcal B}^i} = \sqrt \gamma {B^i}$, which are related to the comoving magnetic field 4-vector b^μ through the relations

$${b^0} = {{W{B^k}{v_k}} \over \alpha},$$

(78)

$${b^i} = {{{B^i}} \over W} + W({B^k}{v_k})\left({{v^i} - {{{\beta ^i}} \over \alpha}} \right),$$

(79)

$${b_i} = {{{B_i}} \over W} + W({B^k}{v_k}){v_i},$$

(80)

$${b^2} = {{{B^i}{B_i}} \over {{W^2}}} + {({B^i}{v_i})^2},$$

(81)

we may rewrite the conservative evolution scheme in the form

$$\begin{array}{*{20}c} {{\bf{U}} = [D,{S_j},\tau, {{\mathcal B}^k}],}\quad\quad\quad\quad\quad\quad\\ {{{\bf{F}}^i} = \alpha \times \left[ {\begin{array}{*{20}c} {D{{\tilde v}^i},}\\ {{S_j}{{\tilde v}^i} + \delta _j^i{P^{\ast}} - {b_j}{{\mathcal B}^i}/W,}\\ {\tau {{\tilde v}^i} + {P^{\ast}}{v^i} - \alpha {b^0}{{B}^i}/W,}\\ {\quad \quad {{\mathcal B}^k}{{\tilde v}^i} - {{\mathcal B}^i}{{\tilde v}^k}}\\ \end{array}} \right]}\\ \end{array},$$

(82)

$${\bf{S}} = \alpha \times \left[ \begin{array}{*{20}c}{0,} \\{\quad {T^{\mu \nu}}\left({{{\partial {g_{\nu j}}} \over {\partial {x^\mu}}} - \Gamma _{\mu \nu}^\lambda {g_{\lambda j}}} \right),\quad} \\ {\quad \alpha ({T^{\mu 0}}{{\partial \ln \alpha} \over {\partial {x^\mu}}} - {T^{\mu \nu}}\Gamma _{\mu \nu}^0),} \\ {\vec 0} \\ \end{array} \right],$$

(83)

where the magnetic field evolution equation is just the relativistic version of the induction equation. An external mechanism to enforce the divergence-free nature of the magnetic field, ${\partial _i}{{\mathcal B}^i}$ must also be applied.