Strong-field scattering of two spinning black holes: Numerical Relativity versus post-Minkowskian gravity

Highly accurate models of the gravitational-wave signal from coalescing compact binaries are built by completing analytical computations of the binary dynamics with non-perturbative information from numerical relativity (NR) simulations. In this paper we present four sets of NR simulations of equal-mass black hole binaries that undergo strong-field scattering: (i) we reproduce and extend the nonspinning simulations first presented in [Damour \textit{et al.}, Phys.Rev.D 89 (2014) 8, 081503], (ii) we compute two suites of nonspinning simulations at higher energies, probing stronger field interactions, (iii) we present a series of \textit{spinning} simulations including, for the first time, unequal-spin configurations. When comparing the NR scattering angles to analytical predictions based on state-of-the-art post-Minkowskian (PM) calculations, we find that PM-expanded scattering angles show poor convergence towards NR data. By contrast, a resummed computation of scattering angles via a spin-dependent, radiation-reacted, effective-one-body potential shows excellent agreement for both nonspinning and spinning configurations.


I. INTRODUCTION
The detection and characterisation of gravitationalwave (GW) observations [1][2][3][4][5] from compact binary coalescences relies on theoretical predictions of the emitted signal.Highly accurate models of coalescing compact binaries are a crucial requirement to perform precise measurements of the properties of black holes (BHs) and neutron stars, to determine their underlying astrophysical distributions and to perform tests of General Relativity in the strong-field regime.Such waveform models are generally built using both analytical and numerical approximations to Einstein's equations for a binary system in quasi-circular orbits.
The aims of this paper are: (i) to extend the results of Ref. [106] by performing higher-initial-energy NR simulations of equal-mass nonspinning BBH scattering 1 ; (ii) to perform NR simulations of the scattering of equal-mass spinning BBHs, both for equal and unequal spins, aligned with the angular momentum; and (iii) to compare the soobtained strong-field numerical scattering angles to analytical predictions based on state-of-the-art PM results 2 .
We denote the masses of the two objects as m 1 and m 2 , and the individual dimensionless spins as χ i = S i /(Gm 2 i ) = a i /(Gm i ), i = 1, 2. In addition, we denote M = m 1 + m 2 , µ = m 1 m 2 /(m 1 + m 2 ), ν = µ/M = m 1 m 2 /(m 1 + m 2 ) 2 .All our simulations will have m 1 = m 2 , so that ν = 1/4.Throughout this paper, we denote the canonical orbital angular momentum simply as L. It is related to the total Arnowitt-Deser-Misner (ADM) angular momentum of the system, J, by (1.1) Unless otherwise stated, we take G = c = 1 and generally work with dimensionless quantities.

II. NUMERICAL RELATIVITY (NR) SIMULATIONS
We performed several sequences of nonspinning and spin-aligned NR simulations modelling the scattering of two equal mass BHs.The numerical simulations were performed using the Einstein Toolkit (ETK) [109], an open source numerical relativity code built on the Cactus framework.The numerical setup of our simulations is broadly similar to that used in [106], though we present the details here for completeness.
We use Bowen-York initial data [110,111] computed using the TwoPunctures thorn [112].As in [106], the BHs are initially placed on the x-axis at a separation of ±X and with initial ADM linear momenta (2.1) where P ADM = | ⃗ P | and b NR denotes an impact parameter.As in [106], we use X = 50M .The Bowen-York initial data determine both the initial total ADM energy E ADM in of the system, and the total initial ADM angular momentum of the system.The latter is given by the simple formula where the initial canonical orbital angular momentum L in is related to b NR and P ADM by [106] L in = 2X|P y | = P ADM b NR . (2.3) We adopt as initial lapse profile α = ψ −2 BL , where ψ BL denotes the Brill-Lindquist conformal factor [7,111,113].In Appendix A, we present a check that the uncertainty linked to the choice of initial lapse profile has a subdominant effect on the inferred scattering angle θ NR compared to the errors in the polynomial fits used to measure the scattering angle.
Time evolution is performed using the W -variant [114] of the BSSNOK formulation [115][116][117] of the Einstein field equations as implemented by the McLachlan [118] thorn.We evolve the BHs using moving punctures gauge conditions [7,119], the lapse is evolved using the 1 + log condition [120], and the shift is evolved using the hyperbolic Γ-driver equation [113].We use eight-order accurate finite differencing stencils with Kreiss-Oliger dissipation [121].Adaptive mesh refinement is provided by Carpet, with the near zone being computed with highresolution Cartesian grids that track the motion of the BHs and the wave extraction zone being computed on spherical grids using the Llama multipatch infrastructure [122].The apparent horizons are computed using AHFinderDirect [123] and the spin angular momenta are calculated using the dynamical horizon formalism provided by the QuasiLocalMeasures thorn [124].
Similarly to [106], we found that the total energy and angular momentum of the system left after the release of the burst of spurious radiation present in the initial data differ from the corresponding ADM quantities computed from the initial data only by negligible fractions of order 10 −5 .
All ETK simulations were managed using Simulation Factory [125] and post-processing of NR data has made used of the open source Mathematica package Simulation Tools [126].

A. Extracting the Scattering Angle
In order to calculate the scattering angle, we follow the prescription detailed in [106].The motion of the BHs is tracked by the Cartesian coordinates of the punctures in the center-of-mass frame.We convert the tracks to polar coordinates (r, φ) (in the x − y plane of motion) for each BH and treat the incoming φ in,i (r) and outgoing φ out,i (r) paths for the i-th BH separately.We fit each of the paths to a polynomial of order n in terms of u = 1/r and extrapolate to find the asymptotic angle φ ∞ in/out,i .The resulting scattering angle for the i-th particle is given by Our choice of fitting window for (u, φ) follows [106], and we assume r ∈ [14,80]M and r ∈ [20,100]M for the incoming and outgoing trajectories respectively.We implement a least-squares fitting method that uses a singular value decomposition (SVD) to drop singular values smaller than 10 −13 times the maximum singular value, as was used in [106].In order to gauge the errors on the scattering angle, we perform a number of sanity checks.The first check is to gauge the impact of the polynomial order n on the resulting scattering angle.As in [106], the preferred polynomial order is taken to be the lowest order for which the SVD method allows for variation in the constant term.The extrapolation error is then estimated by finding the maximum and minimum scattering angle inferred over all polynomial orders between 2 and n − 1.This generally leads to dissymmetric error bounds on the scattering angle: θ +δ+ −δ− .A value of 0 in δ + or δ − occurs when the preferred polynomial is setting the bound.More precisely, δ − = 0 (respectively δ + = 0) occurs when the lower order polynomial fits over-estimate (respectively, under-estimate) the scattering angle relative to the preferred order.When least-square fitting our numerical results to analytical templates, we shall use a symmetrized version of the error bounds, namely Other sanity checks on the robustness of the derived scattering angle include testing the impact of the fitting window, testing the resolution of our numerical simulations, and varying gauge choices in the numerical evolution.We find that the inferred errors on the scattering angle remain subdominant to the choice of polynomial order, in agreement with the conclusions in [106].See Appendix A for further details.

B. Nonspinning Scattering Simulations
Though the main aim of the present work is to simulate the scattering of spinning BBHs, we also performed systematic sequences of nonspinning simulations of BBHs with fixed initial energy and varying initial angular momenta in order to be able to extract direct numerical estimates of the EOB interaction potential (see Sec. V below).For the first sequence, we reproduced and extended the results of Ref. [106], with an initial (adimensionalized) ADM energy ÊADM ≃ 1.02264 and varying initial (rescaled) orbital angular momentum Lin ≡ L in /M 2 .This allowed us both to cross-check the accuracy of our computations and to probe more precisely the boundary between scattering and plunge.We also performed two more sequences at the higher incoming energies, ÊADM in,2 ≃ 1.04033 and ÊADM in,3 ≃ 1.05548.These explore a stronger-field region of the parameter space, with NR impact parameters going down to bNR ≡ b NR /M ≃ 6 (with corresponding minimum isotropic EOB radial coordinate rmin ≃ 2, see below).
We report the results of all our nonspinning simulations in Table I.In Fig. 1 we show the scattering angles versus the impact parameter bNR for all three series of nonspinning systems.For the nonspinning simulations, we have that θ NR = θ 1 = θ 2 .

C. Spinning Scattering Simulations
We present here the results of 35 NR simulations of hyperbolic encounters of (equal-mass) spinning BBHs.As far as we know, it is the first time unequal-spin simulations have been published.In particular, we compute a suite of simulations at fixed energy and (canonical) orbital angular momentum and various spin magnitudes.The main parameters of these simulations are reported in Table II, while the scattering angles, as a function of the rescaled spin variables χ i are shown in Fig. 2.
When unequal spins are considered, the system ceases to be symmetric, inducing a non-zero recoil on the system center-of-mass.The effect of recoil on the scattering has been investigated in Ref. [86].As a consequence of Eq. (3.33) there, the scattering angles measured using the individual trajectories of the two bodies will be different.In the following, we will denote by θ NR the relative scattering angle, which, for equal-mass binaries, is simply given (to linear order in the recoil) by (2.5)

III. NR-DEDUCED PHYSICAL OBSERVABLES
In this section, we discuss some observables that can be directly extracted from the NR simulations presented in Sec.II above by using only minimal analytical assumptions.

A. Critical orbital angular momentum L0 in nonspinning scattering
We define the critical (canonical) angular momentum L 0 of nonspinning binaries as the value marking the boundary between scattering and plunging systems.Several works (both analytical [58,[127][128][129][130] and numerical [104,105,[131][132][133]) have investigated the value of L 0 (and the subtleties in defining it in view of the importance of radiation losses).
Here we follow Ref. [133] in estimating L 0 by fitting the sequence of NR scattering angles at fixed energy using a template expected to capture the singular behavior of θ(L) as L → L + 0 .Namely, we use (recalling the notation L ≡ L/M 2 ) ) where χ3PM ( L; L0 ) is a third order polynomial in 1/ L defined so as to ensure a correct PM expansion up to 3PM, and where a 4 is a 4PM-level fitting parameter (see Sec. VI A of Ref. [133] for details).
For the lower-energy simulations ( Ê1 = 1.02264), we could improve the estimate of L 0 computed in Ref. [133] and get Lfit 0 ( Ê1 ) = 1.0787 ± 0.00028 , a fit 4 ( Ê1 ) = 6.78 ± 0.50 , ( The corresponding estimate of Lfit 0 ( Ê1 ) in [133] is Lfit 0 ( Ê1 ) = 1.0773 .Instead, for the higher-energy ones ( Ê2 = 1.04033 and Ê3 = 1.05548, respectively corresponding to center-ofmass velocities v cm,2 ≃ 0.2757 and v cm,3 ≃ 0.3199), we got In Fig. 3 we display these numerical estimates against v cm .Note that the two new values Lfit 0 ( Ê2 ) and Lfit 0 ( Ê3 ) fill a gap in the previous knowledge of the critical L0 for intermediate energies.Figure 3 also displays the critical angular momentum estimates obtained through the high-velocity simulations of Refs.[104,105].We compare these to analytical estimates based on using the EOB transcription of PM scattering results (w eob nPM).For more details on these, see Sec.IV of Ref. [133].

B. Spin dependence of the scattering angle
We now focus on the spinning simulations and try to determine from our NR results (with minimal analytical assumptions) the spin dependence of the scattering angle at fixed initial energy and angular momentum.The color coding of Fig. 2 makes it apparent that the scattering angle mainly depends on the sum of the spins, i.e. is nearly constant along the lines χ 1 + χ 2 = cst.In particular, the simulations with opposite spins (χ 1 = −χ 2 , i.e. the second diagonal in Fig. 2) all lead to scattering angles θ NR very close (within the percent level) to the nonspinning one at the center.This numerical result corresponds to the well-known analytical fact (see, e.g., [134]) that the leading-PN value of the linear-in-spin interaction potential (spin-orbit) depends on the effective spin S eff = [1+3m 2 /(4m 1 )]S 1 +[1+3m 1 /(4m 2 )]S 2 , while the leading-PN value of the spin-spin interaction potential depends on In the equal-mass case considered here, this means that the two leading-order spin-dependent interactions only depend on the total spin S tot = S 1 + S 2 , and therefore only on χ 1 + χ 2 in the spin-aligned case.We have found that this property remains true to a very good approximation, much beyond the leading-order PN approximation.As we will describe below, we have derived a high-PM-accuracy Hamiltonian incorporating the state of the art analytical knowledge of PM gravity, and we found (as displayed in Table V in Appendix C) that the averaged scattering angle of two (equal-mass) spinning BHs is analytically predicted to depend essentially only on the average spin χ + ≡ 1 2 (χ 1 + χ 2 ).Let us also recall that the spin-orbit interaction term is repulsive (respectively, attractive) for spins parallel (respectively, antiparallel) to the orbital angular momentum L. This fact has a strong effect on the last stable orbits of quasi-circular spinning binaries (e.g., [134,135]).In the present, scattering situation, it implies that the leadingorder spin-dependent contribution to the scattering angle is negative (positive) when L • S eff > 0 (L • S eff < 0), see e.g.[31].This leading-order PN analytical prediction holds also true for the predictions from the PM-accurate EOB Hamiltonian derived below, as can be seen in the results displayed in Table VI in Appendix C. Our numerical results (displayed in Fig. 2) agree remarkably well with the analytical expectation that the (average) scattering angle depends essentially on the average spin, with only a weak dependence on the antisymmetric spin combination In particular, the numerical results displayed in Table II show that, for most of the (half) spin differences explored by our simulations, the dependence of θ NR ≡ 1 2 (θ 1 + θ 2 ) on χ − stays within the numerical uncertainty with which we could determine θ NR .For instance, even when χ − = 0.80, Table II reports that the difference between θ NR (χ 1 = 0, χ 2 = 0) = 221.823+0.762 −0.002 and θ NR (χ 1 = 0.80, χ 2 = −0.80= 221.679+0.489  −0.002 is only 0.144 ± 0.454.In view of this finding, the present limited precision of our numerical results does not allow us to meaningfully extract any reliable estimate of the dependence of the function θ NR (χ + , χ − ) on χ − .[We note in passing that, for the equal-mass systems we are considering, the relative (i.e.average) scattering angle ) is a symmetric function of χ 1 and χ 2 , and therefore an even function of χ − .] On the other hand, if we focus (for maximum accuracy) on the equal-spin subset of our numerical results (χ 1 = χ 2 , so that χ + = χ 1 = χ 2 and χ − = 0), we can fit our results to various possible analytical templates so as to extract a numerical estimate of the univariate function θ NR,eq (χ + ) = θ NR (χ + , χ − = 0).Removing the (χ 1 = χ 2 = −0.17)data point (which stood out as an outlier compared to its neighbouring data points in all our fits; see also Fig. 6), we found that we could fit within numerical errors3 the remaining eighteen equal-spin data points to the following simple five-parameter template: Note that this template incorporates a logarithmic singularity at some critical spin χ crit + .This singularity models the image on the χ + axis of the logarithmic singularity generically expected when crossing the border between scattering and coalescence [133].We find the best-fit values (when measuring angles in radians) to be leading to a satisfactory reduced chi-squared: χ 2 /(18 − 5) ≃ 0.63.
The successful fit of our numerical results to the template (3.8) allows us to meaningfully extract from our NR results an estimate of the critical value of χ + (at least in the equal-spin case) leading to immediate plunge rather than scattering.Indeed, the value χ crit + = −0.2932± 0.0013 formally corresponds to an infinite scattering angle.This result completes the critical values of L 0 displayed in Fig. 3 by yielding the numerical estimate of one point on the surface L crit 0 /E 2 = F ( Ê, χ 1 , χ 2 ) describing the critical initial data for spinning BBHs leading to immediate coalescence, namely the point (3.10)The computations above have been based on the relative scattering angle θ rel = 1 2 (θ 1 + θ 2 ).We leave to future work a NR-based study of the dissymmetry between θ 1 and θ 2 linked to asymmetric-spin recoil effects.

A. Post-Minkowskian expanded scattering angles
In the nonspinning limit, the PM approximation (which is a power series in the gravitational constant G) translates to a power series in 1/ℓ, where ℓ denotes the rescaled orbital angular momentum, The total (nonspinning) scattering angle up to nPM order (included) can then be written as where γ denotes the Lorentz factor between the two incoming worldlines.We recall that γ is equal to the rescaled EOB effective energy, γ = E eob eff /µ, which is related to the (real) total energy of the system by When considering spinning BHs, the PM expansion is modified because of the spin-induced multipolar structure of BHs.More precisely, the multipolar structure of Kerr black holes is polynomial in the Kerr parameters a i = S i /m i (with dimension of length), so that each spin-order, O(S m ) involves a m-th order polynomial in a 1 and a 2 .In the current PM literature on spinning bodies, it is usual to consider the Kerr parameters (or "ring radii") a i , as expansion parameters independent of the basic PM expansion parameters Gm 1 , Gm 2 .In other words, one uses a double expansion in powers of Gm i (PM expansion), and in powers of a i (spin expansion).This leads to a double counting (n, m), in which the PM order n counts the (sum of the) powers of Gm i , while m counts the (sum of the) powers of a i .However, one must remember that we are interested in the dynamics of spinning BH's for which we have the inequality a i ≤ Gm i .Therefore, a term of order (Gm) n a m is of real PM order n + m.This issue must be kept in mind, and will come back in our discussion below, but, for compatibility with the literature, we will continue to refer to a term of order (Gm) n a m as being formally of n-PM order.
The natural dimensionless expansion parameters associated with the PM and spin expansions are Gm i /b and a i /b where b denotes the impact parameter.We will, instead of b, use the (rescaled) angular momentum ℓ ≡ L/(Gm 1 m 2 ), with L = P cm b.We also utilize S i = m i a i as spin variables, so that each spin order will be identified by a corresponding power of S i /ℓ.One then writes the combined PM, and spin, expansion of the scattering angle as (in the spin-aligned case) with the kPM contribution to the scattering angle, θ k (γ, ℓ, S i ), being further expanded in spin powers, namely Here, θ S m k (γ) denotes a homogeneous polynomial in the two spins of order m.For instance As several theoretical papers on spinning scattering use as basic variables covariantly-defined impact parameters, spins and orbital angular momentum (say L cov ), while we work here with the canonically-defined orbital angular momentum (L ≡ L can ), we need to recall the connection between L cov and L ≡ L can .In the spin-aligned case, it is simply given by using the relation [138] where a i = S i /m i and E i = m 2 i + P 2 c.m. .In order to obtain the coefficients θ S m k of the double PM-and spinexpansion as functions of the canonical orbital angular momentum, we transform L cov into L can using Eq.(4.8), re-expand in powers of spin, and neglect higher-order contributions.
In the following, when considering spinning systems, we include: (i) nonspinning contributions up to 4PM (including radiation-reaction effects) [65,66,69,81,86,87]; (ii) quadratic-in-spin terms up to the 3PM order (including radiation-reaction effects) [145]; (iii) cubic and fourth order in spin effects up to the 2PM order [60].In addition, as during the development of this work the conservative [151], and radiative [152], linear-in-spin contributions have been computed at the 4PM level, we also incorporated these new analytical results.Their impact is separately discussed below.

B. w eob resummed scattering angles
In addition to PM-expanded scattering angles, we also compute w eob -resummed PM angles, as introduced in Ref. [133] and generalized here to spinning bodies.
Let us briefly introduce the EOB formalism (see, e.g., Sec.IIB of Ref. [133] for more details).The EOB framework [17,18,153] is a way of mapping the generalrelativistic dynamics of two masses m 1 , m 2 (considered in the center-of-mass system) onto the effective relativistic dynamics of a single body of mass µ = m 1 m 2 /(m 1 + m 2 ) = νM .The "real" center-of-mass Hamiltonian is related to the "effective" one through For scattering motions, the value of the effective energy E eff = H eff is simply related to the Lorentz factor γ = −u 1 • u 2 between the two incoming worldlines by E eff = µγ.This implies that the total ADM energy of the system is related to γ via The effective Hamiltonian H eff is obtained by solving (in P 0 = −H eff ) some relativistic-like mass-shell condition of the general form where Q(R, P µ ) gathers terms more than quadratic in P µ .When considering scattering problems, it is convenient to use a mass-shell condition where the effective metric is the Schwarzschild metric, and where Q(R, P µ ) is expressed as a function of the effective energy −P 0 = µγ.Using isotropic coordinates, and rescaled variables, r = R/M , p α = P α /µ, the PM-expansion of the EOB mass-shell condition formally takes the form of a nonrelativistic energy-conservation law, namely [58,154] where p ∞ = γ 2 − 1 and where the Newtonian-looking potential w eob (r; γ) [more precisely, −w eob (r)] encapsulates the (relativistic, energy-dependent) attractive gravitational interaction.For nonspinning systems, the PM expansion of the w eob potential reads where each order in 1 r corresponds to the same PM order (e.g., the first term, , with w orb 1 (γ) = 2(2γ 2 − 1), describes the 1PM gravitational interaction [57]).
For aligned-spin binaries we can generalize Eq. (4.13) simply by considering a spin-dependent radial potential of the form (up to fourth order in the spins) Here each w S m nPM (r, γ) is a polynomial of order m in spins, whose coefficients admit an expansion in powers of 1 r up to the order 1 rn included.For instance, we have (at spin orders m = 1 and m = 2) rk .
(4.16)The energy-dependent w * k (γ) coefficients entering the EOB potentials are in 1-to-1 correspondence with the θ * k (γ) coefficients of the scattering angle expansion.These relations can be found by expanding, and solving, at each order in G and in spin, the following equation . (4.17) Here rmin denotes the turning point radius.The explicit equations linking the θ * k (γ) coefficients and the w * k (γ) ones are given in Appendix B. When using the existing PM-expanded results for the scattering angle, the relations of Appendix B allow one to compute the explicit values of the PM-expansion coefficients w * k (γ) describing the spin interactions within the EOB formalism.These values will be found in the ancillary file of this paper.Note that the so-constructed EOB potential is a function of the incoming energy and angular momentum of the system that incorporates both conservative and radiation-reaction effects.

A. Numerics/analytics comparison for nonspinning simulations
We start by comparing, in Fig. 4, the results of the nonspinning NR simulations at fixed energy to both the PMexpanded, perturbative, predictions, Eq. (4.4), and their EOB-resummed versions, defined by Eq. (4.18).This comparison extends the one considered in Refs.[95,133], which corresponded to the lower initial energy Ê1 used in [106] (see also Appendix C).The main findings of the latter comparison when considering the PM-expanded perturbative predictions remain true for the new, higher energies considered in this paper, Ê2 = 1.04033 and Ê3 = 1.05548.The sequence of PM-expanded results (dashed lines) monotonically approach, in an increasingly accurate way, the NR results in the high angular momenta domain (i.e., for relatively large impact parameters fields).However, they perform less and less well as the angular momentum decreases.Note, in particular, that, for our minimum L, even the 4PM-accurate prediction differs by a factor ≳ 2 from the NR result.
On the other hand, the sequence of w eob -resummed (PM-based) predicted angles (solid lines) perform much better than their perturbative counterparts.While the 1PM-based w eob -resummed prediction is not better (and even slightly worse) than its perturbative analog, the 2PM-based w eob -resummed fares as well as the 4PM perturbative one, and the 3PM and, especially, 4PM EOBresummed results exhibit a remarkable agreement with NR data for essentially all orbital angular momenta.It should, however, be noted that the 4PM-based EOBresummed predicted angles start to show visible discrepancies with NR results for the lowest L's, especially for the highest energy Ê3 = 1.05548.More precisely, w eob 4PM somewhat overestimates the scattering angle, while w eob 3PM underestimates it.We discuss below how these findings can lead to NR-based ways of improving the knowledge of the exact w eob (r, γ) potential.

B. Extracting from NR results the gravitational interaction potential between nonspinning black holes
Having computed new suites of numerical simulations with higher initial energies, we can extend the strategy used in Sec.VIC of Ref. [133] and directly extract from the sequence of NR scattering angles (at a given energy) the NR radial potential w NR (r) (in our chosen EOB coordinates).This extraction is done by iteratively using Firsov's inversion formula [155,156] to invert the functional relation w(r, γ) → θ(ℓ, γ), Eq. (4.17) (for nonspinning systems), i.e. to deduce the function w NR (r) leading to the function θ NR (ℓ), tabulated in Table I (for each fixed value of γ).
The so-extracted NR radial potential for the lower value γ 1 (E 1 ) ≃ 1.09188 has been displayed in Fig. 7 of Ref. [133].See Appendix C for an updated determination of w NR (r, γ 1 ).We focus here on the determination of the NR potentials w NR (r, γ 2 ) and w NR (r, γ 3 ) corresponding to the two higher-energy sequences of simulations reported in the present work.In addition, we compare the so-extracted NR gravitational potentials to their analytical counterparts, as defined by the PM-based EOB radial potential w eob (r, γ) defined in Eq. (4.13) above (for the nonspinning case).
Instead of plotting the various energy-dependent radial potentials w * (r, γ), it is useful to plot the corresponding "effective radial potentials", completed by the additional centrifugal potential + ℓ 2 r2 .In other words, we shall compare the (Newtonian-like) effective potentials made by adding the (repulsive) centrifugal barrier term ℓ 2 r−2 to the (attractive) gravitational radial potential −w * (r, γ).The comparison between V NR (r, γ, ℓ) and V eob nPM (r, γ, ℓ) (for n = 3 and n = 4) is displayed in the two panels of Fig. 5.The left panel corresponds to the intermediate energy γ 2 (E 2 ) ≃ 1.16456, and uses the smallest NR orbital angular momentum leading to scattering, namely L ≃ 1.15950 (corresponding to bNR = 7.73).The right panel corresponds to the higher energy γ 3 (E 3 ) ≃ 1.22808, and uses the smallest NR orbital angular momentum leading to scattering for that energy, namely L ≃ 1.22500 (corresponding to bNR = 7.00).The NR-extracted potentials are determined only for (isotropic) radii larger than about 1.95 (for E 2 ) and about 1.85 (for E 3 ).
Note that, for both energies, the NR-extracted potential is very close to the 4PM-based EOB potential down to radii r ≃ 2.8.However, in the interval 2 ≲ r ≲ 2.8 V eob 4PM (r) starts to exhibit visible differences with V NR (r).For the intermediate energy E 2 the peak of V eob 4PM lies just below the p 2 ∞ line, meaning that the EOB-4PM potential predicts a plunge instead of a scattering for this angular  FIG. 5. Comparison between the gravitational potential V extracted from NR simulations and its EOB-PM equivalent.The Schwarzschild potential is plotted as a reference.Left panel: Intermediate energy Êin ≃ 1.04033.Right panel: Higher energy Êin ≃ 1.05548.In both cases, the angular momentum is set to the smallest values of angular momentum for which the system scattered.The horizontal lines mark the values of p 2 ∞ for each energy, determining the maximum point up to which it is possible to extract information about the numerical potential.
momentum and this energy.By contrast, the (less attractive) EOB-3PM potential still predicts scattering for this angular momentum and this energy.The situation is similar for the highest energy, E 3 .In this case, the NR potential seems to interpolate between V eob 4PM (r) and V eob 3PM (r).
These results suggest that it would be useful to resum the EOB-PM potentials 4 , so as to improve their performances at small radii.Leaving such an investigation to future work, we shall content ourselves here by exploring 4 We recall that we use here PM-expanded EOB radial potentials w(r) the effect of adding a formal, additional 5PM-level contribution + w5(γ) to the analytically predicted 4PM EOB potential.We can use our determination of w NR (r) to extract a best-fit value for the (energy-dependent) coefficient w 5 (γ).Fitting (for our three energies) w NR (r, γ) to w eob 4PM (r) + w5(γ) r5 yields the following results:  The negative sign (present at all energies) reflects the fact that the EOB-4PM potential is slightly too attractive at small radii (while the EOB-3PM potential is too repulsive for all radii).The strange behaviour at small radii, for the highest energy (right panel of Fig. 5) indicates that the addition of a (repulsive) 5PM term is only physically valid in a small range of radii (say in the range 2 ≲ r ≲ 2.8).Indeed, when one uses a non-resummed potential, w(r) = w1 r + w2 r2 + • • • wn rn , the last term determines the strong-field behavior.In particular, a negative 5PM term entails a (probably unphysical) repulsive core at very small radii.The use of suitably resummed versions of the EOB-PM potentials would probably avoid such undesirable features.

C. Comparing spinning simulations to analytical predictions (without 4PM spin-orbit terms)
We now shift our attention to the suite of 35 NR simulations which computed the scattering angles of systems having many different spins, but with fixed initial energy ( Êin ≃ 1.02264) and orbital angular momentum ( Lin = 1.14560).As we have already seen in Sec.III that the scattering angle mainly depends on the average spin χ + , we actually focus on the 19 simulations with equal spins (varying between χ 1 = χ 2 = −0.30and χ 1 = χ 2 = +0.80)listed in Table II.
In Fig. 6, we compare the equal-mass, equal-spin simulations to two types of PM-informed predictions: (i) standard PM-expanded angles, of the type of Eq. (4.4); and (ii) w eob -resummed ones [see Eq. (4.18)].In this section, we do not take into account the recently derived linearin-spin (spin-orbit) 4PM results of Refs.[151,152].In other words, we only use at the 4PM level (conservative and radiative) orbital effects.Since the analytical PM knowledge of the spin-dependent couplings varies with the spin order, we shall consider, for simplicity, only the following cases: (i) 1PM accuracy in all couplings up to quartic-in-spin terms (labelled as 1PM); (ii) 2PM accuracy in all couplings up to S 4 (labelled as 2PM); (iii) 3PM accuracy in the (S 0 , S 1 , S 2 ) couplings and 2PM accuracy in the (S 3 , S 4 ) ones (labelled as 3PM); and (iv) respectively (4PM,3PM,3PM,2PM,2PM) accuracy in the (S 0 , S 1 , S 2 , S 3 , S 4 ) couplings (labelled as 4PM).
A first conclusion one can draw from Fig. 6 is that PM-expanded results are rather unsatisfactory.Both the spin-averaged baseline and the spin-dependence are inaccurate, even for the highest PM accuracy.The agreement becomes, however, more satisfactory in the high-positivespin domain, which involves (because of the repulsive effect of parallel spins) weaker-field interactions.
By contrast, the w eob -resummed angles show, at each PM order, a systematically better agreement with NR data points, especially for the 4PM accuracy which exhibits a remarkable agreement.

D. Inclusion of 4PM spin-orbit terms
As already mentioned, conservative and radiative 4PM-level linear-in-spin contributions have been computed [151,152] during the development of this work.We transformed the latter scattering-angle results into additional contributions to the 4PM EOB potential w S 4PM .In Fig. 7 we compare the effect of incorporating the radiatively corrected 4PM spin-orbit terms ("4PM so"), to the 4PM-level angles computed in the previous section which only included orbital effects at 4PM ("4PM no so").For simplicity, we do not separately consider here the conservative contribution to the 4PM spin-orbit terms.
Though the inclusion of 4PM spin-orbit terms has a relatively minor effect in the PM-expanded scattering angles (see dashed lines), their inclusion in w S 4PM leads to much larger differences with respect to the corresponding 4PM-level EOB-resummed angles displayed in Fig. 6.
Looking at Fig. 7, it is clear that PM-expanded results remain rather unsatisfactory, and are significantly less accurate than w eob -resummed ones.However, the inclusion of 4PM-level spin-orbit effects in w eob worsens the EOB-NR agreement that was obtained when only including 4PM-level orbital coefficients.
It is possible that the worsening due to the inclusion of 4PM spin-orbit terms is a signal pointing out the need to introduce some resummation of the spin-dependent w S eob potential.[We recall that we are using here PMexpanded w S eob potentials, as indicated in Eq. (4.14).]However, the main cause of the worsening effect of the inclusion of 4PM spin-orbit terms might be the fact, recalled in Sec.IV above, that the nomenclature used in current spin-dependent PM works is physically inappropriate when dealing with spinning BHs.Indeed, because of the Kerr BH limit a i ≤ Gm i , a term of formal n-PM order including the m-th power of spins is actually of (n + m)-PM order.Therefore, the additional spin-orbit contributions computed in [151,152] are actually at the physical 5PM order.Similarly, the formally 2PM-level terms quartic in spins [140] are actually at the physical 6PM order.
In view of the latter remark, we have explored the effect of including an additional spin-dependent term in w eob belonging to the formal 5PM level, but to the physical 6PM level.For simplicity, we only considered a term involving a linear-in-spin contribution of the form This additional term contains two dimensionless parameters: (i) w 5 , parametrizing an orbital 5PM-level contribution, and (ii) w SO  5 , parametrizing a spin-orbit contribution at the physical 6PM level 5 .
In view of our above finding that the purely orbital part of the EOB potential is slightly too attractive and needs to be completed by a w5(γ) r5 term, where w 5 (γ) is an energy-dependent negative coefficient, we fix its value to 5 For simplicity, we do not include here a purely orbital 6PM-level contribution w 6 /r 6 .The horizontal line marks the "energy level" p 2 ∞ defining the turning point of the EOB mass-shell condition, Eq. (4.12).
w fit 5 (γ 1 ), Eq. (5.2), as the incoming energy of our spinning simulations corresponds to Ê1 ≡ Ê(γ 1 ).Concerning the second, spin-orbit-related, parameter w SO 5 in Eq. ( 5.3), we determine it by imposing that the corresponding effective Newtonian-like potential V eob , Eq. (5.1), predicts immediate coalescence for the critical value χ crit = −0.2932± 0.0013, Eq. (3.9), obtained above by fitting numerical results to the simple template of Eq. (3.8).This fixes the two parameters entering Eq. (5.3) to the approximate values w 5 = −1.00± 0.17 , w SO 5 = ± 0.6 . ( The corresponding w eob -resummed scattering angles computed using the (central) values of Eq. (5.4) are displayed in Fig. 7 as a red curve (labelled "4PM so + w 5 + w SO 5 ").The addition of the term ∆w eob , Eq. ( 5.3), leads to a better visual agreement between NR and EOB, at least for negative spins (on the large positive spin side, it slightly worsens the NR-EOB difference).
In order to better understand the physical origin of the differences between the various 4PM-level curves displayed in Fig. 7, we plot in Fig. 8 the corresponding Newtonian-like effective potentials V eob , Eq. (5.1).Let us recall the meaning of the labels in Fig. 7: (i) "4PM no so" refers to a potential including (radiation-reacted) 4PM orbital effects without including any spin-orbit contributions at 4PM; (ii) "4PM so" additionally includes 4PM radiatively-corrected spin-orbit terms [152]; and finally (iii) "4PM so + w 5 +w SO 5 " differs from "4PM so" by including the potential ∆w eob , Eq. ( 5.3), with the (central) values of Eq. (5.4).The dash-dotted horizontal line corresponds to the value of p 2 ∞ , which defines the turning point of the EOB mass-shell condition, Eq. (4.12).
The vertical distance between "4PM no so" (blue line) and "4PM so" (orange line) measures the effect of the (radiation-reacted) 4PM spin-orbit contribution.The inclusion of these terms leads to a more attractive potential.Actually, the corresponding potential is now below the p 2 ∞ "energy level" which means that the system is predicted to no longer scatter but to plunge (in agreement in Fig. 7).
The ∆w eob correction makes the potential (red line) more repulsive.This leads to a turning point around r ≃ 2.5 and to a finite scattering angle.This corrected potential also includes a (probably unphysical) repulsive core around r ≃ 1.
At this stage, the meaning of these results is unclear.One will need to explore more values of the incoming energy, of the incoming angular momentum, and of the spins to assess whether it is indeed necessary to add (energy-dependent) terms, such as Eq. ( 5.3), to the current best analytical knowledge of the EOB potential.

VI. DISCUSSION
In this paper we presented numerical simulations of the scattering of equal-mass BBHs using the Einstein Toolkit [109].
We first computed three sequences of equal-mass, nonspinning BBHs at fixed initial energy and varying angular momentum.The first sequence, with initial energy Êin,1 ≃ 1.02264, reproduces and extends the results of Ref. [106].These served as a useful cross-check against previous simulations performed with different numerical codes.They were also used to improve the estimate of the critical angular momentum L 0 , marking the boundary between scattering and plunging systems.We then computed two other sequences of simulations at higher energies, Êin,2 ≃ 1.04033 and Êin,3 ≃ 1.05548, probing stronger-field regimes.These two series allowed us to compute the critical angular momentum L 0 in previously unexplored regions of center-of-mass velocities, v cm,2 ≃ 0.2757 and v cm,3 ≃ 0.3199.This leaves as regions not yet covered by NR simulations of nonspinning BHs the low velocity regime, v cm ≲ 0.2, and the intermediate velocity one, v cm ≃ 0.3 − 0.6.
We then presented, for the first time, a suite of numerical simulations of equal-mass aligned-spins BHs.For these, we fixed both initial energy and angular momentum to be ( Êin , Lin ) ≃ (1.02264, 1.1456) and varied the individual spins in the range −0.3 < χ i < 0.8.We found (as expected both from NR simulations of coalescing binaries and from analytical predictions) the spin interaction to depend mainly on the total spin.We extracted from numerical data the coefficients of the linear and quadratic spin dependence of the scattering angle at the considered energy and angular momentum.
We compared our numerical results to PM-based predictions for the scattering angles of equal-mass BHs.We used two types of PM-based analytical predic-tions: (i) usual, perturbative PM-expanded results, see Eq. (4.4); and (ii) a transcription of PM results within the EOB formalism using a (spin-dependent, radiation-reacted) radial potential in isotropic (EOB) gauge, w eob (r, L, S 1 , S 2 ), see Eq. (4.18).We found that PM-expanded scattering angles exhibit (both for nonspinning and spinning systems) an acceptable agreement with NR data only for large impact parameters, consistently with the PM framework being a weak-field expansion.On the contrary, a transcription of the PM information into a corresponding (energy-dependent, radiationreacted, spin-dependent) EOB gravitational potential w eob (r, L, S 1 , S 2 ) yields remarkable improvements in the scattering angle comparisons, especially when using the highest available PM accuracy.See Figs. 4 and 6.In the case of spinning simulations, we found that the inclusion of the recently computed 4PM-level spin-orbit contribution worsened the remarkable agreement (displayed in Fig. 6) between NR angles and the EOB-resummed 4PMlevel ones (without 4PM spin-orbit terms).See Fig. 7.
Visible discrepancies between NR results and EOB-PM-predicted ones occur only for the few nonspinning cases where the impact parameter is close to the critical one leading to plunge rather than scattering and, for spinning simulations, when including the 4PM-level spinorbit contributions.See Fig. 5, where we compare (for nonspinning systems) the EOB-PM potential to a corresponding potential directly extracted from the NR scattering data by inverting the relation: potential → scattering angle.We could, however, improve the accuracy of the EOB-PM gravitational potential by adding (an NRfitted) additional (energy-dependent) 5PM-level contribution w fit 5 (γ)/r 5 , see Eq. (5.2).Similarly, in the spinning case, we could improve the NR-EOB match obtained when incorporating the 4PM-level radiation-reacted spinorbit terms by including an additional (formally) 5PMlevel contribution including a spin-orbit term, Eqs.(5.3) and (5.4).The sensitivity of the dynamics to spindependent contributions at the 4PM level is studied in Fig. 8.
The present work can be extended in several directions.On the numerical side, a more thorough exploration of the parameter space is called for, notably by considering different energies, mass ratios, and spin orientations.In addition, it would be useful to complete the computation of the scattering angles by accurately estimating the emitted waveform, and by numerically estimating the corresponding radiative losses of energy, linear momentum and angular momentum.
On the analytical side, one should explore possible resummations of the EOB-PM radial potential, and generalize it to the non-aligned-spin case.It would also be interesting to explore other possible EOB reformulations of perturbative PM results, e.g. using different EOB gauges.Finally, when having in hands NR waveforms, one should compare them to the recent analytical computations of O(G 3 )-accurate, PM-expanded waveforms [157][158][159].
Despite the above-mentioned limitations, we hope that the results presented here will offer a useful starting point for constructing accurate waveform models for spinning BHs in eccentric and hyperbolic motions, by combining information from NR and PM data within an EOB framework.

FIG. 9.
Error on the average scattering angle θ NR,av derived from the numerical relativity simulations for three different resolutions on the finest grid.Here, N = 80 represents the highest resolution, and default choice, whereas N = 64 denotes the lowest resolution.For all three resolutions, the upper errors are set by lower-order polynomial fits that overestimate the scattering angle.The lower error bars are set by high-order polynomials that manifestly agree with the preferred order.

Impact of Gauge Conditions
In order to help understand the impact of gauge conditions on our calculation of the scattering angle, we re-run a canonical test simulation, taken to be b NR = 10.0 and P ADM = 0.11456439, with an alternative lapse profile.The default lapse profile, as noted in earlier, is α = ψ −2 BL .For the alternative prescription, we adopt the lapse profile proposed in [160], which makes use of an approximate helical Killing vector field to minimise the dynamical evolution of the lapse.Writing the lapse in a pre-collapsed form, we have [160] where m i and r i denote the mass and location of the i-th puncture respectively.The lapse is then averaged, α = (1 + α)/2, to ensure that it lies within a range of [0, 1].The initial shift is taken to be β = 0.For this robustness test, we adopt a nonspinning, equal mass binary with P ADM = 0.11456439 and b NR , as shown in Fig. 10.We find that the choice of lapse again has a minor impact on the calculated scattering angle, finding θ NR,av = 221.823+0.762 −0.002 for the default lapse and θ NR,av = 222.089+0.605 −0.002 for the alternative lapse.The differences are smaller than the errors arising from the polynomial fitting procedure.

Appendix B: Mapping between scattering angle and potential coefficients
We here report the relations between the w i coefficients entering the radial potential w nPM , Eq. (4.14), and the θ i coefficients constituting the scattering angle θ nPM , Eq. (4.4).
We start by recalling the equations linking the various energy variables to the γ factor of the system: The relations between the nonspinning coefficients are known, and read In the following, we list the relations containing spin terms relevant to this paper: linear relations up to 4PM, quadratic-in-spin up to 3PM, cubic and quartic-in-spin at 1PM and 2PM.At linear order in spins we get while the quadratic-in-spin terms read The cubic-in-spin and quadratic-in-spin contributions at 1PM order are simply while the 2PM coefficients can be obtained by  4) for the lower-energy simulations of equal-mass nonspinning BBHs.FIG.
12. Gravitational potential V as a function of (centerof-mass, isotropic coordinates) separation (as in Fig. 5) for the lower energy simulations of equal-mass nonspinning BBHs.
The explicit values of the w i coefficients are given in the ancillary file of this paper.

Appendix C: Post-Minkowskian data points
In this appendix, we complement the PM results presented in the main paper.
In Figs.11 and 12, we mimic Figs. 4 and 5 for the nonspinning systems with lower energy, Êin,1 ≃ 1.02264.These figures confirm the outstanding agreement between EOB-PM predictions and NR results found in Ref. [106].
In Tables IV and Tables V, we list the scattering angles for the presented NR simulations together with the corresponding best-performing analytical predictions (not considering 4PM spin-orbit terms), θ w 3PM and θ w 4PM .
Finally, we compare the spin dependence of the NR scattering angles to PM predictions in the small spin range.Using the simplifying fact that we are considering equal-mass systems, and the relative scattering angle θ rel = 1 2 (θ 1 + θ 2 ), we can write, at quadratic order in spins, the simple general template We can, in principle, determine the θ i coefficients for the initial energy and angular momentum of our simulations, ( Êin , Lin ) ≃ (1.02264, 1.1456) by fitting the scattering angles for configurations with small (absolute) values of the spins.We decided to focus on the 6 configurations with |χ 1/2 | ≤ 0.1 and to fit them to the template Eq. (C1) (truncated to quadratic terms).We so obtained the following values: The reduced chi-squared of this fit is χ 2 /(6 − 4) ≃ 1.17 (corresponding to six data points and four degrees of freedom).This value is probably affected by the fact that we are assuming constant initial data, while the energy slightly changes between simulations.Let us emphasize that we could not extract a meaningful value for the θ fit χ 2 − coefficient, consistently with our finding that these effects are subdominant (even for higher spin values).
In Table VI we compare these NR-fitting coefficients to their (PM-expanded and w eob -resummed) analytical analogs.In order to extract w eob 3PM and w eob 4PM (not including 4PM spin-orbit terms), we repeat the same fitting procedure employed for the NR angles (but setting the initial energy to the constant, average, value).For the considered initial data, ( Êin , Lin ) ≃ (1.02264, 1.1456), the PM-expanded predictions are not satisfactory, as could be deduced from Fig. 6.The successive PM orders show a slow convergence towards the numerical results.The EOB-resummed equivalents, instead, are in good agreement with the numerical estimates, with the w eob 4PM values mostly agreeing within one standard deviation with the NR ones.

FIG. 2 .
FIG. 2.NR scattering angles for equal-mass simulations with Êin ≃ 1.02264 and Lin ≃ 1.1456 against individual spin variables.This figure shows (see the color code) that the angles essentially depend only on the sum of the spins.

3 FIG. 3 .
FIG. 3. Comparison between the critical (rescaled) orbital angular momentum L0/E 2 extracted from NR simulations and various analytical predictions.The three leftmost NR data points correspond to our fitted estimates [Eqs.(3.2)-(3.4)].The first point updates the estimate of Ref. [106], while the other two (highlighted by a vertical grey line) explore a new velocity interval.The five points on the right correspond to the high-energy simulations of Refs.[104, 105].We also display, for comparison, some EOB-PM predictions.

FIG. 6 .
FIG.6.Scattering angle comparison between NR data and both the perturbative, and the EOB-transcribed, PM-based predictions for equal-mass, equal-spin simulations.The x-axis represents χ eff = χ1 = χ2.As in Fig.4, we use dashed lines for the PM-expanded scattering angles, while the solid lines correspond to the w eob -resummed ones.Similarly to the nonspinning case, the w eob -resummed predictions are markedly more accurate, especially the 4PM-informed ones.

FIG. 7 . 6 .
FIG.7.Effects due to the inclusion of the 4PM spin-orbit contributions[151,152] to the scattering angle comparison of Fig.6.Dashed lines correspond to PM-expanded scattering angles, while solid lines represent w eob -resummed predictions.The blue curves contain spin-dependent terms up to 3PM and nonspinning contributions at 4PM (as in Fig.6above).The orange curves also incorporate (radiation-reacted) spin-orbit effects at 4PM level.Finally, the red line includes the two additional parameters w5 and w SO 5 of Eqs.(5.3)-(5.4).

5 FIG. 8 .
FIG.8.Comparison of various EOB-PM gravitational potentials V eob in the most negative spin case (χ eff = −0.25),i.e. the leftmost NR point in previous figures.The displayed analytical potentials correspond to the solid lines in Fig.7.The horizontal line marks the "energy level" p 2 ∞ defining the turning point of the EOB mass-shell condition, Eq. (4.12).

FIG. 10 .
FIG.10.Trajectories of each BH using two different gauge choices for the initial lapse.The first choice is α = ψ −2 BL , and the second the averaged lapse defined in Eq. (A1).The dots denote the starting position of each BH.

TABLE I .
Equal-mass, nonspinning NR simulations.The scattering angles have been replaced by dots where the BHs have plunged.An asterisk denotes simulations where there is uncertainty on an eventual plunge.Longer simulations are required to determine whether the system becomes bound after the first encounter.bNR ÊADM

TABLE II
Comparison between NR scattering angles of nonspinning BBH systems and analytical PM predictions, both PMexpanded (dashed lines) and w eob -resummed (solid lines).Left panel: Intermediate energy Êin ≃ 1.04033.Right panel: Higher energy Êin ≃ 1.05548.The w eob -resummed 3PM and 4PM results show a remarkable agreement with numerical down to strong-field regimes.

TABLE III .
A set of equal-mass, nonspinning simulations with an initial momentum P ADM = 0.175 and different resolutions. )

TABLE IV .
NR equal-mass, nonspinning simulations presented in TableI.The last two columns list the w eobresummed scattering angles at 3PM and 4PM as in Figs.4 and 11.The scattering angles are replaced by dots when the black holes (are predicted to) plunge.An asterisk denotes uncertainty about an eventual plunge.

TABLE V .
NR equal-mass, spinning simulations listed in TableII.All systems were computed with fixed initial ADM energy ÊADM in ≃ 1.02264 and orbital ADM angular momentum LADM in ≃ 1.14560.The last two columns list the w eobresummed scattering angles at 3PM and 4PM as in Fig.6.

TABLE VI .
Small-spin dependence of the scattering angle for Êin ≃ 1.02264 and Lin ≃ 1.1456.PM predictions for the different spin-orders (if available) together with values extracted by fitting NR data [Eq.(C2)].