Precision tracking of massive black hole spin evolution with LISA

The Laser Interferometer Space Antenna (LISA) will play a vital role in constraining the origin and evolution of massive black holes throughout the Universe. In this study we use a waveform model (IMRPhenomXPHM) that includes both precession and higher multipoles, and full Bayesian inference to explore the accuracy to which LISA can constrain the binary parameters. We demonstrate that LISA will be able to track the evolution of the spins -- magnitude and orientation -- to percent accuracy, providing crucial information on the dynamics and evolution of massive black hole binaries and the galactic environment in which the merger takes place. Such accurate spin-tracking further allows LISA to measure the recoil velocity of the remnant black hole to better than $100\,\mathrm{km}\,\mathrm{s}^{-1}$ (90\% credibility) and its direction to a few degrees, which provides additional important astrophysical information on the post-merger association. Using a systematic suite of binaries, we showcase that the component masses will be measurable at the sub-percent level, the sky area can be constrained to within $\Delta \Omega_{90} \approx 0.01 \, \rm{deg}^2$, and the binary redshift to less than $0.01$.


I. INTRODUCTION
The Laser Interferometer Space Antenna (LISA) will allow us to observe the mHz gravitational-wave (GW) spectrum.One of the key science goals of the mission is to trace the origin, evolution, and merger history of massive black hole (MBHs) throughout the Universe [1][2][3], where M ∼ 10 4 − 10 9 M ⊙ .These sources will be detected with large signal-to-noise ratios (SNR), allowing us to precisely measure the mass and spin angular momenta of the constituent black holes.
The population properties of the massive black hole binaries (MBHBs) that will be detected by LISA depends on a wide range of possible evolutionary pathways.Observing and characterising this population will provide vital insight into black hole seed formation, galaxy assembly, and the formation and evolution of structure throughout the Universe [1,[3][4][5][6].Whilst there are numerous open issues, there are several mechanisms that are of key importance but are still dominated by large theoretical and observational uncertainties.For example, we have a comparatively poor understanding of the mechanisms that generate the initial black hole seeds at high redshift (z ≳ 15) and the mechanisms that drive formation of the first protogalaxies.This, in turn, leads to a large uncertainty in the anticipated detection rates and the underlying astrophysical population properties [1,2,7].
MBHs are known to reside in the centres of galaxies, e.g.[8][9][10][11][12][13][14][15][16][17], and there is evidence that some galaxies contain a MBH, even at the early cosmological times probed at high redshift [18][19][20].Coupled with the hierarchical assembly of galaxies through repeated mergers [21][22][23][24], it is anticipated that LISA will observe numerous MBHBs across a wide range of cosmic epochs.However, there is huge uncertainty on our understanding of the astrophysical mechanisms that produced the initial seeds from which they grew.Did they originate from the collapse of heavy population III stars forming in low-metallicity environments (light seeds) [25][26][27] or from the collapse of large proto-galactic disks (heavy seeds) [28][29][30][31][32][33]?How did these black holes (BHs) subsequently grow?Accretion is an inevitable process and plays a key role in the growth of MBHs [4,34,35].However, there are other mechanisms that can drive the growth of MBHs, such as repeated hierarchical mergers [5,[36][37][38].Observations of MBHBs by LISA will provide us with a unique opportunity to tackle these questions and to constrain the evolutionary pathways responsible for the MBHs observed today.
For a binary of MBHs to merge, efficient mechanisms that dissipate the angular momentum of the system are required.Consequently, our current understanding of massive BH evolution can be broadly divided into distinct phases that are characterized by the astrophysical mechanisms driving angular momentum loss in the binary.The initial phase is primarily driven by dynamical friction caused by the gravitational interaction between the BHs and the surrounding gas and stellar material.Within this phase, the BHs are driven from a separation of approximately ∼ O(10 −2 Mpc) to ∼ O(1pc) [39].Below these separations, the dominant mechanism is hardening [39][40][41].Hardening refers to the interaction of fastorbiting stellar objects and the MBHB, resulting in energy being removed from the system.If the galactic host is sufficiently gas rich, a high-density circumbinary disk can form, leading to additional hardening of the binary through viscous dissipation [42][43][44].At scales below the hardening separation, e.g.[45], GW emission becomes the dominant mechanism driving the inspiral of the binary [46].However, depending on the physical properties of the BHs, and the environment of their galactic hosts, the astrophysical mechanisms at play can be widely dif-ferent.
A key tracer of the astrophysical formation mechanism will be the spin magnitude and orientation of the progenitor MBHs [1,2,24,[47][48][49][50][51], though these can be extremely challenging to measure [52].The orientation of the spin, in particular, is considered a key observable to discriminate between astrophysical formation pathways.When the BH spins are not aligned with the orbital angular momentum, they induce relativistic spin precession, a phenomena that drives the time evolution of the spins and the orbital plane, leaving characteristic imprints in the emitted GW signal, see Sec.III.
Of particular interest is the role of the BH spins during the phase of disc-driven migration.Notably, the spin orientations of the BHs by the time they reach the GW driven regime will be sensitive to whether the galactic environment is gas-rich or gas-poor.
If the environment is gas-rich, the inspiral of the binary can lead to the formation of a cavity in the circumbinary disc with accretion onto the BHs leading to the formation of minidiscs around each BH [53][54][55][56][57].If the viscous minidiscs are misaligned with respect to the BH spin, then Lense-Thirring precession (see [58]) will cause the inner disc to warp and dynamically relax, eventually reaching a coplanar configuration.This phenomenon is known as the Bardeen-Petterson effect [48,[59][60][61][62][63].However, the efficiency of Lense-Thirring precession is limited to the Bardeen-Petterson radius, beyond which the outer disc may maintain its original misalignment.A consequence, first noted in [64], is that the outer disc will drive the secular alignment of the BH spin with the total orbital angular momentum of the disc itself.This process occurs on a longer timescale than the warping of the minidisc but on a faster timescale than the growth of the BH, i.e. we can treat the mass and spin of the BH as being fixed [65,66].As such, BH spins that are preferentially aligned with the orbital angular momentum are thought to be a clear tracer for the Bardeen-Petterson effect which, in turn, would suggest a gas-rich galactic host [67][68][69].Conversely, MBH mergers in gas-poor environments are not expected to undergo any such alignment and we would expect the spin orientations to be generically oriented.If the spins are preferentially aligned, we are de-facto reducing the amount of precession in the binary relative to the generically oriented configurations.Spin precession is therefore a key observable in unraveling astrophysical formation pathways and will provide valuable insights into the environment of the galactic host.
A caveat to the discussion above can occur due to the torque on the disc induced by a companion BH.These torques can lead to the break-up of the disc and potentially result in BH spins that are misaligned, even in gasrich environments [69][70][71][72].A consequence is that this could introduce degeneracies between the environment and BH spin orientations [69].
Massive binaries, such as those discussed above, will have large signal-to-noise ratios in LISA, even at high redshift.For these binaries, we anticipate precision con-straints on the masses and spins of the constituent BHs.Early pioneering work, see [73], demonstrated the importance of spin-precession in breaking parameter degeneracies and in improving the precision to which the binary parameters could be constrained.This work was subsequently expanded in a series of papers, e.g.[74][75][76][77], that explored the impact of higher-order PN corrections, especially higher-order spin-orbit and spin-spin terms.Two particularly important conclusions are that the modulations in the GW signal induced by spin-precession allow one to accurately measure the BH spins and that precession also breaks the degeneracy between the binary orientation and its position on the sky, which could have important implications for identifying potential EM counterparts, e.g.[78][79][80].
LISA is sensitive to the entire sky, with the sensitivity depending on the source location and polarization, e.g.[81].The detector orbits the Sun such that its barycentere varies on a timescale of approximately one year.For long duration signals, the detector's translational motion induces a Doppler shift in the observed GW signal, which is sensitive to the angular position of the source [81].LISA also precesses around the normal to the ecliptic, leading to amplitude and phase modulations in the observed signal due to time-dependent variations in the antenna response pattern.The modulations of h + and h × depend on both the sky-location and orientation of the binary [81].Moreover, information on the angular position is encoded in the relative amplitudes and phases of the polarizations.Incorporating higher multipoles into the waveform model used to analyse the data will break parameter degeneracies, notably the degeneracy between the distance and the inclination.By accurately modelling higher multipoles, a more precise determination of the angular position and distance, e.g.see [82,83] for the context of MBHB in LISA.
When the spin angular momenta S i are misaligned with the orbital angular momentum L, precession occurs in the orbital plane and the spin vectors [84,85].This results in characteristic amplitude and phase modulations in the gravitational-wave signal, leading to several important observations.For example, in simple precession, the relative orientations of L and S i precess about the total angular momentum J [84], causing time-dependent changes in the binary's inclination ι, the polarization angle ψ, and the BH spins.The precession cycle's timescale is longer than the orbital timescale but typically shorter than the LISA timescale, depending on the masses and spins of the black holes.As a consequence, precessioninduced modulations enrich the gravitational-wave signal with valuable information about the binary and allow us to further break parameter degeneracies [73][74][75][86][87][88].The importance of precession in this context was first highlighted in the pioneering work of [73] and subsequently refined in [74][75][76][77]89].
We find broad agreement in the precision to which binary parameters are measured compared to the inspiralonly analyses, e.g.[1,[73][74][75].We demonstrate the ability of LISA to accurately constrain the spin orientation hours before merger and can provide tight constraints on the the kick velocity and orientation of the remnant BH.

A. Time-Delay Interferometry
Here, we broadly follow [82,98] to produce the LISA time-delay interferometry (TDI) response.In this framework, we apply a delay d(t) to a signal h(t) followed by a modulation F (t) which, in the notation of [98], can be written as The Fourier transform of the signal can then be written in terms of a transfer function T Following [98], the transfer function can be decomposed into a set of signal-dependent timescales that correspond to different physical effects.The first timescale of interest is which can be shown to reduce to the radiation-reaction timescale when the stationary phase approximation (SPA) holds [82] where ω is the orbital phase of the binary.This naturally leads to a time-frequency correspondence of where t c denotes the coalescence time of the binary.The effect of higher-order phase corrections to the transfer function can be written in terms of the time-frequency map as [98,99] T where G denotes a frequency-dependent kernel [98].This expansion is related to the shifted uniform asymptotics expansion (SUA) introduced in [99,100].In this study, we focus on the leading order approximation, i.e. p = 1, and defer a detailed analysis of the impact of higher-order corrections to the transfer function to future work.
The second key timescales arises from amplitude corrections beyond the leading order [98] which leads to a transfer function of the form [98] T Here we incorporate amplitude corrections up to second order, i.e. p = 2.
In practice, the TDI response is constructed in terms of a frequency-dependent transfer function acting on a set of waveform modes The modes are generated using the precessing, highermultipole inspiral-merger-ringdown model IMRPhe-nomXPHM [90][91][92].The term T A,E,T (f, t ℓm (f )) denotes the TDI transfer function and {A, E, T } the 1G TDI channels, which are related to the original 1G Michelson TDI variables {X, Y, Z} by [101,102] One of the key advantages of using the {A, E, T } variables is that they form a set of noise-orthogonal variables [101].As such, we assume that the noise in each channel is uncorrelated resulting in a diagonalized noise matrix that simplifies the resulting analysis [102].Note that the original variables introduced in [101] were constructed from the Sagnac variables {α, β, γ} and are therefore slightly different from the definition adopted here and in recent work, e.g.see the discussion in [93,103].Finally, we note that in order to generate the {X, Y, Z} TDI variables, we adopt the rigid adiabatic approximation (RAA) [98,104], in which where L = 2.5 × 10 9 m is the mean LISA armlength.

B. Bayesian Inference
A central aim of Bayesian inference is to reconstruct the posterior distribution p(θ|d) for the source parameters θ given data d.From Bayes' theorem we have where L denotes the likelihood of the data given the parameters θ, π the prior distributions for θ, and Z the evidence.We perform Bayesian inference using a coherent analysis of the full LISA TDI output, d = {d j ; j = A, E, T }, using the Balrog code.In particular, we employ nested sampling [105,106] as implemented by Dynesty [107,108].The likelihood of the data d given a set of binary parameters θ is [109] ln where h j denotes the TDI output for channel j produced by the GW signal.The noise weighted inner product is defined by where ã(f ) denotes the Fourier transform of the time series a(t) and S j (f ) denotes the noise power spectral density (PSD) of the j-th TDI channel.We use the noise spectral densities provided by the ESA Science Requirements Document [110,111] and include the unresolved galactic confusion noise according to the analytical fit of [112].As in [83], we adopt a low-frequency cut-off of f low = 10 −4 Hz, corresponding to the pessimistic scenario in which LISA only retains sensitivity over the frequency window specified in the mission design requirements [110,111].Any sensitivity at lower frequencies will help improve the constraints on measured parameters, so our results can be interpreted as a lower bound on the LISA science performance.Finally, we adopt the zero-noise approximation such that ñ(f ) = 0, noting that the noise still enters the likelihood in Eq. ( 17) through the PSD.We also assume that the PSD is constant over the observation duration, though in reality we would expect the PSD to be non-stationary and to vary slowly in time, e.g.[113].

III. MEASURING PRECESSION A. Black Hole Spins and Precession
Spin-induced precession is a unique and exciting phenomena that arises in the general-relativistic two-body problem [84,85,[114][115][116].If the black hole spins are misaligned with the orbital angular momentum, L, spinorbit and spin-spin couplings induce precession.The spin-orbit couplings drive a time-dependent evolution of the orbital plane, i.e. the direction of L, whereas the relativistic spin-spin couplings, which first enter at 2PN, induce a torque on the spin-angular momenta that results in the nutation of L, e.g.[84,85].The time evolution of L can be decomposed into terms corresponding to precession and a term driven by radiation reaction Here Ω L denotes the spin precession frequency in terms of a contribution arising from each of the BH spins where the contribution to the precession frequency from each spin can be written at next-to-leading order in terms of the effective aligned spin [116][117][118] as [117,119] where q = m 2 /m 1 ≤ 1 is the mass ratio of the binary.The dynamical effects that govern quasi-circular precessing binaries evolve on three distinct timescales.The orbital timescale is governed by t orb ∝ (r/M ) 3/2 , the precession timescale t prec ∝ (r/M ) 5/2 [84,85], and the radiation reaction timescale by t rad ∝ (r/M ) 4 .At large separations, r ≫ M , the timescales obey a hierarchy In the context of LISA, we expect MBH of total mass ∼ 10 4 −10 7 M ⊙ to be observable for O(days) to O(years), meaning that the binaries can undergo a comparatively large number of precession cycles compared to current ground-based observations of stellar mass black holes [120][121][122].Therefore, we can anticipate that MBHB observations with LISA will provide valuable opportunities for precision measurements of black hole spins and a comprehensive understanding of their evolution over multiple precession cycles.

B. Constraining Binary Parameters and Early Warning
In this study, we focus on a systematic series of investigations centred on two MBHBs that correspond to a low-spin (LS) and a high-spin (HS) configuration respectively.For each configuration, we consider three different inclination angles for the binary θ JN and three different spin tilt angles cos θ i = ⃗ χ i • L. The three spin tilts chosen correspond to a series of binaries in which the primary spin is approximately i) aligned with L, ii) perpendicular to L, and iii) anti-aligned with L. See Tables I and II in App.A for further details.Throughout this paper, we adopt HS2 as a canonical binary configuration, corresponding to a high-spin binary with χ 1 = 0.8, χ 2 = 0.6, ι = 42 degrees, θ 1 = 21 degrees and θ 2 = 31.4degrees.Throughout this work, we keep the redshift of the binaries the same, z = 3, as well as the detector-frame (i.e.redshifted) masses, corresponding to 6.4 × 10 6 M ⊙ and 1.6 Even at a redshift of z ≈ 3, the set of binaries considered here have large signal-to-noise-ratio ρ ≈ O(10 3 ) resulting in excellent constraints on the BH masses and spins, irrespective of the orientation of the binary θ JN or the orientation of the primary BH spin θ 1 .As in other recent studies of MBHBs, e.g.[82,83,124], we find that the individual masses are resolved to ∆m i ≈ O(0.1%), as shown in Fig. 1.
For both the high-and low-spin binaries, we broadly find that the spin magnitude of the primary is measured to an error on the order of ∆χ 1 ≈ O(0.1%), the secondary spin magnitude to ∆χ 2 ≈ O(1%), the primary tilt to within ∆θ 1 ≈ O(1%), and the secondary tilt to ∆θ 2 ≈ O(5%), as shown in Fig. 2 for the high-spin series and Fig. 7 in the Appendix for the low-spin series.A full list of recovered parameters and their 90% credible intervals can be found in Tab III.As we are able to resolve the spin magnitudes and orientations to exquisite precision, there is no need to rely on effective spin parameterizations, such as those proposed in [119,125].
Finally, we provide an assessment of the precision to which LISA can constrain the sky localization of the binary system when using an IMR waveform model that incorporates spin precession.For all binaries, the measurement accuracy for the ecliptic longitude is approximately ∆l ≈ O(0.05%), while the ecliptic latitude is constrained to ∆b ≈ O(1%).Furthermore, the inclination angle is determined to a precision of ∆θ JN ≈ O(0.1%), and the redshift to approximately ≈ O(0.1%).For the luminosity distance, we find that the distance to the binaries can be measured to with ∆d L ≈ O(0.01Gpc) which, assuming a Planck18 cosmology [123], translates to ∆z ≈ O(0.01).
These measurement accuracies correspond to a 90% credible sky area of ∆Ω 90 ≈ O(0.01deg 2 ) and a 90% confidence sky localization volume of ∆V 90 ≈ 7 × 10 −4 Gpc 3 .In contrast to [83], which used an aligned-spin model, our findings indicate an improvement in the angular resolution of LISA due to the incorporation of spin precession, in agreement with previous studies [73][74][75].The complete sky localization posteriors are shown in Fig. 3.The accuracy to which we can constrain the sky location and distance plays a critical role in determining the multimessenger prospects, e.g see the recent detailed study in [126] which outlines the number and properties of EMbright signals that will be observable in both LISA and a subset of EM facilities.
Finally, we note that the sky localization we report above is based on the analysis of the complete data from when it first enters the LISA band at f low ≈ 10 −4 Hz through to merger.In reality, the binary could be detectable before merger, depending on how the SNR accumulates.If a binary can be detected in sufficiently low-latency, then a preliminary skymap could be produced using partial data, e.g.adapting the methodology of [127], followed by a skymap generated using the complete data.As shown in Fig. 4, the SNR rapidly accumulates and passes a detection threshold, arbitrarily taken to be ρ det ≈ 12, on the order of 150 hours be- FIG. 2. Constraints on the spin magnitude and orientations for the high-spin series.The top row corresponds to the primary (heavier) BH and the bottom row the secondary (lighter) BH.The left-most panels depict a configuration where the primary spin is approximately aligned with L, the middle panels a highly-precessing configuration, and the right-most panels where χ1 is approximately anti-aligned with L. For each configuration, we consider three orientations that approximately correspond to: face-on (blue), edge-on (orange), and face-off (green) binaries.fore merger.This would provide ample time to request a protected LISA science data period to ensure the merger is recorded based on the current lead time estimate of ≈ 2 days [128].Note, however, that we do not make any attempt to tackle challenges pertaining to the the low-latency detection and localization of these binaries on the ≈ O(hrs) timescale.

C. Constraining Spin Evolution
In gas-rich environments, the binary can undergo a gas driven migration phase until the viscous timescales become smaller than the gravitational radiation-reaction timescale.At later stages of the inspiral, the binary evolution will be driven by GW emission.For the binaries considered here, the initial separation r i = r(f low ) is within the GW-dominated phase and we can neglect diskmigration effects on the evolution.In addition, we note that if one assumes that the total angular momentum of the binary is preferentially aligned with the circumbinary disk, then the spin tilt angles θ i are equivalent to those inherited from the gas-driven migration, see recent work in [69,[129][130][131][132].
We therefore evolve the binary from an initial separation r i to a final separation of r f = 5M using an orbit- FIG. 4. Growth of the signal-to-noise ratio (SNR) as a function of frequency or, equivalently, the time to merger (in hours).The horizontal dashed line denotes ρ det = 12, which we use as a proxy for the SNR at which a binary can be detected.For the canonical binary considered here, the system is visible up to ≈ 150 hours before merger.
averaged post-Newtonian evolution [133,134].The spin tilt angles of each BH are constrained to 1 degree or better which, when combined with the stringent mass constraints at the 0.01% level, allows for a precise determination of the evolution of the individual BH spins.
LISA will have the capability to accurately ascertain the magnitudes and orientations of the spins at a comparable level of precision, down to approximately 1 hour before the merger, as depicted in Fig. 5.As we shall see in the next section, detailed knowledge of the spin orientations at merger is vital if we are to accurately infer the recoil velocity and orientation.Finally, in the bottom panel of Fig. 5 we show the contributions of each spin Ω i to the precession frequency Ω L , as defined in Eqns.( 21) - (24).Precession has been measured in binary pulsars with typical precession frequencies being on the order of 10 −10 Hz, e.g.[135][136][137][138][139]. For the MBHs considered here, the precession frequencies fall in the mHz regime.

D. Constraining Remnant Properties
As is well known, GWs carry energy and angular momentum away from the system, driving the inspiral of the binary.The orbit shrinks until the two BHs eventually plunge and form a single, perturbed BH.The emission of GWs transfers linear momentum away from the binary, causing a displacement to its center of mass in the opposite direction [140][141][142][143].However, the loss of linear momentum through the merger is asymmetric and imparts a net linear velocity, i.e. a recoil or a kick, onto the remnant BH [144][145][146][147][148][149][150].The magnitude and orientation of the recoil has important astrophysical implications, with the first evidence for large recoil velocities being recently reported [150].This mechanism can displace or even eject BHs from the galactic center [151], potentially influencing the evolution of the MBH and its host galaxy [152,153].Recoils approaching 10 3 km s −1 could have a velocity that exceeds the escape velocity of the galactic host [151,154], suppressing the fraction of galaxies that host a central MBH [155,156], which has implications for the estimated event rates in LISA [3,157,158].
If a MBH is surrounded by a circumbinary disk, the perturbations to the gas induced by a recoiling BH can potentially give rise to an EM counterpart [159][160][161][162][163][164][165][166].The detectability of the EM counterpart is expected to depend significantly on the kick velocity of the recoiling BH and, potentially, to the orientation relative to the circumbinary disk.For the signal to be observable, a large recoil velocity on the order of O(10 3 km s −1 ) would be necessary.Kick velocities closer to O(10 2 km s −1 ) are unlikely to produce an observable EM counterpart, based on current theoretical models [165].The EM emission spectra will be highly reliant on the specific manner in which the disk undergoes shocks [166].Joint GW-EM inference will allow us to build a comprehensive understanding of the astrophysical processes at play in MBHB mergers.LISA will be able to accurately measure the masses and spins of the BHs as it approaches the merger, providing vital information on the geometry of the binary.
In order to understand the accuracy with which we can infer the remnant properties, we evolve the spins from f low through to the merger, as above, and use the NR calibrated fits of [167] to estimate the mass M f of the remnant BH, the magnitude of the kick velocity |v f |, and the orientation of the kick velocity θ v f .For the canonical binary we have focused on here, we find a moderate kick velocity of |v f | ≈ 575 km s −1 with an orientation of θ v f ≈ 165 deg, as shown in the top right panel of Fig. 6.The top left panel of Fig. 6 highlight the precision to which LISA will be able to resolve the binary geometry in terms of the orientations of Ĵ, χ1 , χ2 , and vf .As discussed in [167], these quantities are defined in an instantaneous frame in which the total orbital angular momentum Lorb points along ẑ at a reference time of t = −100M .Whilst a kick velocity of v f ≈ 575 km s −1 is likely to be too small to result in an EM-bright signal [166], the detection and accurate characterisation of the BH recoil through GW observations alone will allow us to place important constraints on the underlying astrophysical processes.See also the discussion in [66], which explores differential misalignment and the implications for the kick velocity.As discussed in [66], kick velocities on the order of ≈ 1500km s −1 can be produced in highly spinning binaries, with larger kick velocities only being anticipated for comparable mass BHs due to the dynamical alignment of the BH spins.The bottom panel Hz down to a separation of r = 5M , where we terminate the PN integration.On the top axis, we show the approximate time to merger in days, noting that we are considering a MBHB with a total redshifted mass of Mz = 8 × 10 6 M⊙ or, equivalently, a source-frame mass of Msrc = 2 × 10 6 M⊙.As noted in [83], a binary of this mass is expected to emit GWs in the LISA band for τ ≈ 6 days.Bottom panel: Posterior distributions on the evolution of the spin precession frequencies Ωi in Hz as defined in Eq. (21).Towards the end of the inspiral, the primary BH has a precession frequency O(20 mHz) and the secondary BH O(0.5 mHz).
of Fig. 6 shows the posterior distributions on the kick velocity as a function of the tilt of the primary spin at t = −100M for all configurations in the HS series.For the binaries in which the primary spin is essentially orthogonal to L, the approaches a maximum on the order of |v f | ≈ 1200 km s −1 .

IV. CONCLUSIONS AND FUTURE WORK
In this work, we have revisited the Bayesian parameter estimation of MBHBs in LISA using a state-of-the-art inspiral-merger-ringdown waveform model that incorporates both precession and higher multipoles [90][91][92].A FIG. 6.The top panels refer to the canonical "HS2" binary configuration.Top left: Posterior distributions for normalized total angular momentum and the dimensionless spin vectors of each BH.All quantities are defined in the orbital angular momentum frame L = ẑ at t = −100M [167].Top right: Posterior distributions for the magnitude of the kick velocity v f , the polar angle of the kick orientation, the remnant source-frame BH mass M f , and the dimensionless BH spin χ f .The bottom panels show the effect of different geometrical spin configurations for the "high spin binary" on the posterior distributions for the magnitude of the kick velocity v f .Here θ1 is the polar angle of the primary spin defined at a time t = −100M , as used to estimate the recoil velocity.
key focus was to investigate LISA's capability to accurately constrain the masses and spins of MBHs using the complete GW signal.Given the precision to which LISA can resolve the individual spin magnitudes and orientations, we highlighted how LISA will be able to accurately track the evolution of the BH spins through to the merger.An important consequence is that LISA will be able to accurately predict the magnitude and orientation of the recoil velocity, a key observable in exploring potential EM-bright counterparts to MBH mergers.In addition to the intrinsic parameters, we also demonstrated that, irrespective of the binary geometry, LISA will be able to constrain the sky location and orientation of the binary to sub-percent level, event out to a redshift of z ≈ 3, at least for the subset of binaries considered here.A detailed study exploring a larger population of MBHBs is underway to explore the broader parameter space and the concomitant implications for joint GW-EM analyses.
The detection rate of MBHBs is anticipated to be sensitive to the astrophysical distribution of spin magnitudes and orientations, especially for heavy binaries with M > 10 7 M ⊙ .Detecting and inferring the properties of these binaries will play a critical role in understanding the astrophysical origin and evolution of MBHs throughout the Universe.This also has important implications in exploring the complex interplay between the MBH population observed by LISA and the stochastic GW background observed by PTAs, for which significant evidence has now been observed [168][169][170][171][172][173][174].In a recent paper, [175] demonstrated the implications of observing this PTA background on the anticipated detection rate of MBHBs in LISA.Although [175] focused on a systematic suite of non-spinning binaries, the framework can be expanded to incorporate MBHBs with generically oriented spins, as explored here.
Detailed modeling of spin-precession plays a crucial role in resolving parameter degeneracies and enhancing the precision with which we can measure binary parameters [73,74,[87][88][89].It also plays a vital role in exploring MBH formation channels and constraining the astrophysics that drives their mergers.However, an important issue that remains unaddressed here is that of waveform systematics.The stringent accuracy requirements for future waveform models, necessary to reach this level of precision, pose immense challenges for source modelling.Recent developments, such as those in [176][177][178][179][180][181][182], have shown progress in enhancing the accuracy of the current generation of waveform models.Yet, achieving the accuracy and data analysis demands of LISA still remains a significant endeavor.In future work, we will explore the impact of waveform systematics on our ability to constrain the astrophysics behind MBH mergers.Moreover, we note that the GW signal can be influenced by environmental effects, such as dynamical friction or torquing from angular momentum exchange between accretion discs and the binary [183][184][185][186]. Detailed phenomenological modeling of environmental effects opens a promising avenue to discern the astrophysical processes at play in galactic nuclei but could introduce systematic biases into our analysis [187].
The results presented here are a first step towards integrating waveform models with precession and higher multipoles into a global fit pipeline.Notably, [97] has demonstrated the application of aligned-spin phenomenological waveform models to address the MBHB aspect of the global fit.In our case, we have employed IMRPhe-nomXPHM to perform Bayesian inference on individual MBHBs.Finally, we note that IMRPhenomXPHM can be effectively combined with the heterodyned likelihood [188], even when including precession and higher multipoles [189].In forthcoming work, we will build upon the findings presented here to demonstrate how MBHBs can effectively constrain the underlying astrophysical population [190] and facilitate tests of General Relativity [191].I. Summary table for the low spin (LS) binaries, characterized by χ1 = 0.2 and χ2 = 0.4.We consider three configurations, corresponding to θ1 = {21 • , 78 • , 169 • } degrees respectively.For each configuration, we consider three inclination angles ι = {10 • , 42 • , 80 • } to gauge the impact of the binary orientation on our ability to resolve the binary parameters.Here 10 degrees corresponds to a near face-on configuration and 80 degrees to an approximately edge-on configuration.The masses correspond to the source-frame values and a dashed line denotes that the value is equivalent to the row above.

Appendix B: Summary of Bayesian Inference
Here we list the exact configurations used for the lowspin (LS) and high-spin (HS) binaries.For each set of binaries we fix the sky-location, distance, component masses, spin magnitudes, and spin tilt of the secondary BH.We then systematically vary the spin tilt of the primary BH from an approximately aligned configuration through to an approximately anti-aligned configuration.Spins that are misaligned with the orbital angular momentum give rise to relativistic spin couplings that drive the precession of the orbital plane and the spins themselves [84,85].By tilting the primary spin we are effectively changing the amount of precession in the binary.In addition to varying the primary spin magnitude, we also explore three binary orientations for each tilt configuration to explore the impact of observing the binary from a near face-on configuration through to a near edge-on configuration.Due to the fact that the binaries undergo multiple precession cycles after entering the observable LISA band, we find that the orientation has little impact on our ability to constrain the BH spins.

SUPPLEMENTARY PLOTS
Figure 7 is the low-spin counterpart of Fig. 5, demonstrating that even for low spins, LISA can accurately resolve the magnitudes and orientations of BH spins to a comparable level of accuracy.

CORNER PLOTS
The complete joint posteriors for our canonical binary configuration, HS2, are shown in Fig. 8.

5 FIG. 1 .
FIG. 1. Corner plot showing the joint posteriors on the source frame component masses.The mass of the primary BH is measured to within O(10 3 M⊙) and the mass of the secondary to O(10 2 M⊙).

t
[hours to merger]

FIG. 5 .
FIG. 5. Tracking precession of binary HS2.Top panel:Posterior distributions on the spin evolution from the initial binary separation defined by the starting frequency of f low = 10 −4 Hz down to a separation of r = 5M , where we terminate the PN integration.On the top axis, we show the approximate time to merger in days, noting that we are considering a MBHB with a total redshifted mass of Mz = 8 × 10 6 M⊙ or, equivalently, a source-frame mass of Msrc = 2 × 10 6 M⊙.As noted in[83], a binary of this mass is expected to emit GWs in the LISA band for τ ≈ 6 days.Bottom panel: Posterior distributions on the evolution of the spin precession frequencies Ωi in Hz as defined in Eq.(21).Towards the end of the inspiral, the primary BH has a precession frequency O(20 mHz) and the secondary BH O(0.5 mHz).
As in Tab.I, but now for the high spin (HS) binaries, where χ1 = 0.8 and χ2 = 0.6.

TABLE III .
Summary of the 90% credible intervals on the recovered parameters for all injections considered in this paper.Here HS denotes the high-spin series and LS the low-spin series.The component masses mi correspond to the source-frame values, χi the dimensionless spin magnitudes, θi the spin-tilt angles, and ϕ12 the difference in the azimuthal spin angles between the two BH spins.JN = 10• θ JN = 42 • θ JN = 80 • FIG.7.As per Fig.2but for the low-spin systematic series. θ