GRAVITATIONAL MEMORY IN BINARY BLACK HOLE MERGERS

Accelerating massive bodies generate gravitational waves (GWs), a potentially rich source of astrophysical information which a number of large experiments have been designed to mine over the next decade. The principle sources for the interferometric detectors (LIGO, Virgo, LISA) are orbiting bodies, for which the dominant signal is oscillatory. However, it has long been known that GWs will also contain non-oscillatory features, resulting from a net change in time derivatives of the multiple moments of the system (Zel'Dovich & Polnarev 1974; Braginskii & Thorne 1987). A nonlinear contribution also results from the interaction of the waves with themselves (Payne 1983; Christodoulou 1991; Blanchet & Damour 1992). Evaluating these modes involves an integral in time, requiring knowledge of the entire past history of the spacetime. As such, they have been given the name "memory" or "hereditary" components of the GW.

The nonlinear memory is sourced by emitted gravitational energy. For a relativistic binary inspiral, the typical profile is of a slow growth over time which sees a rapid increase during the late-inspiral and merger, and reaches a constant value as the merger remnant rings down and ceases to emit. The non-oscillatory nature of the memory suggests that it will not be a dominant feature in interferometric detectors, for which drifts in the background metric are factored out. However, the step induced during merger may be observable in pulsar timing measurements, as has been noted in a recent set of papers (Pshirkov et al. 2009; van Haasteren & Levin 2009; Seto 2009).

Post-Newtonian (PN) calculations provide an accurate estimate of the memory to within a few orbits of the merger (Kennefick 1994), however they fail to model the important merger phase where the effect is largest. In a series of papers, Favata (2009a, 2009b, 2009c, 2010) has included an estimate of the merger contribution for binary black holes (BHs) by using an effective-one-body model which was tuned using the results of numerical simulations. Recent progress in numerical models of BH spacetimes provides the opportunity to measure the effect directly, though as outlined in Favata (2009c), the problem is challenging due to the low amplitude of the memory and systematic error in standard techniques of wave measurement in numerical relativity.

In this Letter, we apply a newly developed evolution code (Pollney et al. 2009a, 2009b) and wave extraction techniques (Reisswig et al. 2009a, 2010) to carry out fully relativistic binary BH simulations to directly measure the nonlinear GW memory through inspiral, merger, and ringdown in a gauge-invariant and mathematically unambiguous notion. We concern ourselves with equal-mass, spinning binaries. For these models, the memory is contained in the ₋₂Y_ℓ0 spin-weighted spherical harmonic components of the GW strain, h⁺. These modes, previously studied in the context of head-on collisions (Anninos et al. 1995), exhibit two dominant effects. The first is a non-oscillatory term which rises as the inspiral progresses and is associated with the memory. During the merger, an oscillatory signal is superposed onto these modes, induced by the ringdown of the BH remnant. After the ringdown has subsided, the modes are offset from their original value. We find that for the dominant (ℓ, m) = (2, 0) memory mode, this offset is largest in the case of the merger of aligned, maximally spinning BHs, as might be expected given that these are the strongest gravitational emitters (Reisswig et al. 2009b). Interestingly, we observe the strongest oscillatory ringdown signal in the case of lowest total spin (maximally spinning BHs anti-aligned with the orbital angular momentum). As such, and as opposed to the observability of the dominant (ℓ, m) = (2, 2) mode, these models produce a more visible (ℓ, m) = (2, 0) mode to interferometric detectors. We conclude by estimating the strain offset which would be visible to a pulsar timing array.

We integrate the BSSNOK formulation of the vacuum Einstein equations (see, e.g., Alcubierre 2008) numerically, using the Llama evolution code (Pollney et al. 2009a, 2009b) within the Cactus framework (Goodale et al. 2003). The evolution equations are discretized via finite differences in space on a grid composed of multiple patches with locally adapted coordinates in the wave zone. Adaptive mesh-refinement is implemented via the Carpet grid driver (Schnetter et al. 2004) on the central patch surrounding the BHs. Importantly, for the measurement of small non-oscillatory features such as the memory, the artificial grid outer boundary is causally disconnected from the wave measurements.

We determine initial data for binary BHs by the conformal puncture data method (Brandt & Brügmann 1997) and evolve them using the "moving puncture" approach (Baker et al. 2006; Campanelli et al. 2006) whereby the choice of gauge (Alcubierre et al. 2003) prevents the spacetime slicing from encountering the curvature singularity (Hannam et al. 2007).

Crucial for determining the memory, we compute the GW signal of the source using the method of Cauchy-characteristic extraction (CCE), which computes unambiguous and coordinate invariant signals at future null infinity, $\mathcal {J}^+$, corresponding to a detector far removed from the source (Bishop et al. 1999; Babiuc et al. 2005; Winicour 2005; Reisswig et al. 2009a, 2010). By this technique, metric data is collected on a world tube at finite radius, and a formulation of the full Einstein equations along null hypersurfaces is used to transport the signal to $\mathcal {J}^+$.

We measure two quantities that encode the gravitational signal at $\mathcal {J}^+$. The first is the Newman–Penrose scalar ψ₄, the Weyl curvature component with slowest falloff in asymptotically flat spacetimes (Newman & Penrose 1962). Alternatively, we measure the Bondi "news" (Bondi et al. 1962; Sachs 1962), defined by

$$\begin{equation} \mathcal {N}= -\Delta \bar{\sigma }, \end{equation} \tag{ 1 }$$

where Δ = l^a∇_a is a derivative operator defined on an outgoing null geodesic, l^a, and σ is the Newman–Penrose shear scalar (Newman & Penrose 1962). These scalars, ψ₄ and $\mathcal {N}$, can be evaluated directly from the local curvature on the sphere at $\mathcal {J}^+$and in the naturally defined gauge of an asymptotic observer (Bondi et al. 1962; Sachs 1962). To calculate the GW strain h which is observed by a detector, these quantities need to be numerically integrated in time:

$$\begin{equation} h = h_+ - ih_\times = \int ^t_{-\infty } dt^{\prime } \mathcal {N}\, = \int ^{t}_{-\infty } dt^\prime \int ^{t^\prime }_{-\infty } dt^{\prime \prime } \psi _4\,. \end{equation} \tag{ 2 }$$

We expand h in terms of a basis of spin-2 spherical harmonics ₋₂Y_ℓm. The dominant non-oscillatory hereditary component is contained in the (ℓ, m) = (2, 0) mode,

$$\begin{equation} r\,h_{20} \equiv \int _\Omega d\Omega ({}_{-2}\bar{Y}_{20}) h, \end{equation} \tag{ 3 }$$

where Ω is the sphere at $\mathcal {J}^+$and r is the distance to the source.

Since the simulations are necessarily of finite length, the time integration of $\mathcal {N}$, Equation (2), leaves a constant to be determined, corresponding to the initial value of the signal at the start time of the simulation. For purely oscillatory signals (e.g., h₂₂) we can fit this by adjusting the post ringdown value to zero or by some averaging procedure over multiple inspiral wavelengths. This is not possible in the case of the monotonically growing memory. Instead, we approximate an initial value of the h₂₀ mode by matching to a PN calculation, using the 3PN expansion for h^mem derived in Favata (2009c). For spinning binaries we incorporate 2.5PN spin contributions to the frequency evolution (Blanchet et al. 2006). The initial offset of the numerical waveforms is performed by shifting its amplitude so that it fits against the PN models, over an interval from t = −300 M to t = −200 M. For instance, the PN solution for h₂₀ in the non-spinning case is plotted as a dotted line in the lower panel of Figure 1. In Table 1, we report both the overall amplitude offset of the ringdown signal in the h₂₀ mode matched to the PN inspiral as well as the offset from the point t = −300 M, which is computed entirely from numerical data.

Zoom In Zoom Out Reset image size

In Figure 1, we show the (ℓ, m) = (2, 0) mode during late-inspiral and merger for the case of an equal-mass, non-spinning binary using data from the model presented in Pollney et al. (2009b) and evaluated at $\mathcal {J}^+$ in Reisswig et al. (2009a). The upper panel shows the development of the (ℓ, m) = (2, 0) component of the news, $\mathcal {N}_{20}$, whereas the lower shows its integral, h₂₀. The time coordinate places t = 0 M at the position of the amplitude peak of the dominant (2, 2) mode (not shown), corresponding roughly to the time of merger. Two main features are present. The first is a slow monotonic growth of the $\mathcal {N}_{20}$ mode during the inspiral and plunge phase. This corresponds to the non-oscillatory memory contribution. At t ≃ −20 M, there is a local maximum in $\mathcal {N}_{20}$, and shortly afterward a prominent oscillatory signal sets in. This signal represents the binary merger and subsequent BH ringdown. It is not strictly an aspect of the GW "memory," as it is not dependent on the entire past history of the spacetime, however it is superposed on the growing memory signal.

The lower plot shows the strain, determined via Equation (2). The superposition of the memory and ringdown signals results in a notable kink in the waveform near t = 0 M and exponentially decaying oscillation before reaching the steady-state amplitude of

$$\begin{equation} h^{\rm tot}_{20}=0.097\pm 2\times 10^{-3}. \end{equation} \tag{ 4 }$$

We measure approximately fourth-order convergence in $\mathcal {N}_{20}$ when comparing three different grid resolutions for the same model, and we estimate an error on the order of 1 × 10⁻⁴ in $\mathcal {N}_{20}$. The much larger error of 2 × 10⁻³ in h^tot₂₀ is an upper limit of the estimate on the total error composed of integration error (2) and error in the PN fit for all of the binary models considered in this Letter.

Perturbative results provide a prediction for the observed quasi-normal modes (QNM) of $\mathcal {N}_{20}$ (Berti et al. 2009). We perform a least-squares fit of the numerical results to a function of the form

$$\begin{equation} f(t) = A \exp (-t / \tau) \cos (\omega _{20} t + \phi), \end{equation} \tag{ 5 }$$

over the ringdown portion, t ∈ [30, 80], with fitting parameters, A, and w₂₀, τ, ϕ, and arrive at a frequency

$$\begin{equation} M_{\rm f}\omega _{20} = (0.3940 \pm 2)\times 10^{-4}\,, \end{equation} \tag{ 6 }$$

using the measured mass of the remnant of M_f = (0.951764 ± 20) × 10⁻⁶ (Pollney et al. 2009b). Berti et al. (2009) provide a tabulated value of Mω^lit.₂₀ = 0.393245 for the (2, 0) mode of a BH with the measured dimensionless spin a_f = (0.686923 ± 1) × 10⁻⁵. Thus, we find a difference of

$$\begin{equation} \big|M_{\rm f}\omega _{20}-M\omega _{20}^{\rm lit.}\big|\approx 8\times 10^{-4}\,, \end{equation} \tag{ 7 }$$

representing an error of less than 0.2%. We list the final spin and the associated QNM frequencies of the binary systems for all simulations performed in Table 1.

Finally, we have also examined m = 0 modes for ℓ = 4 and ℓ = 6. For this model, the overall amplitudes of the final memory are h^tot₄₀ = 1.7 × 10⁻³ and h^tot₆₀ ≈ 5.0 × 10⁻⁴. However, these modes are difficult to distinguish from the numerical error, and we conclude that the higher order modes contribute less than 2% to the total memory.

We have performed simulations using the code infrastructure described in Section 2 in the two-dimensional parameter space of non-precessing equal-mass binaries with spin that has also been studied in Reisswig et al. (2009b) and Rezzolla et al. (2008). We have focused on initial data configurations along the two main axes in the space of individual dimensionless BH spins (a₁, a₂) aligned with the orbital angular momentum. The first set of models has equal but opposite spins, a₁ = −a₂; the second has identical spins, a₁ = a₂. The model details are listed in Table 1.

Figure 2 plots the evolution of the h₂₀ modes for the two sequences. The upper panel shows models for which the spins of the individual BHs are anti-aligned, a₁ = −a₂, so that the spacetimes have the same total angular momentum. The nearly identical lines are consistent with previous observations that the (ℓ, m) = (2, 2) modes also do not vary appreciably along this direction of parameter space (Vaishnav et al. 2007; Reisswig et al. 2009b), as well as the radiated energies (Reisswig et al. 2009b), at least not to within the error bars as stated in the latter reference.

Zoom In Zoom Out Reset image size

The lower plot shows the h₂₀ modes for configurations with equal and aligned spins, a₁ = a₂. The amplitude of the memory offset increases with higher total spin. This is expected as the memory is sourced by the emitted GWs which are known to be more energetic in the higher spin models (Reisswig et al. 2009b; Lousto et al. 2010). The total offsets, h^tot₂₀, are listed in Table 1.

We find that the amplitude of the ringdown is largest for models with lowest total spin, a₁ = a₂ = −0.8. In these cases, spin–orbit coupling increases the final angle of impact so that the merger is more nearly head-on. The ringdown is suppressed in the high-spin a₁ = a₂ = +0.8 model. Fits to the QNM frequencies are listed in Table 1 and show good agreement with perturbative results. However, for simulations with a₁ = a₂ > 0 the weakness of the ringdown contributes to larger variation in the estimates, particularly if the fitting interval is varied, though still within 5% of the analytic values.

The memory is well estimated by a quartic polynomial in the total spin a = (a₁ + a₂)/2,

$$\begin{eqnarray} h^{\rm tot}_{20}(a) r &=& 0.0969 + 0.0562\,a + 0.0340\,a^2 + 0.0296\,a^3 \nonumber \\ && + 0.0206\,a^4\,, \end{eqnarray} \tag{ 8 }$$

determined by a least-squares fit to the measured data. Figure 3 plots h^tot₂₀ as a function of spin for the aligned models, a₁ = a₂, as well as the leading order multipole moment estimate resulting from the radiated energy (see Favata 2009a; Equation (5))

$$\begin{equation} h^{\rm mem}_+ \simeq \frac{\eta M h^{\rm mem}}{384\pi R} \sin ^2\theta (17 + \cos ^2\theta), \end{equation} \tag{ 9 }$$

with

$$\begin{equation} h^{\rm mem} \simeq \frac{16\pi }{\eta }\left(\frac{\Delta E_{\rm rad}}{M}\right). \end{equation} \tag{ 10 }$$

The measured radiated energy, ΔE_rad, can also be estimated using the quadratic fit developed in Reisswig et al. (2009b), verified here with the plotted data points corresponding to the new and more accurate simulations (see Table 1). The radiated energy provides a reasonably good estimate of the memory offset, though with a slight underprediction. The measured memory values also appear to have a somewhat stronger dependence on the total spin. Extrapolating to the extremal a₁ = a₂ = 1 case, we estimate a maximum h^tot₂₀ ≃ 0.24 in the case of aligned spins. Fits to ΔE_rad for generic spin configurations (Lousto et al. 2010) suggest that this will indeed be the maximum value for arbitrary spins. Additional simulation in the high-spin regime would be required to establish the accuracy of the extrapolation.

Zoom In Zoom Out Reset image size

The m = 0 modes are small (on the order of 10%) compared to the dominant ℓ = m components of the GW signal. However, they exhibit both non-oscillatory growth and ringdown features, which present some interesting aspects for detection and identification of these modes. We discuss the prospects for observation in both interferometers, and pulsar timing arrays, in the following sections.

5.1. Interferometers

The signal-to-noise ratio (S/N) within a given detector, ρ, in GW searches based on matched filtering is given by Flanagan & Hughes (1998):

$$\begin{equation} \rho ^2 \equiv \left(\frac{{\rm S}}{{\rm N}}\right)^2_{\rm matched} = 4 \int _0^\infty \frac{|\tilde{h}(f)|^2}{S_h(f)} df \,, \end{equation} \tag{ 11 }$$

where S_h(f) is the sensitivity curve of the instrument. Restricting the integral to the (ℓ, m) = (2, 0) mode, we can determine its contribution to the overall S/N, testing against the proposed advanced-LIGO (Shoemaker 2010) and LISA noise curves.³

Table 1 lists the maximum S/N attained over a range of masses for a given model. For advanced-LIGO, the optimal total mass is M_BH ≃ 290 M_☉ (independent of the model and invariant under the distance), i.e., a 145 M_☉ + 145 M_☉ binary, for which we list the S/N at a reference distance of 300 Mpc (the S/N scales linearly with the distance). The LISA results refer to an optimal (redshifted) total mass of M_BH = 5.35 × 10⁶ M_☉ at a luminosity distance d = 15.8 Gpc corresponding to a redshift z = 2, where models suggest at least one expected merger event per year (Sesana et al. 2005). By assuming a minimum S/N of ρ = 8 for detection, it will be possible to observe the (ℓ, m) = (2, 0) mode in LISA out to the given distance. They are unlikely to be visible with significant event rates within the planned advanced ground-based detectors, though future generation experiments may have prospects.

For both detectors, and as opposed to the dominant (ℓ, m) = (2, 2) mode (Reisswig et al. 2009b), we note that the highest S/N occurs for the cases with spins anti-aligned to the orbital angular momentum which is due to their stronger ringdown. In the aligned cases, the weaker ringdown causes a much flatter high-frequency falloff (Figure 4), and this effect overrides the non-oscillatory step within these detectors.

Zoom In Zoom Out Reset image size

An important aspect of the matched-filtering algorithm is the construction of templates against which signals are compared. We examine the influence of neglecting the (2, 0) contribution by computing the mismatch between a signal, h(θ = π/3, ϕ = π/2), containing all of the (2, m) modes, with one in which the (2, 0) mode is left out. For advanced-LIGO, the mismatch for the model a₁ = a₂ = +0.8 is largest in the mass range 250 M_☉–400 M_☉, but negligible, with a magnitude of 10⁻¹¹. For the a₁ = a₂ = −0.8 model, the mismatch grows to be on the order of 10⁻⁵. Thus, although the h₂₀ mode induces a notable offset in the wave signal when viewed in the time domain, it makes little difference to the detection algorithms.

5.2. Pulsar Timing Arrays

An alternative proposal for measuring GWs may be more sensitive to the non-oscillatory step-function nature of the memory. Precise, distant clocks will experience a residual in time-of-arrival measurements, which can be associated with GWs (see, e.g., Detweiler 1979). Several such projects are in development (Hobbs et al. 2010; Manchester 2008; Jenet et al. 2009; Lazio 2009). While their principal aim is to examine the stochastic background, recently a number of papers have pointed out the potential for observing non-oscillatory step-function burst signals, such as the memory (Pshirkov et al. 2009; van Haasteren & Levin 2009; Seto 2009).

Primary candidates for detection are supermassive BH (SMBH) binary mergers. While there is a great deal of ongoing work understanding the nature and evolution of such systems, if there is a tendency for the bodies to align, then the results obtained here are directly applicable. For instance, Dotti et al. (2009) indicate that in gas-rich mergers, the spins are aligned on a short timescale, with spin magnitudes in the range 0.6–0.8. Assuming a m₁ = m₂ = 10⁸ M_☉ binary source at a distance of 1 Gpc which has aligned spins with a magnitude of a₁ = a₂ = 0.7, the results of Section 4 lead to an angle-averaged memory offset of

$$\begin{equation} \langle h\rangle _\theta = 1.6 \times 10^{-16} \left(\frac{m}{10^8 \,M_\odot } \right) \left(\frac{1\,{\rm Gpc}}{R} \right). \end{equation} \tag{ 12 }$$

(The corresponding constants for zero and maximally spinning bodies are 0.9 × 10⁻¹⁶ and 2.0 × 10⁻¹⁶, respectively.) The results of Pshirkov et al. (2009) and van Haasteren & Levin (2009) indicate that a burst will be observable over a 10 year observation period at amplitude h ≃ 2 × 10⁻¹⁵, suggesting that the merger of a pair of 6 × 10⁸ with a₁ = a₂ = 0.7 would be observable at 1 Gpc.

If we assume rather conservative estimates of SMBH merger event rates, e.g., Enoki et al. (2004) and Sesana et al. (2005) (Berti 2006 for a review), we would have a merger rate of 0.1 yr⁻¹ at redshifts z ≲ 1 for ∼10⁸  M_☉ binary sources. Hence, we can expect to detect approximately one such event over a period of 10 years.

The authors thank Marc Favata, Ian Hinder, Sascha Husa, and Christian D. Ott for helpful input. This work is supported by the Bundesministerium für Bildung und Forschung and the National Science Foundation under grant numbers AST-0855535 and OCI-0905046. D.P. has been supported by grants CSD2007-00042, FPA-2007-60220, and FPA2010-16495 of the Spanish Ministry of Science. Computations were performed on the Teragrid (allocation TG-MCA02N014), the LONI network (http://www.loni.org), at LRZ München, the Barcelona Supercomputing Center, and at the Albert-Einstein-Institut.