High-accuracy mass, spin, and recoil predictions of generic black-hole merger remnants

Vijay Varma,

∗

Davide Gerosa,

†

François Hébert,

‡

Leo C. Stein,

1, 2,

and Hao Zhang

1, 3,

1

TAPIR 350-17, California Institute of Technology, 1200 E California Boulevard, Pasadena, CA 91125, USA

2

Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA

3

Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

(Dated: September 26, 2018)

We present accurate fits for the remnant properties of generically precessing binary black holes,

trained on large banks of numerical-relativity simulations. We use Gaussian process regression to

interpolate the remnant mass, spin, and recoil velocity in the 7-dimensional parameter space of

precessing black-hole binaries. For precessing systems, our errors in estimating the remnant mass,

spin magnitude, and kick magnitude are lower than those of existing fitting formulae by at least an

order of magnitude. In addition, we also model the remnant spin and kick directions. Improvement

is also reported for aligned-spin systems. Being trained directly on precessing simulations, our fits

are free from ambiguities regarding the initial frequency at which precessing quantities are defined.

As a byproduct, we also provide error estimates for all fitted quantities, which can be consistently

incorporated into current and future gravitational-wave parameter-estimation analyses. Our model(s)

are made publicly available through a fast and easy-to-use Python module called

surfinBH

Introduction

– As two black holes (BHs) come together and

merge, they emit copious gravitational waves (GWs) and

leave behind a BH remnant. The strong-field dynamics

of this process are analytically intractable and must be

simulated using numerical relativity (NR). However, from

very far away, the merger can be viewed as a scattering

problem, depicted in Fig. 1. The complicated dynamics

of the near zone can be overlooked in favor of the gauge-

invariant observables of the in- and out-states: the initial

BH masses and spins, the outgoing GWs, and the final

BH remnant. This final BH is fully characterized by its

mass, spin, and recoil velocity; all additional complexities

(“hair”) of the merging binary are dissipated away in

GWs [1–3].

All GW models designed to capture the entire inspiral-

merger-ringdown (IMR) signal from BH binary coales-

cences need to be calibrated to NR simulations (e.g.,

[

–

]). In particular, the BH ringdown emission is cru-

cially dependent on the properties of the BH remnant

— properties obtained from NR simulations. Accurate

modeling of the merger remnant is therefore vital for

construction of accurate IMR templates.

Besides waveform building, accurate knowledge of the

remnant properties is also instrumental to fulfill one of

the greatest promises of GW astronomy: testing Ein-

stein’s general relativity (GR) in its strong-field, highly

dynamical regime. Current approaches to test the Kerr

hypothesis attempt to measure the properties of the inspi-

ralling BHs from the low frequency part of the GW signal,

then use NR fits to predict the corresponding remnant

mass and spin; this final-state prediction is compared to

the properties inferred from the high frequency part of the

∗

vvarma@caltech.edu

†

Einstein Fellow; dgerosa@caltech.edu

‡

fhebert@caltech.edu

lcstein@olemiss.edu

zhangphy@sas.upenn.edu

Interaction

region

f

f

f

Gravitational waves to

I

FIG. 1. Quasi-circular binary BH merger problem viewed as

a scattering process via a “Feynman” diagram. Time flows to

the right. All quantities are well defined in the asymptotically

flat region far from the interaction (merger).

GW signal [

]. Inaccuracies in remnant models there-

fore directly propagate to the final fundamental-physics

test.

The importance of building fits for the remnant proper-

ties was realized soon after the NR breakthrough [

–

]

and has been periodically revisited by several groups since

then [

–

]. There are two important shortcomings in all

existing fitting formulae. First, they enforce analytic an-

sätze (with NR-calibrated coefficients) that are physically

motivated, but lack a rigorous mathematical justification.

Therefore, current fits can be prone to systematic errors,

especially in regions of parameter space where the intu-

ition used to design the formulae become less accurate.

Second, current expressions for remnant mass and spins

are calibrated on aligned-spin simulations and therefore

fail to fully capture the rich physics of precessing systems

(but see e.g. [

] where a non-generic subspace of precess-

ing configurations is considered). For example, current

LIGO/Virgo parameter-estimation pipelines [

] rely

on ad-hoc corrections to partially account for precession

effects [

]. Aligned fits applied to precessing systems

are inevitably ambiguous, as the outcome will depend on

where

(in time, separation, or frequency) the spins are

defined and inserted into the fits (e.g., [43]).

In this

Letter

we tackle both these issues for the first

arXiv:1809.09125v1 [gr-qc] 24 Sep 2018

time. We construct surrogate models that fit the rem-

nant properties from a large sample of generic, precess-

ing, quasi-circular binary BH simulations performed with

the Spectral Einstein Code (SpEC) [

]. Surrogates are

trained directly against the NR simulations, using Gaus-

sian process regression (GPR) without any phenomenolog-

ical ansätz, and achieve accuracies comparable to those

of the NR simulations themselves. In their regime of

validity, the models presented here are at least an order

of magnitude more accurate than previous fits.

In particular, we present two models:

surfinBH7dq2

: a fit trained against precessing sys-

tems with mass ratios

≤

and dimensionless spin

magnitudes

,χ

≤

surfinBH3dq8

: an aligned-spin model trained

against systems with mass ratios up to

≤

and

(anti-)aligned spin magnitudes

,χ

≤

Both these models can be easily accessed using the pub-

licly available Python module

surfinBH

[

], and are

ready to be incorporated in both waveform constructions

and GW parameter-estimation studies.

Fitting procedure

– We construct fits for the BH rem-

nant mass

, spin vector

, and recoil kick vector

as functions of the binary mass ratio

and spin vec-

tors

. Our fits for

surfinBH7dq2

(

surfinBH3dq8

)

map a 7-

(3-)

dimensional input parameter space to a 7-

(4-)

dimensional output parameter space. The fits are per-

formed in the coorbital frame at

−

100

, with

= 0

the peak of the total waveform amplitude (cf. Ref. [

]

for details). The coorbital frame is defined such that

the

-axis lies along the direction of the orbital angular

momentum, the

-axis runs from the smaller BH to the

larger BH, and the

-axis completes the triad.

All fits are performed using GPR [

]; details are pro-

vided in the supplemental material. Notably, GPR natu-

rally returns estimates of the errors of the fitted quantities

across the parameter space.

The values of spins, masses, and kicks used in the

training process are extracted directly from the NR sim-

ulations. We use the simulations presented in Ref. [

]

for

surfinBH7dq2

and those of Ref. [

] for

surfinBH3dq8

Both spins and masses are evaluated on apparent hori-

zons [

]; the dimensionful spin

solves an eigenvalue

problem for an approximate Killing vector, and the mass

is determined from the spin and area

following the

Christodoulou relation

irr

)

, where

irr

is the irreducible mass. The masses

are determined close to the beginning of the simula-

tion at the “relaxation time” [

], whereas the spins

≡

are measured at

−

100

. The

remnant mass

and spin

are determined long after

ringdown, as detailed in [

]. All masses are in units of

the total mass

at relaxation. The remnant

kick velocity is derived from conservation of momentum,

−

rad

[

]. The radiated momentum flux

rad

is integrated [

] from the GWs extrapolated to future

null infinity [

]. Before constructing the fits,

and

are transformed into the coorbital frame at

−

100

Besides the GPR error estimate, we further address the

accuracy of our procedure using “k-fold” cross validations

with

= 20

. First, we randomly divide our training

dataset into

mutually exclusive sets. For each set, we

construct the fits using the data in the other

−

sets and

then test the fits by evaluating them at the data points in

the considered set. We thus obtain “out-of-sample” errors

which conservatively indicate the (in)accuracies of our

fits. We compare these errors against the intrinsic error

present in the NR waveforms, estimated by comparing

the two highest resolutions available. We also compare

the performance of our fits against several existing fitting

formulae for remnant mass, spin, and kick which we denote

as follows: HBMR ([

] with

= 3

), UIB

[

], HL [

], HLZ [

], and CLZM ([

]

as summarized in [

]). To partially account for spin

precession, fits are corrected as described in Ref. [

] and

used in current LIGO/Virgo analyses [

]: spins are

evolved from relaxation to the Schwarzschild innermost

stable circular orbit, and final UIB and HL spins are

post-processed adding the sum of the in-plane spins in

quadrature.

Aligned-spin model

– We first present our fit

surfinBH3dq8

which is trained against 104 aligned-spin simulations [

]

with

≤

and

−

≤

,χ

≤

. Symmetry implies

that the kick lies in the orbital plane while the final spin

is orthogonal to it [

]. We therefore only fit for four

quantities:

, and

Figure 2 shows the out-of-sample errors of

surfinBH3dq8

. Our fits are as accurate as the NR

simulations used in the training process.

per-

centile errors lie at

∆

∼

−

∆

∼

−

, and

∆

∼

−

. The kick direction is predicted with an

accuracy of

∼

radians, which is the inherent accuracy

of the NR simulations. Our errors for the remnant mass

and kick magnitude are comparable to the most accurate

existing fits. On the other hand, for the final spin, our

procedure outperforms all other formulae by at least a

factor of

Precessing model

– We now present

surfinBH7dq2

, a rem-

nant model trained on 890 simulations [

] of generic,

fully precessing BH binaries with mass ratios

≤

and

spin magnitudes

,χ

≤

. Out-of-sample errors are

shown in Fig. 3.

percentiles are

∼

−

for

mass,

∼

−

for spin magnitude,

∼

−

radians

for spin direction,

∼

−

for kick magnitude, and

∼

radians for kick direction. As in the aligned-spin

case above, our errors are at the same level as the NR res-

olution error, thus showing that we are not limited by our

fitting procedure but rather by the quality of the training

dataset. Our fits appear to outperform the NR simula-

tions when estimating the spin direction, which suggests

this quantity has not fully converged in the NR runs.

Figure 3 shows that our procedure to predict remnant

mass, spin magnitude, and kick magnitude for precessing

∆

[

]

3dq8

HBMR

UIB

∆

[

001

]

3dq8

CLZM

HLZ

−

∆

−

cos

−

(

) [radians]

FIG. 2. Errors in predicting remnant mass, spin, kick magnitude, and kick direction for non-precessing binary BHs with mass

ratios

≤

, and spin magnitudes

,χ

≤

. The direction error is the angle between the predicted vector and a fiducial

vector, taken to be the high-resolution NR case and indicated by a

?

. The square markers indicate median values while the

triangle markers indicate

percentile values. Our model

surfinBH3dq8

is referred to as 3dq8. The black histogram shows the

NR resolution error while the dashed histograms show errors for different existing fitting formulae.

∆

[

]

7dq2

HBMR

UIB

∆

[

001

]

7dq2

CLZM

∆

cos

−

(

) [radians]

−

cos

−

(

) [radians]

−

cos

−

(

) [radians]

−

[

001

]

7dq2

3dq8

FIG. 3. Errors in predicting the remnant mass, spin magnitude, spin direction, kick magnitude, and kick direction for precessing

binary BHs with mass ratios

≤

, and spin magnitudes

,χ

≤

. The square markers indicate the median values while the

triangle markers indicate the

percentile values. Our model,

surfinBH7dq2

is referred to as 7dq2. The black histogram shows

the NR resolution error while the dashed histograms show errors for different existing fitting formulae. In the bottom-right

panel we show the distribution of kick magnitude vs. error in kick direction.

[

]

−

Error

−

∆

[

]

∆

cos

−

(

)

∆

[

001

]

cos

−

(

)

FIG. 4. Left panel: Errors for

surfinBH7dq2

in predicting

remnant properties when the spins are specified at an orbital

frequency of

=10

, for different total masses. Right panel:

Errors for

surfinBH7dq2

when extrapolating to higher mass

ratios, with the spins specified at

−

100

. The labels on the

horizontal axis indicate the range of mass ratios being tested.

Note that the distributions in these plots are normalized to

have a fixed height, not fixed area.

systems is more precise than all existing fits by at least

an order of magnitude. These existing fits presented sig-

nificantly lower errors when applied to aligned binaries

(cf. Fig. 2), which suggests that they fail to fully capture

precession effects despite the augmentation of Ref. [

Some impact of precession effects on the final spin and

recoil is expected, since both of these quantities have been

found to depend strongly on the in-plane orientations of

the spins of the merging BHs [

]. More surpris-

ingly, we find that spin precession significantly affects the

energy radiated as well, which was expected to depend

mostly on the aligned-spin components via the orbital

hang-up effect [55–57].

The largest errors in the kick direction can be of order

∼

radian. The bottom-right panel of Fig. 3 shows the

joint distribution of kick magnitude and kick direction

error for both

surfinBH7dq2

and

surfinBH3dq8

, showing

that errors are larger at low kick magnitudes. Our error in

kick direction is below

∼

radians whenever

&

−

Regime of validity

– The errors in Fig. 3 are obtained by

evaluating fits using input spins specified at

−

100

i.e., where the GPR interpolation is performed. The

input spins can also be specified at earlier times; this

case is handled by two additional layers of time evolution.

Given the spins at an initial orbital frequency

, we first

evolve the spins using a post-Newtonian (PN) approxi-

mant — 3.5PN SpinTaylorT4 [

–

]— until the orbital

frequency reaches a value of

018

rad

. At this point,

we are in the range of validity of the (more accurate)

NRSur7dq2 approximant [

], which we use to evolve

the spins until

−

100

. Thus, spins can be specified

at any given orbital frequency and are evolved consis-

tently before estimating the final BH properties. This is

a crucial improvement over previous results, which, being

calibrated solely to non-precessing systems, suffer from

ambiguities regarding the separation/frequency at which

spins are defined.

The left panel of Fig. 4 shows the errors when the spins

are specified at an orbital frequency

= 10

. These

errors are computed by comparing against 20 long NR

simulations [

] with mass ratios

≤

and generically

oriented spins with magnitudes

,χ

≤

. None of

these simulations were used to train the fits. Longer PN

evolutions are needed at lower total masses, and the errors

are therefore larger. These errors will decrease with an

improved spin evolution procedure. Note, however, that

our predictions are still more accurate (and, crucially,

unambiguous) than those of existing fitting formulae (cf.

Fig. 3).

Finally, the right panel of Fig. 4 shows the the per-

formance of

surfinBH7dq2

when extrapolating to more

extreme mass ratios. We compare against 175 (225) NR

simulations [

] with

≤

(

≤

), and generically

oriented spins with magnitudes

,χ

≤

specified at

−

100

. The error distribution broadens, but our fits

still provide a reasonable estimate of the final remnant

properties even far out of the training parameter space.

Detailed results on extrapolation accuracy are provided

in the supplemental materials.

Conclusion

– We have presented two highly accurate sur-

rogate models for the remnant properties of BH binaries.

surfinBH7dq2

(

surfinBH3dq8

) is trained against 890 (104)

NR simulations with mass ratios

≤

(

≤

) and pre-

cessing (aligned) spins with magnitude

,χ

≤

. Both

models use GPR to provide fits for the remnant mass,

spin, and kick velocity (both magnitudes and directions).

Our findings are implemented in a public Python module

named

surfinBH

(details are provided in the supplemental

materials).

For aligned spins, errors in

surfinBH3dq8

are com-

parable to existing fitting formulae for the final mass

and kick magnitude, while the spin is predicted about 5

times more accurately. For precessing systems, errors in

surfinBH7dq2

for final mass, spin magnitude, and kick

magnitude are lower than all existing models by at least

an order of magnitude. Crucially, our fits are free from

ambiguities regarding the time/frequency at which pre-

cessing quantities are specified. This is a point of major

improvement over previous models, which all fail to fully

capture precession effects.

Is this increased accuracy necessary? For current events

like GW150914, the estimated error in the remnant prop-

erties are

∆

∼

and

∆

∼

[

]. These mea-

surements are currently dominated by statistical errors,

as the systematics introduced by existing fits used in the

analysis are

∆

∼

−

and

∆

∼

−

(see

percentile values in Fig. 3). Because statistical errors

scale approximately linearly with the detector sensitivity

[

], we estimate that systematic errors in current models

for

will start dominating over statistical uncertainties

at signal-to-noise ratios which are

∼

times larger than

that of GW150914. This will happen sooner rather than

later, with current interferometers expected to reach their

design sensitivity in a few years [

], and future instru-

ments already being scheduled [

] or planned [

Our fits, being an order of magnitude more accurate (see

Fig. 3), introduce systematic errors which are expected

to be relevant only at SNRs

∼

times larger than that

of GW150914. As shown above, errors are largely domi-

nated by the underlying NR resolution, not by our fitting

procedure. The inclusion of self-force evolutions alongside

NR in the training dataset might also be exploited to

improve extrapolation performance at

; we leave

this to future work.

Moreover, the GPR methods employed here naturally

provide error estimates along with the fitted values (some

results are provided in the supplemental material). This

constitutes a further key application of our results: when

performing, e.g., consistency tests of GR [

], sys-

tematic uncertainties introduced by remnant fits can be

naturally incorporated into the statistical analysis and

marginalized over (cf. Ref. [

] for a similar application

of GPR).

As GW astrophysics turns into a mature field, increas-

ingly accurate tools such as those presented here will

become crucial to uncover more hidden secrets in this new

field of science.

Acknowledgments

– We thank Jonathan Blackman,

Stephen Taylor, David Keitel, Anuradha Gupta, and

Serguei Ossokine for useful discussions. We made use of

the public LIGO Algorithm Library [

] in the evaluation

of existing fitting formulae and to perform PN evolutions.

We thank Nathan Johnson-McDaniel for useful discus-

sions, comments on the manuscript, and for sharing his

code to evaluate the HLZ kick fits. V.V. and F.H. are

supported by the Sherman Fairchild Foundation and NSF

grants PHY–1404569, PHY–170212, and PHY–1708213

at Caltech. D.G. is supported by NASA through Einstein

Postdoctoral Fellowship Grant No. PF6–170152 awarded

by the Chandra X-ray Center, which is operated by the

Smithsonian Astrophysical Observatory for NASA under

Contract NAS8-03060. L.C.S. acknowledges support from

NSF grant PHY–1404569 and the Brinson Foundation.

H.Z. acknowledges support from the Caltech SURF Pro-

gram and NSF Grant No. PHY–1404569. Computations

were performed on NSF/NCSA Blue Waters under allo-

cation NSF PRAC–1713694 and on the Wheeler cluster

at Caltech, which is supported by the Sherman Fairchild

Foundation and by Caltech.

[1]

W. Israel, Communications in Mathematical Physics

8

245 (1968).

[2] B. Carter, PRL

26

, 331 (1971).

[3]

M. Heusler,

Black hole uniqueness theorems, Cambridge

University Press, 1996.

(1996).

[4]

M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa,

F. Ohme, G. Pratten, and M. Pürrer, PRL

113

, 151101

(2014), arXiv:1308.3271 [gr-qc].

[5]

S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer,

X. J. Forteza, and A. Bohé, PRD

93

, 044007 (2016),

arXiv:1508.07253 [gr-qc].

[6]

S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme,

X. J. Forteza, and A. Bohé, PRD

93

, 044006 (2016),

arXiv:1508.07250 [gr-qc].

[7]

A. Buonanno, Y. Pan, H. P. Pfeiffer, M. A. Scheel, L. T.

Buchman, and L. E. Kidder, PRD

79

, 124028 (2009),

arXiv:0902.0790 [gr-qc].

[8]

A. Bohé, L. Shao, A. Taracchini, A. Buonanno, S. Babak,

I. W. Harry, I. Hinder, S. Ossokine, M. Pürrer, V. Ray-

mond, T. Chu, H. Fong, P. Kumar, H. P. Pfeiffer,

M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace,

M. A. Scheel, and B. Szilágyi, PRD

95

, 044028 (2017),

arXiv:1611.03703 [gr-qc].

[9]

S. Babak, A. Taracchini, and A. Buonanno, PRD

95

024010 (2017), arXiv:1607.05661 [gr-qc].

[10]

S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, and

M. Tiglio, PRX

4

, 031006 (2014), arXiv:1308.3565 [gr-qc].

[11]

J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley,

D. A. Hemberger, P. Schmidt, and R. Smith, PRD

95

104023 (2017), arXiv:1701.00550 [gr-qc].

[12]

J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley,

C. D. Ott, M. Boyle, L. E. Kidder, H. P. Pfeiffer, and

B. Szilágyi, PRD

96

, 024058 (2017), arXiv:1705.07089

[gr-qc].

[13]

A. Ghosh, N. K. Johnson-McDaniel, A. Ghosh, C. Kant

Mishra, P. Ajith, W. Del Pozzo, C. P. L. Berry, A. B.

Nielsen, and L. London, CQG

35

, 014002 (2018),

arXiv:1704.06784 [gr-qc].

[14]

B. P. Abbott

et al.

(LIGO Scientific Collaboration

and Virgo Collaboration), PRL

116

, 221101 (2016),

arXiv:1602.03841 [gr-qc].

[15] F. Pretorius, PRL

95

, 121101 (2005), gr-qc/0507014.

[16]

M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo-

chower, PRL

96

, 111101 (2006), gr-qc/0511048.

[17]

M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D.

Matthews, and H. P. Pfeiffer, PRD

79

, 024003 (2009),

arXiv:0810.1767 [gr-qc].

[18]

F. Herrmann, I. Hinder, D. M. Shoemaker, P. La-

guna, and R. A. Matzner, PRD

76

, 084032 (2007),

arXiv:0706.2541 [gr-qc].

[19]

M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Mer-

ritt, PRL

98

, 231102 (2007), gr-qc/0702133.

[20]

J. A. González, M. Hannam, U. Sperhake, B. Brügmann,

and S. Husa, PRL

98

, 231101 (2007), gr-qc/0702052.

[21]

J. A. González, U. Sperhake, B. Brügmann, M. Hannam,

and S. Husa, PRL

98

, 091101 (2007), gr-qc/0610154.

[22]

M. Campanelli, C. Lousto, Y. Zlochower, and D. Merritt,

ApJ

659

, L5 (2007), gr-qc/0701164.

[23]

L. Rezzolla, E. Barausse, E. N. Dorband, D. Pollney,

C. Reisswig, J. Seiler, and S. Husa, PRD

78

, 044002

(2008), arXiv:0712.3541 [gr-qc].

[24]

L. Rezzolla, P. Diener, E. N. Dorband, D. Pollney, C. Reis-

swig, E. Schnetter, and J. Seiler, ApJ

674

, L29 (2008),

arXiv:0710.3345 [gr-qc].

[25] M. Kesden, PRD

78

, 084030 (2008), arXiv:0807.3043.

[26]

W. Tichy and P. Marronetti, PRD

78

, 081501 (2008),

arXiv:0807.2985 [gr-qc].

[27]

C. O. Lousto and Y. Zlochower, PRD

77

, 044028 (2008),

arXiv:0708.4048 [gr-qc].

[28]

E. Barausse and L. Rezzolla, ApJ

704

, L40 (2009),

arXiv:0904.2577 [gr-qc].

[29]

Y. Pan, A. Buonanno, M. Boyle, L. T. Buchman, L. E.

Kidder, H. P. Pfeiffer, and M. A. Scheel, PRD

84

, 124052

(2011), arXiv:1106.1021 [gr-qc].

[30]

E. Barausse, V. Morozova, and L. Rezzolla, ApJ

758

, 63

(2012), [Erratum: ApJ, 2014, 786, 76], arXiv:1206.3803

[gr-qc].

[31]

C. O. Lousto, Y. Zlochower, M. Dotti, and M. Volonteri,

PRD

85

, 084015 (2012), arXiv:1201.1923 [gr-qc].

[32]

C. O. Lousto and Y. Zlochower, PRD

87

, 084027 (2013),

arXiv:1211.7099 [gr-qc].

[33]

J. Healy, C. O. Lousto, and Y. Zlochower, PRD

90

104004 (2014), arXiv:1406.7295 [gr-qc].

[34]

Y. Zlochower and C. O. Lousto, PRD

92

, 024022 (2015),

arXiv:1503.07536 [gr-qc].

[35]

F. Hofmann, E. Barausse, and L. Rezzolla, ApJ

825

, L19

(2016), arXiv:1605.01938 [gr-qc].

[36]

D. Gerosa and M. Kesden, PRD

93

, 124066 (2016),

arXiv:1605.01067 [astro-ph.HE].

[37]

X. Jiménez-Forteza, D. Keitel, S. Husa, M. Hannam,

S. Khan, and M. Pürrer, PRD

95

, 064024 (2017),

arXiv:1611.00332 [gr-qc].

[38]

J. Healy and C. O. Lousto, PRD

95

, 024037 (2017),

arXiv:1610.09713 [gr-qc].

[39]

J. Healy and C. O. Lousto, PRD

97

, 084002 (2018),

arXiv:1801.08162 [gr-qc].

[40]

B. P. Abbott

et al.

(LIGO Scientific Collaboration

and Virgo Collaboration), PRX

6

, 041015 (2016),

arXiv:1606.04856 [gr-qc].

[41]

B. P. Abbott

et al.

(LIGO Scientific Collaboration and

Virgo Collaboration), PRL

118

, 221101 (2017), [Erratum

PRL, 2018, 21, 129901], arXiv:1706.01812 [gr-qc].

[42]

N. K. Johnson-McDaniel, A. Gupta, P. Ajith,

D. Keitel, O. Birnholtz, F. Ohme, and S. Husa,

dcc.ligo.org/T1600168/public.

[43]

M. Kesden, U. Sperhake, and E. Berti, PRD

81

, 084054

(2010), arXiv:1002.2643 [astro-ph.GA].

[44]

L. E. Kidder, M. A. Scheel, S. A. Teukolsky, E. D. Carlson,

and G. B. Cook, PRD

62

, 084032 (2000), gr-qc/0005056.

[45]

V. Varma

et al.

pypi.org/project/surfinBH,

doi.org/10.5281/zenodo.1418525.

[46]

C. E. Rasmussen and C. K. I. Williams,

Gaussian Pro-

cesses for Machine Learning, by C.E. Rasmussen and

C.K.I. Williams. ISBN-13 978-0-262-18253-9

(2006).

[47]

V. Varma, S. Field, M. A. Scheel, and J. Blackman,

(2018), in preparation.

[48]

G. Lovelace, R. Owen, H. P. Pfeiffer, and T. Chu, PRD

78

, 084017 (2008), arXiv:0805.4192 [gr-qc].

[49] M. Boyle

et al.

, (2018), in preparation.

[50]

D. Gerosa, F. Hébert, and L. C. Stein, PRD

97

, 104049

(2018), arXiv:1802.04276 [gr-qc].

[51]

M. Ruiz, M. Alcubierre, D. Núñez, and R. Takahashi,

General Relativity and Gravitation

40

, 1705 (2008),

arXiv:0707.4654 [gr-qc].

[52]

M. Boyle and A. H. Mroué, PRD

80

, 124045 (2009),

arXiv:0905.3177 [gr-qc].

[53]

L. Boyle, M. Kesden, and S. Nissanke, PRL

100

, 151101

(2008), arXiv:0709.0299 [gr-qc].

[54]

E. Berti, M. Kesden, and U. Sperhake, PRD

85

, 124049

(2012), arXiv:1203.2920 [astro-ph.HE].

[55]

M. Campanelli, C. O. Lousto, and Y. Zlochower, PRD

74

, 041501 (2006), gr-qc/0604012.

[56]

C. O. Lousto and Y. Zlochower, PRD

89

, 021501 (2014),

arXiv:1307.6237 [gr-qc].

[57]

M. A. Scheel, M. Giesler, D. A. Hemberger, G. Lovelace,

K. Kuper, M. Boyle, B. Szilágyi, and L. E. Kidder, CQG

32

, 105009 (2015), arXiv:1412.1803 [gr-qc].

[58]

A. Buonanno, Y. Chen, and M. Vallisneri, PRD

67

104025 (2003), gr-qc/0211087.

[59]

M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroué, H. P.

Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky,

PRD

76

, 124038 (2007), arXiv:0710.0158 [gr-qc].

[60]

S. Ossokine, M. Boyle, L. E. Kidder, H. P. Pfeiffer,

M. A. Scheel, and B. Szilágyi, PRD

92

, 104028 (2015),

arXiv:1502.01747 [gr-qc].

[61]

V. Varma, S. Field, M. A. Scheel,

et al.

, (2019), in

preparation.

[62] M. Vallisneri, PRD

77

, 042001 (2008), gr-qc/0703086.

[63]

B. P. Abbott

et al.

(VIRGO, KAGRA, LIGO Scientific),

LRR

21

, 3 (2018), arXiv:1304.0670 [gr-qc].

[64]

P. Amaro-Seoane

et al.

(LISA Core Team), (2017),

arXiv:1702.00786 [astro-ph.IM].

[65] M. Punturo

et al.

, CQG

27

, 194002 (2010).

[66]

B. P. Abbott

et al.

(LIGO Scientific Collaboration

and Virgo Collaboration), CQG

34

, 044001 (2017),

arXiv:1607.08697 [astro-ph.IM].

[67]

C. Cahillane, J. Betzwieser, D. A. Brown, E. Goetz, E. D.

Hall, K. Izumi, S. Kandhasamy, S. Karki, J. S. Kissel,

G. Mendell, R. L. Savage, D. Tuyenbayev, A. Urban,

A. Viets, M. Wade, and A. J. Weinstein, PRD

96

, 102001

(2017), arXiv:1708.03023 [astro-ph.IM].

[68]

LIGO Scientific Collaboration and Virgo Collaboration,

git.ligo.org/lscsoft/lalsuite.

Supplemental Material

Gaussian process regression

– We construct fits in this

work using Gaussian process regression (GPR) [

]

as implemented in

scikit-learn

[S3].

The starting point is a training set of

observations,

{

(

))

= 1

,...,n

}

, where

denotes an in-

put vector of dimension

and

(

)

is the corresponding

output. In our case,

is mass ratio and spins of the

merging binary, and

(

)

is the remnant property we are

fitting. Our goal is to use

to make predictions for the

underlying

(

)

at any point

∗

that is not in

A Gaussian process (GP) can be thought of as a prob-

ability distribution of functions. More formally, a GP

is a collection of random variables, any finite number of

which have a joint Gaussian distribution [

]. A GP is

completely specified by its mean function

(

)

and co-

variance function

(

x,x

′

)

, i.e.

(

)

∼GP

(

)

(

x,x

′

))

Consider a prediction set of

∗

test inputs and their

corresponding outputs (which are unknown):

{

(

∗

(

∗

))

= 1

,...,n

∗

}

. By the definition of a GP,

outputs of

and

(respectively

{

(

)

}

∗

{

(

∗

)

}

) are related by a joint Gaussian distribution

[

∗

]

(

[

∗

])

(S1)

where

∗

denotes the

∗

matrix of the covariance

(

x,x

∗

)

evaluated at all pairs of training and prediction

points, and similarly for the other

matrices.

Eq. (S1) provides the Bayesian prior distribution for

∗

The posterior distribution is obtained by restricting this

joint prior to contain only those functions which agree

with the observed data points, i.e. [S1]

(

∗

|TS

) =

(

∗

−

, K

∗

−

∗

−

∗

)

(S2)

The mean of this posterior provides an estimator for

(

)

∗

, while its width is the prediction error.

Finally, one needs to specify the covariance (or ker-

nel) function

(

x,x

′

)

. In this

Letter

we implement the

following kernel

(

x,x

′

) =

exp





−

∑

(

−

′

)





x,x

′

(S3)

where

x,x

′

is the Kronecker delta. In words, we use

a product between a squared exponential kernel and a

constant kernel, to which we add a white kernel term to

account for additional noise in the training data [

GPR fit construction involves determining the

+ 2

hyperparameters (

and

) which maximize the

marginal likelihood of the training data under the GP

prior [

]. Local maxima are avoided by repeating the

optimization with 10 different initial guesses, obtained by

sampling uniformly in log in the hyperparameter space

described below.

Before constructing the GPR fit, we pre-process the

training data as follows. We first subtract a linear fit and

the mean of the resulting values. Datapoints are then

normalized by dividing by the standard deviation of the

resulting values. The inverse of these transformations is

applied at the time of the fit evaluation.

For each dimension of

, we define

∆

to be the

range of the values of

and consider

∈

∆

]

. Larger length scales are unlikely

to be relevant and smaller length scales are unlikely to be

resolvable. The remaining hyperparameters are sampled

∈

[10

−

]

and

∈

[10

−

]

. These choices

are meant to be conservative and are based on prior ex-

ploration of the typical magnitude and noise level in our

pre-processed training data.

Input parameter space

– Fits for

surfinBH3dq8

are pa-

rameterized using

= [

log

(

)

χ,χ

]

, where

is the spin

parameter entering the GW phase at leading order [

–

]

in the PN expansion,

113(1 +

)(

qχ

)

−

(

)

113(1 +

)

−

(S4)

and

is the “anti-symmetric spin”,

(

−

)

(S5)

For

surfinBH7dq2

use

[

log

(

)

,χ

χ,χ

,χ

]

Subscripts

and

refer to components specified in the coorbital

frame at

−

100

. We empirically found these

parameterizations to perform more accurately than the

more intuitive choice

= [

q,χ

,χ

]

In the main text we describe how we evolve spins

given at earlier times to

−

100

, using PN and NR-

Sur7qd2. Is it worth noting that the NR spins used to

train

NRSur7qd2

had some additional smoothening fil-

ters applied to them (see Eq. 6 in [

]). This introduces

additional systematics when evolving spins from times

−

100

. We verified that the resulting errors on our

fits are subdominant.

Extrapolation erorrs

– The right panel of Fig. 4 shows

the errors in remnant quantities when extrapolating

surfinBH7dq2

to mass ratios beyond its training range

(

≤

). These errors are computed using the spins at

−

100

. If the spins are given at earlier times, we

expect larger extrapolation errors as this also involves

extrapolation of the NRSur7dq2 waveform model (which

was also trained in the

≤

space). Figure S1 shows the

extrapolation errors when the spins are specified at at

orbital frequency

= 10

for a total mass

= 70

computed by comparing against the same NR simulations

as in Fig. 4. Errors are comparable to or lower than those

of existing fits for

≤

. For

< q

≤

, our errors for

the remnant spin magnitude can become larger, but the