{'fraction_non_alphanumeric': 0.04101948602089805, 'fraction_numerical': 0.017897486585710252, 'mean_word_length': 4.299102132435466, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mapping the asymptotic inspiral of precessing binary black holes to their merger remnants

Multiple approaches are required to study the evolution of black-hole binaries. While the post-Newtonian approximation is sufficient to describe the early inspiral (even from infinitely large orbital separation), only numerical relativity can capture the full complexity of the dynamics near merger. We combine multi-timescale post-Newtonian integrations with numerical-relativity surrogate models, thus mapping the entire history of the binary from its asymptotic configuration at past-time infinity to the post-merger remnant. This approach naturally allows us to assess the impact of the precessional and orbital phase on the properties - mass, spin, and kick - of the merger remnant. These phases introduce a fundamental uncertainty when connecting the two extrema of the binary evolution.


Introduction
Binary black holes (BHs) are prominent sources of gravitational waves (GWs). Ten [1] or more [2,3] BH mergers have been detected during the first two LIGO/Virgo observing runs. Data from dozens of other candidates [4] resulting from the third run are currently being analyzed. With more interferometers becoming operational both on the ground [5,6] and in space [7], the observed population of BH binaries is expected to grow dramatically in the near future.
Some of the events observed so far show clear evidence for at least one BH with nonzero spin [8][9][10][11]. Systems where both BHs are spinning represent the most general and complex scenario: the three angular momenta of the binary (i.e. the spins of the two BHs and the orbital angular momentum) are all coupled to each other and precess about the total angular momentum of the system [12]. Meanwhile, the emitted GWs dissipate energy and angular momentum, causing the orbital separation to shrink [13] and ultimately driving the BHs to merger. GWs also carry linear momentum, which results in a recoil (or kick) imparted to the center of mass of the system [14]. The final Kerr BH left behind following a merger is fully characterized by its mass M f , spin χ f , and proper velocity v f .
Precessing BH binaries evolve on three different timescales. Let us denote the binary total mass with M , the orbital separation with r, and use geometrical units c = G = 1. The BHs orbit about each other with period t orb /M ∝ (r/M ) 3/2 , the angular momenta precess on a timescale t pre /M ∝ (r/M ) 5/2 [12], and GW radiationreaction takes place on the timescale t RR /M ∝ (r/M ) 4 [13]. In the early inspiral where the dynamics is well described by the post-Newtonian (PN) approximation, one has r M and therefore t orb t pre t RR . This timescale separation implies that BH binaries complete a very large number of orbits and precessional cycles before they merge. By the time they become detectable, BHs will have lost memory of the precessional and orbital phase with which they formed. Astrophysical formation models aimed at describing BH binaries from assumptions on their progenitors cannot predict the orbital and precessional phase with which GW sources enter the detector sensitivity band.
The same reasoning applies to the properties (final mass, spin and kick) of the merger remnant which are known to depend significantly on both the precessional [28,29] and the orbital dynamics [30,31]. These quantities are set by two phases, describing orbit and precession, which are essentially random.
The timescale hierarchy of spinning binary BHs has been exploited to develop orbit- [32,33] and precession-averaged [34,35] PN approaches. Together, these schemes are capable of evolving sources from their asymptotic configuration at r → ∞ (f → 0, where f is the GW frequency) down to the small separations where the PN approximation ceases to be an accurate physical description. From those small separations, estimates of the remnant parameters can be quickly obtained by evaluating surrogate models trained on numerical-relativity (NR) simulations [36][37][38]. These models have been shown to reach accuracies comparable to those of the underlying NR data. In particular, Ref. [38] recently presented a surrogate model for the spin dynamics and remnant properties of generically precessing binary BHs covering mass ratios from 1:1 to 1:4 and dimensionless spin magnitudes up to 0.8.
In this paper, we combine the multi-timescale approach of Ref. [34,35] with the NR surrogate model of Ref. [38]. We design a procedure to put together precessionaveraged, orbit-averaged, and NR-surrogate evolutions. This allows us to connect the asymptotic parameter space of precessing BH binaries to the properties of the merger remnant, thus mapping the entire general relativistic two-body problem from past to future time infinity. We highlight how, for a given astrophysical configuration, the predictions of the remnant properties suffer from theoretical uncertainties due to the sampling of the orbital and precessional phases.
The paper is organized as follows. In Sec. 2 we illustrate our procedure for numerically evolving BH binaries; in particular we describe the implementation of multiple layers of PN/surrogate techniques as well as the random sampling of the phases. In Sec. 3 we present our results. We first focus on individual sources and assess the impact of phase sampling on the merger remnant properties. We then consider BH binary populations to determine the impact of spin precession on the statistical distributions of the remnant parameters. We draw our conclusions in Sec. 4.

Multiple layers of black-hole evolution
Let us consider spinning BH binaries on quasi-circular orbits. Let m 1 and m 2 be the masses of the two BHs, M = m 1 + m 2 the total mass, q = m 2 /m 1 ≤ 1 the mass ratio, η = m 1 m 2 /M 2 the symmetric mass ratio, S i = m 2 i χ i (i = 1, 2) the BH spin magnitudes (with 0 ≤ χ i ≤ 1), and L = m 1 m 2 r/M the magnitude of the Newtonian orbital angular momentum. We write vectors in boldface, denote the corresponding unit vectors by hats, and their magnitudes as e.g. L = |L|. Let us also define the following angles describing the spin directions

Early inspiral
Much like the inequality t orb t pre can be exploited to average over the orbital motion, the relation t pre t RR allows averaging over spin precession. The computational cost of orbit-averaged integrations blows up as r/M → ∞ where infinitely many precession cycles need to be tracked. On the contrary, precession-averaged evolutions are capable of modeling the binary inspiral from its asymptotic conditions at past time infinity. The cost one pays for precession (orbital) averaging is the ability to track the precessional (orbital) phase.
The relative orientations of the three angular momenta L, S 1 , and S 2 at a given separation r (or equivalently L) are fully specified by the following three parameters: which is conserved on both t RR and t pre at 2PN order [39,33] ; • the magnitude of the total angular momentum which is constant on t pre but slowly varies on t RR ; • the magnitude of the total spin which oscillates between two extrema S + ≡ max ξ,J (S) and S − ≡ min ξ,J (S) on the shorter timescale t pre .
The BH binary inspiral can thus be modelled by averaging over S (the only parameter that varies on t pre ) and solving a single ordinary differential equation for J (the only other parameter that varies on on t RR ). In practice, the parameter S sets the phase of the coupled precessional motion of L, S 1 , and S 2 .
Since J ∼ L ∝ √ r diverges at large separations, it is convenient to introduce the auxiliary variable κ ≡ J 2 − L 2 2L (5) * This same parameter is more often indicated as χ eff in the GW literature.
which instead asymptotes to a constant The integration domain can be compactified by setting u = 1/2L such that initializing an evolution at u = 0 corresponds to integrating over the entire past history of the binary. Performing a precession-averaged integration at 1.5PN reduces to solving [35] (see also [40]) where dS/dt is obtained directly from the 2PN orbit-averaged spin-precession equations, cf. Eq. (26) in Ref. [35]. This scheme is implemented in the public code precession [41]. A BH binary spin configuration at r/M → ∞ can be equivalently described using either ξ and κ ∞ , or the asymptotic projections of the two spins along the orbital angular momentum

Late inspiral
As shown extensively in Ref. [35], precession-and orbit-averaged PN techniques are in very good agreement for separations r 50M . In the late inspiral, the timescales t pre and t RR become more comparable and the precession-averaged approach loses accuracy. We transition from precession-to orbit-averaged integrations at an orbital separation r pre→orb = 1000M . This value is well within the validity of the precession-averaged approach and allows keeping the computational cost of the orbit-averaged evolution under control. The robustness of this choice is investigated in Sec. 3.5.
Initializing an orbit-averaged evolution requires the mutual orientations of the spins and the orbital angular momentum, which can be specified by either (θ 1 , θ 2 , ∆Φ) or equivalently (ξ, J, S). Only two of these parameters are provided by the previous precession-averaged treatment which averages over S. We estimate a value S pre→orb of S at the transition threshold statistically by drawing a random sample from the probability distribution function p(S) ∝ |dS/dt| −1 between S + and S − . Intuitively, it is more (less) likely to find the binary with a precessional phase S if the "velocity" dS/dt is smaller (larger).

Plunge and merger
The mass, spin and recoil velocity of the merger remnant are estimated using the NR surrogate models of Ref. [38]. This includes the NRSur7dq4 model for the spin dynamics, and the NRSur7dq4Remnant model for the remnant properties. Both models are trained on NR simulations with mass ratios q ≥ 0.25, spin magnitudes χ 1 , χ 2 ≤ 0.8, and generic spin directions. Surrogate remnant predictions outperforms alternative fitting formulae by orders of magnitude [37,38] and have recently been shown to be critical to model and extract BH kicks from GW signals [31,42].
The surrogate fits are provided as a function of the initial orbital frequency ω, mass ratio q, and six Cartesian components of the spins S 1 and S 2 . The orbit averaged scheme used in the previous step only captures the relative orientation of the spins and the orbital plane, which corresponds to five spin degrees of freedom (i.e. χ 1 , χ 2 , ξ, J, and S). The additional parameter needed is the orbital phase ϕ, describing the location of the BHs on their orbits with respect to the spins or, equivalently, the orientation of the in-plane components of the spins with respect to the vector connecting the two BHs.
The PN orbit-averaged evolution is carried over until the orbital frequency reaches ω orb→sur = 0.025 rad/M . This falls well within the validity of the NRSur7dq4 model, which only includes ∼ 20 orbits before the merger; † tests are reported in Sec. 3.5. The conversion between PN separation r and orbital frequency ω is performed using the 2PN expression [43] At the transition frequency, we sample a value ϕ orb→sur of the orbital phase uniformly in [0, 2π]. First, the NRSur7dq4 spin dynamics is evolved from the initial orbital frequency ω orb→sur to a time t = −100M before the waveform amplitude peak. Second, the NRSur7dq4Remnant model is evaluated to predict the remnant mass, spin, and recoil velocity. The Gaussian-process regression algorithm [44,37] implemented in Ref. [38] returns mean values M f , χ f and v f and 1-σ errors σ M f , σ χ f and σ v f of respectively the remnant mass, spin, and kick. Unless otherwise specified, we estimate these quantities by sampling values M f , χ f and v f from the corresponding normal distributions, thus incorporating the fitting errors of the NR surrogate model.

Summary
Our procedure is summarized in Fig. 1. We map the entire history of the binary from its asymptotic configuration at r/M → ∞ to its merger product. We start by initializing a precession-averaged PN evolution at u = 0 (r = ∞) for a binary with parameters q, χ 1 , χ 2 , θ 1∞ , and θ 2∞ . We integrate Eq. (7) from the initial configuration to an orbital separation r pre→orb = 1000M . At this separation we sample the precessional phase

i u A A H I J j 4 I J z U A P X o A 4 a A I N 7 8
A i e w Y v 1 Y D 1 Z r 9 b b r L V g 5 T P 7 4 I + s 9 y + 3 n J o 2 < / l a t e x i t > r = 1000M < l a t e x i t s h a 1 _ b a s e 6 4 = " In the following, we characterize the remnant properties using 90% intervals and define where x = {M f , χ f , v f } is any of the remnant properties and x n indicates the n-th percentile of its marginalized distribution.

Uncertainty of the black hole remnant
We start by considering a single binary at r/M → ∞ and studying the impact of the two phases.
In the examples below we selected a source with q = 0.6, χ 1 = χ 2 = 0.8, θ 1∞ = π/3, and θ 2∞ = π/5 but the trends we highlight are generic. We first evolve this same binary 10 5 times repeatedly sampling over all three sources of uncertainty: the precessional phase S pre→orb , the orbital phase ϕ orb→sur , and the surrogate Gaussian errors. Figure 2 shows the resulting distribution of the BH remnant properties. GWs dissipate energy and angular momentum relatively gradually during the inspiral; on the other hand, the vast majority of the linear momentum is emitted during the last few orbits and plunge [30,31]. The exact configuration with which the binary approaches merger, therefore, plays a more crucial role for which are even more extreme: possible values of the final recoil range from v f ∼ 0 to ∼ 0.01.

Contributions to the remnant uncertainty
We now dissect the various contributions entering these uncertainties. Figure 3 shows the remnant properties of the same asymptotic configuration considered previously, but now assuming a fixed value of the precessional phase S pre→orb = 0.34M 2 . We then vary ϕ orb→sur from 0 to 2π. The solid curves in Fig. 3 show the median values M f , χ f and v f returned by the NR surrogate; circles mark a sample of remnants extracted from the corresponding Gaussian distributions.
The orbital-phase uncertainty on the remnant properties is due to the breaking down of the PN timescale hierarchy. The series of binaries depicted in Fig. 3 is computed for different phases ϕ orb→sur and fixed orbital frequency ω orb→sur . For t orb /t pre → 0 and t orb /t RR → 0, this sequence is equivalent to the same binary moved along a quasi-adiabatic orbit, and should just result in identical remnants. Deviations from these limits imply that, although these sources are identical up to r pre→orb , they evolve into different configurations and thus result in different merger remnants.
The mean values of the remnant properties returned by the NR surrogate are displayed with a solid curve; circles indicate a statistical sample extracted from the corresponding Gaussian distributions.
ϕ orb→sur ranges from 0 to 2π. An approximate periodicity with period equal to π is also present, such that the curve in Fig. 3 appears as two, largely overlapping loops. A period equal to twice the orbital frequency corresponds to that of the dominant GW emission. The separation between the two loops is typically smaller than the Gaussian errors of the surrogate model, and thus largely disappears when these are considered (scatter points). For binaries with q = 1 and S 1 = S 2 , a transformation ϕ → ϕ + π can be absorbed into a change of reference frame and thus results in identical values of M f , χ f , and v f . The two-loop structure of Fig. 3  . Final remnant properties as a function of the precessional phase S pre→orb measured at the transition between precession-and orbit-averaged PN evolution. From left to right, we show the mass-spin, mass-kick and spin-kick planes. For this evolution we assume q = 0.6, χ 1 = χ 2 = 0.8, θ 1∞ = π/3, and θ 2∞ = π/5. Each row is generated similarly to Fig. 3, by keeping S pre→orb fixed while varying the orbital phase ϕ orb→sur from 0 to 2π. The value of S pre→orb increases from top to bottom from S − to S + as reported on the color bar. In particular, the panels are produced with S pre→orb 0.32, 0.34, 0.37, 0.39, 0.41.

thus collapses into a single loop (in practice, this is
only true within the errors of the surrogate model because this expected symmetry was not built in explicitly in Ref. [38]). For non-symmetric binaries, the π periodicity is broken resulting in a more complex phenomenology. In general, we find that for binaries with more extreme mass ratios q 0.25, the loop of the median values becomes tighter and is eventually "filled" as the surrogate errors are considered. We refer to Refs. [45,46] for a careful analysis of the expected symmetries. Figure 4 shows how the remnant properties vary with the precessional phase S pre→orb . While the precise shape of the curves depends on the sampled value, the general behavior we described remains valid: the remnant properties present an approximate period equal to that of the dominant GW mode, which is half the orbital period. The region covered by the loops mildly changes with S, such that the weighted superpositions of all the loops in Fig. 4 returns the full distribution we first showed in Fig 2. Each individual contribution is weighted by the probability density function p(S) ∝ |dS/dt| −1 , which is strongly peaked at the two extrema S − and S + [41]. The loops with the larger and smaller values of S pre→orb (top and bottom row in Fig. 4) provide the largest contributions to the overall distribution and can indeed be glimpsed when carefully inspecting Fig. 2.

Black-hole binary populations
We now consider BH binary populations and investigate the statistical properties of their remnants. In Fig. 5, we focus on a population with fixed q = 0.8 and χ 1 = 0.6 and χ 2 = 0.7. First, we map the entire (cos θ 1∞ , cos θ 2∞ ) plane into the (M f , χ f ) parameter space (gray region in the right panel of Fig. 5). Final mass and final spin are strongly anticorrelated, with large values of M f corresponding to small values of χ f , and vice-versa. This is due to the orbital hang-up effect [39,47,48]. Binaries with spins which are both aligned to the orbital angular momentum (θ 1∞ θ 2∞ 0, corresponding to larger positive values of ξ in the upper right corner of the left panel of Fig. 5) inspiral for longer: they dissipate more energy via GWs and result in lighter remnant. Their angular momenta add constructively (|L + S 1 + S 2 | L + S 1 + S 2 ) resulting in a larger final spin. On the contrary, binary with spins which are largely anti-aligned plunge from afar but the spins and angular momentum partially cancel each other, thus producing remnants which are heavier and more slowly rotating. We then consider a subsample of 25 equispaced sources from the same population, shown as coloured dots in the left panel of Fig. 5. In particular, binaries are coloured according to the projected effective spin ξ [Eq. (2)]. We evolve each of these systems 10 4 times, sampling over S pre→orb , ϕ orb→sur and the surrogate Gaussian errors. The corresponding distributions in the remnant parameter space are shown in the right panel of Fig. 5 with the same color scheme. The larger these regions, the more uncertain the remnant properties are. The size of the regions correlates with ξ: remnant uncertainties tend to be larger for |ξ|∼ 0 and smaller for values at the extrema |ξ|∼ (χ 1 + qχ 2 )/(1 + q). This is because the former systems have larger precessional amplitudes |S + − S − | while for the latter the range where S pre→orb can be sampled is more limited [35]. Figure 6 shows the uncertainties on the remnant properties for a wider set of BH binary populations. In particular, we consider three cases with fixed values q = 0.8 and χ 1 = χ 2 = 0.3, q = 0.6 and χ 1 = χ 2 = 0.8, q = 0.95 and χ 1 = χ 2 = 0.8, as well as a fourth scenario where we distribute the mass ratio according to p(q) ∝ q 6.7 in [0. 25,1] and the spin magnitudes uniformly in [0.1, 0.8]. The latter is motivated by current GW detections which favour equal-mass systems [49] (but see Ref. [10]). In all four cases, we consider both isotropically distributed spins (top panels) and partially aligned spins (bottom panels) at past time infinity, which correspond to distributing cos θ i∞ either uniformly in [−1, 1] or with a truncated normal distribution centered on 1 with dispersion 0.1, respectively. We extract 1000 samples from each of these populations. These asymptotic configurations are then evolved 1000 times, sampling the precessional and orbital phases, to obtain corresponding distributions of the remnant properties. We then construct the estimators δM f , δχ f and δv f of Eq. (11) which are shown in Fig. 6.
Uncertainties on M f and χ f tend to be smaller for populations where a large fraction of the binaries possess either (i) mass ratio q close to unity, (ii) small magnitudes χ i , or (iii) small spin misalignments θ i∞ . Binaries with q → 1 result in more similar remnants because the spin magnitude S asymptotes to a constant [50]. Systems with either q → 0, χ i → 0, or θ i∞ → 0 are all well described by a single-spin approximation [33,51] where S is also approximately constant. Precession effects are less important for the distribution of v f because the contribution of the orbital-phase sampling is more dominant. In general, for the populations studied here we find δM f 0.004, δχ f 0.03, and δv f 1.5.

Relation to the spin morphologies
We now exploit our evolutionary procedure to highlight correlations between the remnant properties and the precessional dynamics at small separations. We make use of the concept of spin morphology [35,34,50,52,53]. The dynamics of BH binaries can be unambiguously classified into three different classes according to the qualitative behavior of the spin angle ∆Φ [Eq. (1)]. In particular, ∆Φ can either librate about 0 and never reach ±π (L0), circulate through the full range [−π, π] (C), or librate about ±π and never reach 0 (Lπ). All binaries belong to the circulating morphology at large separation, but GW radiation reaction will generically cause transitions towards the two librating classes before merger. In the following, we correlate the properties of the merger remnant to the spin morphology evaluated at ω orb→sur , i.e. right before initializing the NR surrogate.
The spin morphology before merger provides a useful mapping between the two asymptotic configurations of the evolution -the large-separation limit and the postmerger remnant. In Fig. 7 we illustrate this point showing a population of 10 5 sources with q = 0.95, χ 1 = 0.8, χ 2 = 0.7, and isotropic spin directions at large separations. Binaries in each each of the three morphologies tend to cluster in well defined regions of both the (θ 1∞ , θ 2∞ ) and the remnant parameter spaces. ‡ If ∆Φ librates about 0 (π), in-plane spin components are closer to alignment (anti-alignment) with each other and produce, on average, final BHs with masses and spins which are above (below) the hang-up anti-correlation (top-right panel of Fig. 7) In the left panel of Fig 8, we consider the same population and focus on the distributions of the recoil velocities. We find that ∼ 35% (∼ 100%) of the binaries with v f ≥ 1000 km/s (3000 km/s) belong to the ∆Φ ∼ π morphology before merger. This can be also be seen in bottom panels of Fig. 7, where red points preferentially populate the regions with large values of v f . The role of the spin morphology can be traced back to the very first "superkick" NR simulations [54,55] where BHs were initialized with |∆Φ| = π and θ i = π/2. In this setup, frame dragging of the the two BHs can act constructively, boosting the recoil [30]. Binaries with |∆Φ| = π and partial alignment θ i π/2 result in even larger recoils because more energy is emitted [56]. Conversely, binaries which librate about ∆Φ = 0 before merger present lower kicks of O(100) km/s and are responsible for the bulk of the population. Circulating sources present a rather flat distribution and allow for most values of v f but the largest ones. In the right panel of Fig 8, we repeat the same exercise but assume a mass ratio distribution p(q) ∝ q 6.7 [49] with q ∈ [0.25, 1] and uniform spin magnitudes Compared to the previous case, BHs have, on average, lower spins which implies weaker recoils [57] and a larger fraction of circulating binaries [35]. The general trends we highlighted above remain valid: the largest (smallest) kicks are largely due to binaries which were librating about ∆Φ = π (∆Φ = 0) before merger.

Transition thresholds
Finally, we investigate the robustness of our approach with respect to the two transition thresholds r pre→orb and ω orb→sur . Figure 9 shows how median and 90% interval of the three remnant quantities vary for a population of 1000 binaries with q uniformly We also show five random sources from the same population (gray curves in Fig. 9). The median levels of the remnant mass and spin do not depend strongly on the threshold values, though greater cyclic variation is present in the median levels of the remnant kick velocity. Over the full population this variation is not significant compared to the reported dispersion.

Conclusions
We presented an evolutionary procedure aimed at capturing the entire history of coalescing BH binaries. We combined multiple layers of PN integrations [34,35] and surrogate models trained on NR simulations [38]. Our scheme is robust with respect to the choice of the transition thresholds between the various techniques we employ (Sec. 3.5).  We map asymptotic binary configurations at r/M → ∞ to the properties of the BH merger remnant: final mass M f , final spin magnitude χ f , and recoil velocity v f . Crucially, this mapping depends on the sampling of two quantities: the precessional phase S pre→orb and the orbital phase ϕ orb→sur . Physically, this is because the binary completes infinitely many orbits and precession cycles before merger. Orbital and precessional phases impose a fundamental limit on the accuracy of models connecting the properties of astrophysical BH binary progenitors to the parameters of their merger remnants.
We first concentrated on individual sources and exploited our procedure to assess the impact of these uncertainties on the estimates of the remnant properties. Any given configuration at r/M → ∞ is mapped into a discrete region in the remnant parameter space (see Fig. 2). The relative errors that can be inferred from our results are δM f ∼ 0.1% for the remnant masses and δχ f ∼ 1% for the remnant spins. The orbital phase has a major impact on the kick velocity, resulting in δv f ∼ O(1).
We then disentagled the effects of the three different contributions to our uncertainties: the sampling of S pre→orb , ϕ orb→sur , and the NR surrogate errors. The orbital phase plays a dominant role and traces suggestive loop structures in the remnant parameter space (Fig. 3). The precessional phase determines the location and the precise shape of these loops. The remnant uncertainties are found to correlate with the projected effective spin ξ of the binary, with smaller (larger) values of |ξ| presenting large (smaller) dispersion. Uncertainties are more prominent for sources far from the equal-mass and single-spin limits.
Finally, we employed our evolutionary technique to study how the remnant properties correlate with the precessional dynamics at small separations. Merger remnants carry robust information about the precessional spin morphology in the late inspiral. Remnants originating from binaries with different morphological phases at small separations tend to cluster into different regions of the parameter space (Fig. 7). Final masses and final spins are strongly anti-correlated: binaries with spins close to (anti-)alignment produce remnants with lower (larger) masses and larger (lower) spins. BH binaries in the ∆Φ∼ 0 (∆Φ∼ π) librating morphology result in merger remnants that are systematically above (below) the M f -χ f correlation. Kicks are also strongly affected, with the vast majority of large recoils originating from binaries in the ∆Φ∼ π morphology.