Parameter estimation of stellar-mass black hole binaries with LISA
Alexandre Toubiana, Sylvain Marsat, Stanislav Babak, John Baker, Tito Dal Canton
PParameter estimation of stellar-mass black hole binaries with LISA
Alexandre Toubiana,
1, 2
Sylvain Marsat, Stanislav Babak,
1, 3
John Baker, and Tito Dal Canton APC, AstroParticule et Cosmologie, Universit´e de Paris,CNRS, Astroparticule et Cosmologie, F-75013 Paris, France Institut d’Astrophysique de Paris, CNRS & Sorbonne Universit´es, UMR 7095, 98 bis bd Arago, 75014 Paris, France Moscow Institute of Physics and Technology, Dolgoprudny, Moscow region, Russia Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD 20771, USA Universit´e Paris-Saclay, CNRS / IN2P3, IJCLab, 91405 Orsay, France
Stellar-mass black hole binaries (SBHBs), like those currently being detected with the ground-basedgravitational-wave (GW) observatories LIGO and Virgo, are also an anticipated GW source in the LISA band.LISA will observe them during the early inspiral stage of evolution; some of them will chirp through the LISAband and reappear some time later in the band of 3 rd generation ground-based GW detectors. SBHBs couldserve as laboratories for testing the theory of General Relativity and inferring the astrophysical properties of theunderlying population. In this study, we assess LISA’s ability to infer the parameters of those systems, a crucialfirst step in understanding and interpreting the observation of those binaries and their use in fundamental physicsand astrophysics. We simulate LISA observations for several fiducial sources, setting the noise realization tozero, and perform a full Bayesian analysis. We demonstrate and explain degeneracies in the parameters of somesystems. We show that the redshifted chirp mass and the sky location are always very well determined, withtypical errors below 10 − (fractional) and 0 . . The luminosity distance to the source is typically measuredwithin 40 − O (1%). The error onthe time to coalescence improves from O (1 day) to O (30 s) as we observe the systems closer to their merger.We introduce an augmented Fisher-matrix analysis which gives reliable predictions for the intrinsic parame-ters compared to the full Bayesian analysis. Finally, we show that combining the use of the long-wavelengthapproximation for the LISA instrumental response together with the introduction of a degradation function athigh frequencies yields reliable results for the posterior distribution when used self-consistently, but not in theanalysis of real LISA data. I. INTRODUCTION
Debuting with the first detection of a stellar-mass black holebinary (SBHB) in 2015 [1], the LIGO / Virgo collaboration hasissued the first catalog of the gravitational wave (GW) sourcesidentified during the first and second observing runs (O1 andO2) [2–7] and three noteworthy detections from the third ob-serving run (O3) [8–10]. These observations of GW in the10–1000 Hz band, which include eleven SBHBs mergers andtwo binary neutron star (BNS) mergers, have inaugurated theera of GW astronomy and opened a new window to the uni-verse, allowing to infer for the first time the properties of thepopulation of compact binaries [11, 12] and providing newtests of general relativity (GR) [13–15].The Laser Interferometer Space Antenna (LISA) [16],scheduled for launch in 2034, will observe GWs in a di ff er-ent frequency band (the mHz band) and, therefore, comple-ment ground-based detectors. The strongest anticipated GWsources in the LISA data will be massive black hole binaries(MBHBs), with total mass in the range 10 − M (cid:12) [17],and galactic white dwarf binaries (GBs). The latter are sonumerous that they will form a stochastic foreground signaldominating over instrumental noise in addition to a smallernumber ( ∼ ) of individually resolvable binaries [18]. SB-HBs with a total mass as large as those observed by LIGOand Virgo could also be detected by LISA during their earlyinspiral phase, long before entering the frequency band ofground-based detectors and merging [19]. SBHBs in the mHzband can be at very di ff erent stages of their evolution, rang-ing from almost monochromatic sources to chirping sources which leave the LISA band during the mission [19]. We focuson resolvable SBHBs, one of the best candidates for multi-band observations [19, 20].Although SBHBs are not LISA’s main target, the scientificpotential of multiband observations with LISA and groundbased detectors is considerable. These could be used toprobe low-frequency modifications due to deviations from GR[19, 21–26] or to environmental e ff ects [27–31], to facilitateelectromagnetic follow up observations [27], or simply to im-prove parameter estimation (PE) over what is possible withground-based interferometers alone [19, 24]. More precisemeasurements would improve the testing of competing as-trophysical formation models. Several scenarios have beensuggested for the formation of SBHBs, such as stellar evo-lution of field binaries versus dynamical formation channels[32, 33]. Moreover, the possibility that these black holes areof primordial origin [34] cannot be completely discarded. Thevarious possible formation channels typically predict di ff er-ent distributions for the parameters of SBHBs, especially ec-centricity and spin orientation / magnitude [35–37] providingdiscriminating power in astrophysical model selection. A re-cent study has suggested that already a few “special” events,binaries with high primary mass and / or spin, would have ahuge discriminating power [38]. The specific impact of obser-vations of SBHBs with LISA on astrophysical inference wasconsidered in [39–43].LISA will observe MBHBs somewhere between a few daysand few months before their merger, which corresponds tothe end of inspiral when the binary is evolving rapidly [17].On the contrary, GBs are slowly evolving, almost monochro-matic sources and they will remain in band during the whole a r X i v : . [ g r- q c ] J u l LISA mission [44]. Resolvable SBHB signals will fall in be-tween these two behaviours: they are long lived sources butare not monochromatic and some SBHBs can chirp and leavethe LISA band. In addition, all resolvable SBHBs are clus-tered at the high end of LISA’s sensitive band. Thus, SBHBswill produce very peculiar signals of great diversity. In thiswork we do not address the question of how to detect thosesources, although it has been argued to be challenging [45].Instead, we assume that we have at our disposal an e ffi cientdetection method, and concentrate on inferring the parame-ters of the detected signals. We also consider only one signalat the time, while in reality the SBHB signals will be super-posed in the LISA data stream. The PE study presented herewill also be valuable in building search tools as we discuss atthe end of the paper.Most previous studies of PE for SBHBs with LISA reliedon the Fisher approach, and used simple approximations toLISA’s output response to GW signals. While a quick ande ffi cient method for forecasting studies, the Fisher approachmight not be suited for the systems with low signal to noiseratio (SNR) and non-Gaussian parameter distributions [46].SBHBs signals are long lived and emit at wavelengths com-parable to LISA’s size, this implies that we need to properlytake into account the complete LISA response. As a result, theuse of the long-wavelength approximation (also called low-frequency approximation) [47] might not hold and could se-riously bias the PE [48, 49]. In this work we consider thefull LISA response as described in [50, 51] and perform a fullBayesian analysis in 0-noise of all the systems we consider.We also provide a comparison to Fisher-matrix-based PE andbriefly comment on the impact of using the long-wavelengthapproximation.Despite the on-going e ff ort to infer the astrophysical forma-tion channel of the SBHBs observed by LIGO / Virgo, a hugeuncertainty still remains. The situation will improve as wedetect more signals. We expect that 3 rd generation detectors(Einstein Telescope [52–54], Cosmic Explorer [55]) will beoperational in parallel with LISA, with SNR figures reachinghundreds or thousands, thus significantly improving the PEover current observations. Given the current uncertainty in thepopulation of SBHBs, we cannot reliably specify the proper-ties of most detectable sources [19, 39, 45, 56–59]. In ourstudy we focus on a fiducial system consistent with the pop-ulation of currently observed SBHBs instead of working witha randomized catalog of sources. From there, we perform asystematic scan of the parameter space by varying few param-eters at a time, investigating their qualitative impacts on thePE. We consider quasi-circular binaries consisting of spinningblack holes, with the spins aligned or anti-aligned with the or-bital angular momentum, merging no later than 20 years fromthe beginning of observations. We start with a GW150914-like system [1] and explore the parameter space by varyingat most three parameters at a time. For each of these sys-tems we infer the posterior distribution assessing correlationsbetween parameters and the accuracy in measuring each pa-rameter across the parameter space.The paper is organised as follows. In section II we describehow we generate GW signals and our tools to perform PE. Then we give details on how we choose all the systems onwhich we perform PE in section III . Detailed description andanalysis of PE results are given in section IV. There we alsoprovide a comparison to a slightly modified version of Fishermatrix analysis and assess the validity of long-wavelength ap-proximation for the PE. Finally, in section V we discuss thescientific opportunities o ff ered by LISA observations of SB-HBs in the light of our results. II. ANALYSIS METHODA. Bayesian framework
Data measured by LISA ( d ) will consist in a superpositionof GW signals ( s ) and a noise realisation ( n ): d = s + n .The instantaneous amplitude of a GW signal is much lowerthan noise making their detection very challenging. We usematched filtering as a main technique for detection, the mainidea is to search for a specific pattern (GW template) in thedata [60]. It is done by correlating the data with a set ofGW templates in frequency domain (˜ h ( f , θ ) which are func-tions of source parameters θ ). This correlation is given by thematched-filter overlap( d | h ) = R e (cid:32)(cid:90) + ∞ ˜ d ( f ) ˜ h ∗ ( f ) S n ( f ) d f (cid:33) , (2.1)where S n ( f ) is the power spectral density (PSD) of the de-tector noise, assumed to be stationary. In this work, we usethe LISA “Proposal” noise level as given in [16]. We do notdiscuss detectability of SBHBs in this paper, we assume thatall sources discussed here could be detected and we concen-trate on the parameter extraction / estimation. We should notethat the detection itself could be a challenge, at least for thetraditional method of template banks [45].We work in a Bayesian framework for the parameter es-timation, treating the set of parameters of the source, θ , asrandom variables. The Bayes theorem tells us that given theobserved data d , the posterior distribution p ( θ | d ) is given by: p ( θ | d ) = p ( d | θ ) p ( θ ) p ( d ) . (2.2)In the right hand side of this equation p ( d | θ ) is the likelihood, p ( θ ) is the prior distribution and p ( d ) is the evidence. As thenoise and GW signal models will be fixed in this study, the ev-idence can be seen as a normalisation constant which does notneed an explicit calculation. Assuming noise to be stationaryand Gaussian, the likelihood is given by: L = p ( d | θ ) = exp (cid:34) −
12 ( d − h ( θ ) | d − h ( θ )) (cid:35) . (2.3)In order to speed up the computation, we set the noise real-ization to zero as n =
0, so that d = s . The addition of noiseto the GW signal is not expected to drastically a ff ect the PE,leading at most to a displacement of the centroid of the pos-terior distribution within the confidence intervals (CI) (withthe probability defined by CI). Thus, the analysis of the pos-terior distribution itself should remain representative in pres-ence of noise (still assuming Gaussianity). We consider onlyone source at a time and we neglect all possible systematicerrors due to signal mismodelling: s = h ( θ ) with θ the pa-rameters of the GW source. Under these simplifications, thelog-likelihood is given by (up to a normalization constant in L ): log L = −
12 ( h ( θ ) − h ( θ ) | h ( θ ) − h ( θ )) . (2.4)Since LISA will only observe the inspiral phase of thesebinaries, we expect that the dominant 22-mode is su ffi cientand we neglect the contribution of all other subdominant har-monics. We use the model called PhenomD [61, 62] to gen-erate ˜ h , ± and compute the LISA response to generate thetime delay interferometry (TDI) observables A, E and T (seee.g. [63]) as described in [50, 51]. The three TDI observablesconstitute independent data sets, therefore the log-likelihoodis actually a sum of three terms like Eq. (2.4), one per TDIobservable.Similarly to the treatment of galactic binaries, weparametrise the sources by their initial frequency and phaseat the start of observations, instead of the time to coalescenceand phase at coalescence (more suitable in LIGO / Virgo dataanalysis or for MBHBs with LISA). We define the initial timeas the moment LISA starts observing the system. A systemis characterised by (i) 5 intrinsic parameters: the masses ( m and m ), the GW frequency at which LISA starts observingthe system ( f ) and spins ( χ and χ ); and (ii) 6 extrinsic pa-rameters: the position in the sky defined in the solar systembarycenter frame (SSB) ( λ and β ), the polarisation angle ( ψ ),the azimuthal angle of the observer in the source frame ( ϕ ),the inclination of the orbital angular momentum with respectto the line of sight ( ι ) and the luminosity distance to the source( D L ).We introduce a set of sampling parameters for whichwe believe the posterior distribution should be a sim-ple function, i.e. close to either a uniform or Gaus-sian distribution. The choice of sampling parame-ters is made based on the post-netwonian (PN) expres-sions for GW [64–66]. We use the set of parameters θ = ( M c , η, f , χ + , χ − , λ, sin( β ) , ψ, ϕ, cos ι, log ( D L )), where M c = m / m / ( m + m ) / is the chirp mass, η = q (1 + q ) is the sym-metric mass ratio with q = m m > χ + = m χ + m χ m + m is the e ff ective spin (often denoted χ e ff in theliterature) and χ − = m χ − m χ m + m is an antisymmetric spin com-bination.For the simulated data we assume a sampling rate of 1 Hz sothe Nyquist frequency is f Ny = . f upto f max = min( f Ny , f T obs ) where f T obs is the frequency reachedby the system after the observation time T obs . We considertwo mission durations: T obs = T obs = ffi cient way to explore the parameter space. We do this bymeans of a Markov Chain Monte Carlo (MCMC) algorithm[67]. More specifically we designed a Metropolis-HastingsMCMC (MHMCMC) [68] for this purpose that we presentnext. A less costly alternative would be to use the parameterestimation based on the Fisher information matrix. We willshow how one can modify the Fisher matrix to make it a robustPE tool in the following subsections. We exploit the metric in-terpretation of Fisher matrix in the MHMCMC, thus, we startby reviewing some basics on the Fisher matrix approach anddelegate comparison with Bayesian PE to subsection IV C. B. Fisher matrix
In the Fisher matrix approach, the likelihood is approxi-mated by a multivariate Gaussian distribution [46]: p ( d | θ ) ∝ e − F ij ( θ ) ∆ θ i ∆ θ j , (2.5)where ∆ θ = θ − θ and F is the Fisher matrix given by: F i j ( θ ) = ( ∂ i h | ∂ j h ) (cid:12)(cid:12)(cid:12) θ , (2.6)where the partial derivative ∂ i denotes the derivative with re-spect to θ i . Here again we actually have a sum of three terms,one for each TDI observable. We assume this to be implicitin the following. The inverse of F is the Gaussian covariancematrix of the parameters, it gives an estimate of the error oneach parameter. This approach has been extensively used inthe studies of LISA’s scientific capability thanks to its simplic-ity. However for systems with low SNR as the ones we con-sider, the Fisher approximation might not be valid [46] and weneed to perform a full Bayesian analysis.The Fisher matrix has an alternative interpretation: it canbe seen as a metric on the parameter space associated with thedistance defined by the inner product (2.1): || h ( θ + δθ ) − h ( θ ) || = ( h ( θ + δθ ) − h ( θ ) | h ( θ + δθ ) − h ( θ )) (cid:39) ( ∂ i h δθ i | ∂ j h δθ j ) (cid:12)(cid:12)(cid:12) θ = F i j ( θ ) δθ i δθ j . (2.7)We exploit this property in our MHMCMC sampler. C. Metropolis-Hastings MCMC
We sample the target distribution p ( θ | d ) by means of aMarkov chain, generated from a transition function P ( θ, θ (cid:48) )satisfying the detailed balance condition: p ( θ | d ) P ( θ, θ (cid:48) ) = p ( θ (cid:48) | d ) P ( θ (cid:48) , θ ) . (2.8)We build the transition function from a proposal function π such that P ( θ, θ (cid:48) ) = π ( θ, θ (cid:48) ) a ( θ, θ (cid:48) ) where a ( θ, θ (cid:48) ) is the accep-tance ratio defined as: a ( θ, θ (cid:48) ) = min(1 , p ( θ (cid:48) | d ) π ( θ (cid:48) ,θ ) p ( θ | d ) π ( θ,θ (cid:48) ) ) if π ( θ, θ (cid:48) ) (cid:44)
00 otherwise . (2.9)It is easy to verify that P satisfies the detailed balance condi-tion for any choice of π . In practice, a jump from a point θ to θ (cid:48) is proposed using the function π ( θ, θ (cid:48) ) and the new pointis accepted with probability a ( θ, θ (cid:48) ). If the point is not ac-cepted the chain remains at θ . In both cases, we update thecurrent state of the chain. By repeating this procedure, we ob-tain a sequence of samples of the target distribution. We seefrom the expression of the acceptance ratio that points withhigher posterior are more likely to be accepted, thus the chainwill tend to move towards regions of higher posterior explor-ing part of the parameter space compatible with the observeddata. In theory, the chain should converge independently ofthe chosen proposal function and of where the chain is startedbut it could take an infinite time. Since for PE we are in-terested in high posterior regions we start the chain from theinjection point, i.e. the maximum likelihood point. The max-imum likelihood point coincides with the maximum posteriorpoint if all priors are flat, but it is not true in general and themaximum posterior point can depend on the adopted prior dis-tribution. Even though the posterior does not depend on theproposal, the convergence, e ffi ciency and resolution of tailsof distribution very strongly depends on a particular choice,the best proposal should closely resemble the target (poste-rior) distribution. Thus, most of the work goes into buildingan e ffi cient proposal function. Note that for a symmetric pro-posal ( π ( θ, θ (cid:48) ) = π ( θ (cid:48) , θ )), the acceptance ratio is simply givenby the ratio of the posterior distributions. This specific case iscalled Metropolis MCMC [69] and is the one we consider.Our runs are done in two steps: we first run a short MCMCchain ( (cid:39) points) to explore the parameter space and thenuse the covariance matrix of the points obtained from thischain to build a multivariate Gaussian proposal that we usein a longer chain. During the first stage, called burn-in, weuse a block diagonal covariance matrix. We split the set of pa-rameters in 3 groups: intrinsic parameters ( M c , η , f , χ + , χ − ),angles except the inclination ( λ , sin( β ), ψ , ϕ ) and inclination-distance (cos ι , log ( D L )). Each block is computed invertingthe Fisher matrix of that group of parameters. This separa-tion was based on the intuition, well verified in practice, thatthe stronger correlations are within these groups of parame-ters and is intended to avoid numerical instabilities that mayarise when dealing with full Fisher matrices. Notice that bymaking this choice we do not discard possible correlations be-tween parameters of di ff erent groups, we are simply not takingthem into account when proposing points based on the Fishermatrix. If those correlations exist, they should appear in theresulting covariance matrix that we use to build a proposalfor the main chain. Failing to include existing correlationscould reduce the e ffi ciency of our sampler in its exploratory,or burn-in, phase; however the splitting can easily be adaptedif needed.We rotate the current state vector θ to the basis of the covari-ance matrix’s eigenvectors. In this basis the covariance matrixis diagonal, formed by the eigenvalues of the covariance ma-trix in the original basis. Because for some parameters thedistribution is very flat, the eigenvalues of the covariance ma-trix predicted by Fisher can be very large, reducing the ef-ficiency of the sampler. This is usually the case for poorly constrained but bounded parameters like cos ι and spins. Toavoid this issue we truncate the eigenvalues of the (cos ι, D L )matrices and define an e ff ective Fisher matrix accounting forthe finite extent of the prior on spins: F e ff = F + F p . We take F p χ + ,χ + = F p χ − ,χ − = σ with σ = .
5. We motivate this choicein Sec. IV C.In order to improve e ffi ciency of the sampler in the eventof complicated correlations between intrinsic parameters, weexploit the metric interpretation of Fisher matrix and occa-sionally recompute the covariance matrix for the first groupof parameters with a given probability. By doing so we mightviolate the balance equation, but it is done only during theburn-in stage (exploration of the parameter space) and thosepoints are discarded later from the analysis.We test the convergence of the chains by running multiplechains with di ff erent random number generator seeds, check-ing that they all give similar distributions and computing theGelman–Rubin diagnosis [70] for all the parameters. Potentialscale reduction factors below 1.2, as the ones we get, indicatethat the chains converged [70]. For each chain we accumulate10 − independent samples (by thinning the full chain bythe autocorrelation length) which takes 4 − III. SETUPSA. Systems
We start by considering a system with masses and spinscompatible with the first detected GW signal (GW150914)and label it
Fiducial system. Its parameters are given in Ta-ble I along with its SNR assuming LISA mission lifetimes T obs = T obs = m = (1 + z ) m s where z is thecosmological redshift and subscripts s denote parameters inthe source frame. We adopt the cosmology reported by thePlanck mission (2018) [72]. Note that T obs is the mission du-ration, not the time spent by the system in the LISA band andwe assume an ideal 100% duty cycle. The initial frequencyis derived from the time to coalescence from the beginning ofLISA observation that we fix to 8 years for the Fiducial sys-tem, thus for T obs = Fiducial system is observed fora fraction of its inspirals, while for T obs = f refers to the Fiducial system.We explore the parameter space of SBHBs by changing afew parameters of the
Fiducial system at a time. We list allthe systems we consider in the following subsections, speci-fying what are the changes with respect to the
Fiducial systemand the corresponding labels. For all systems we considertwo possible mission duration quoted above (unless some-thing else is specified).In Table II we show the considered systems and their re-spective SNR. Note that we chose to use the LISA proposalnoise level [16], which does not include a 50% margin in- m (M (cid:12) ) 40 m (M (cid:12) ) 30 m , s (M (cid:12) ) 36 . m , s (M (cid:12) ) 27 . t c (yrs) 8 f (mHz) 12 . χ . χ . λ (rad) 1 . β (rad) π/ ψ (rad) 1 . ϕ (rad) 0 . ι (rad) π/ D L (Mpc) 250 z . T obs (yrs) 4 10SNR 13 . . Fidu-cial . The masses and spins of this system are compatible withGW150914 [71]. The initial frequency is computed such that thesystem is merging in 8 years from the start of LISA observations.We consider two possible durations of the LISA mission: 4 and 10years (in the latter case, the signal stops after 8 years at coalescence).Subscripts s denote quantities in the source frame, bare quantities arein the detector frame. The sky location is given in the SSB frame. troduced to form the “Science requirements” SciRDv1 [73].The SNRs would thus be significantly more pessimistic withSciRDv1. From the point of view of the PE, using one or theother noise model amounts to a constant rescaling of the noisePSD S n , with the same e ff ect as rescaling the distance to thesource.
1. Intrinsic Parameters
Unless specified we take t c = t c . Changing t c (or equivalently f ) amounts in shifting the GW signal infrequency and also defines its frequency bandwidth (withinthe chosen observation time). • Time left to coalescence at the beginning of LISA ob-servations:
Earlier : t c =
20 years,
Later : t c = • Chirp mass keeping the mass ratio unchanged:
Heavy : M c = . M c , f , q = q f , D L =
445 Mpc
Light M c = M c , f . , q = q f , D L =
150 Mpc • Mass ratio, keeping the chirp mass unchanged: q3 : q = m m = M c = M c , f q8 : q = m m = M c = M c , f • Spins:
SpinUp : χ = . χ = . SpinDown : χ = − . χ = − . SpinOp12 : χ = . χ = − . SpinOp21 : χ = − . χ = . Heavy and
Light systems we scaled the distance sothat the SNR remains the same as for the
Fiducial system inthe case T obs = ff ects the SNR so we do not change the distance for thosesystems. Since the Earlier system merges in 2 years, increas-ing the observation time from 4 to 10 years has no impact.
2. Extrinsic parameters
Changes in extrinsic parameters do not a ff ect the time to co-alescence so all systems below have the same initial frequencyas the Fiducial system. • Position on the sky (in the SSB frame):
Polar : β = π − π , λ = λ f Equatorial : β = π , λ = λ f • Inclination:
Edge-on : ι = π − π , D L =
150 Mpc, T obs = • Distance:
Close : D L =
190 Mpc, T obs = Far : D L =
350 Mpc, T obs = Very Far : D L =
500 Mpc, T obs = Edge-on system tosustain a reasonably high SNR. For the same reason, we useonly T obs = ff ect of varying thenoise level and the duty cycle. For the Close system we onlyconsider the T obs = Far and
Very Far systems we only consider the T obs = B. Prior
Regarding the Bayesian analysis, we take our fiducial priorto be flat in m and m with m ≥ m , flat in spin magnitudebetween − T obs = T obs = Fiducial . . Earlier . . Later . / Heavy . . Light . . q3 . . q8 . . SpinUp . . SpinDown . . SpinOp12 . . SpinOp21 . . Polar . . Equatorial . . Edgeon / . Close . / Far / . Very Far / . ff erent systems are derived fromthe Fiducial system, varying a few parameters at once. We use the
Full response. since only 2 ϕ (for a 22-mode waveform) and 2 ψ intervene,we restrict to an interval of π . We obtain the prior probabil-ity density function (PDF) in terms of the sampling param-eters by computing the Jacobian of the transformation from( m , m , χ , χ , D L ) to ( M c , η, χ + , χ − , log ( D L )) which gives: p f ( θ ) = N M c η − / D L √ − η if 0 . ≤ η ≤ . , . (3.1)Just like the evidence in Eq. (2.2), N acts only as a normal-isation constant and thus it is of no importance for us. Thelower limit for η was set according to the maximum mass ra-tio up to which PhenomD is calibrated ( q =
16) [61, 62]. Therange of chirp mass, initial frequency and distance are ordersof magnitude larger than the posterior support so they do nota ff ect the posterior. We label this prior as Flatphys and use itby default unless we specify some other choice, for example,we will consider two additional priors: • Flatmag : uniform prior for the spins orientation andmagnitude • Flatsampl : flat prior in M c , η and log ( D L ).In the Flatmag case we start from a full 3D spin prior, uniformin [0,1] for the spins amplitude and uniform on the sphere for the spins orientation. We then consider only the spin pro-jections on the orbital momentum, thus ignoring the in-planespin components. The resulting prior is p ( χ i ) = − ln( | χ i | ).The Flatmag
PDF is: p ( θ ) = p f ( θ ) p ( χ ) p ( χ ) where p f ( θ ) isgiven in Eq. (3.1). This is the prior generally used by theLIGO / Virgo collaboration [2].The
Flatsampl
PDF is given by: p ( θ ) = N η if 0 . ≤ η ≤ . , . (3.2)This prior has no astrophysical motivation, we will use it tocompare Fisher-based-PE to our full Bayesian inference inSec. IV C.We find it instructive to illustrate how the non trivial priorslook like. As we will show later in section IV, the chirp masscan be constrained by Bayesian analysis to a fractional errorof 10 − , so we can impose a narrow constrain on the prior. Thechirp mass is non-trivially coupled to other parameters (as wewill show in great details in the following sections), and con-straining it to the narrow interval introduces non-linear slic-ing in other parameters. Note that the imposed interval (10 − in relative terms) is still much broader than the typical mea-surement error. In Fig. 1 we display the Flatphys , Flatmag and
Flatsampl prior distributions for η , χ + and χ − obtainedby restraining the chirp mass to the specified interval. Theremarkable features of our fiducial prior, the Flatphys prior,are the double peak at η = .
25 and η = χ + and χ − priors with almost zero support at ex-treme values. The Flatmag is singled out by the strong peakat χ + , − =
0. As we will discuss in section IV, these non-trivialshapes of the priors can strongly a ff ect the resulting posteriordistributions in some cases. C. LISA response
We briefly review how the LISA response is computed andrefer to [51] for a more extensive discussion. We recall thatLISA is composed of three spacecraft linked by lasers acrossarms of length L = . A , E and T are time-delayed linear combinations of the single-link ob-servables y slr that measure the laser frequency shift due toan incoming GW across link l between spacecrafts s and r .We consider only the dominant 22-mode of the waveform,and following [51] we exploit a mode symmetry (for non-precessing systems) between h and h , − to write the signalin terms of h only. The single-link observables can then bewritten with a transfer function:˜ y slr = T slr ˜ h . (3.3)Denoting the amplitude and phase of the mode 22 as ˜ h = A ( f ) e − i Ψ ( f ) , working at leading order in the separation oftimescales in the formalism of [50] the transfer functions are FIG. 1. Comparison between the
Flatphys (blue),
Flatmag (green)and
Flatsampl (orange) priors for η , χ + and χ − . given by (we set c = T slr = G slr ( f , t f ) (3.4) G slr ( f , t ) = i π f L (cid:2) π f L (1 − k · n l ) (cid:3) · exp (cid:104) i π f (cid:16) L + k · ( p Lr + p Ls ) (cid:17)(cid:105) · exp(2 i π f k · p ) n l · P · n l (3.5) t f = − π d Ψ d f , (3.6)where k is the unit gravitational wave propagation vector, n l ( t )is the link unit vector pointing from the spacecraft s to r , p ( t )is the position vector of the center of the LISA constellation inthe SSB frame, p Lr ( t ) is the position of spacecraft r measuredfrom the center of the LISA constellation and P is the polar-isation tensor defined in [51]. We dropped the t dependence inEq. (3.5) for more clarity. The global factor exp(2 i π f k · p ) isthe Doppler modulation in GW phase and the n l · P · n l termis the projection of the GW tensor on the interferometer axeswhich is associated with the antenna pattern function. Notethat both the Doppler modulation in phase and the antennapattern are time dependent due to LISA’s motion, moreover,they depend on the sky position of the source, so that the an-nual variation in the phase and amplitude allows us to localizethe GW source.All our results are obtained using the full LISA responsebut we also assess the impact of using the long-wavelengthapproximation, a simplified version of the LISA response[47]. In this approximation, LISA is somewhat similar to twoLIGO / Virgo-type detectors moving around the sun with theangle between arms being π/
4. It is obtained by taking the 2 π f L (cid:28) G slr ( f , t ) = i π f L i π f k · p ) n l · P · n l (3.7)The sinc function appearing in Eq. (3.5) leads to a dampingof the signal amplitude at high frequencies but in the long-wavelength approximation it is replaced by 1, leading to unre-alistically high SNRs. To compensate for this, inspired by thecomputation of the sky averaged sensitivity [74] we introducea degradation function that multiplies the GW amplitude: R ( f ) = + . π f L ) . (3.8)To explore the validity of this approximation for SBHBs,we will compare the PE for the Fiducial , Polar and
Equato-rial systems using the full response and the long-wavelengthapproximation labeled
Full and LW respectively. We will onlyuse the leading order in the separation of timescales in theframework of [50], keeping in mind that corrections couldbe needed in general, in particular for almost-monochromaticsignals. IV. PARAMETER ESTIMATION OF SBHBS
In order to test the performance of our MHMCMC sam-pler we compared it to our well tested parallel temperingMCMC code ptmcmc . The quality of agreement betweentwo distributions p and p can be quantified by computingthe Kullback-Leibler (KL) divergence [75]: D KL = (cid:88) θ p ( θ ) log (cid:32) p ( θ ) p ( θ ) (cid:33) . (4.1)The KL divergence is zero if two distributions are identical.We computed D KL for the marginalised distributions of eachparameter obtained with two samplers using the Flatsampl prior, and assuming 4 and 10 years of observation. Apart from ψ and ϕ , all divergences were below 0 . .
01 for 10 years of observation showinga very good agreement between samplers. For ψ and ϕ , lesswell determined in general, we get slightly higher values (upto (cid:39) Flatphys prior andthe
Full response. In our discussion, we use redshifted masses(rather than source frame) because they are directly inferredfrom the observed data. We give the full “corner plot” [76]for our fiducial system, comparing results for the two obser-vation times in Fig. 2, this plot shows pair-wise correlationbetween parameters and the fully marginalised posterior foreach parameter. The inset on the right top of the figure showsposterior distributions for ( m , m , χ , χ ). https: // github.com / JohnGBaker / ptmcmc FIG. 2. Inferred parameter distribution for the
Fiducial system, both in the T obs = T obs = It would be di ffi cult to represent the posterior distributionsfor all possible variations (deviations from the Fiducial ) dis-cussed above. Instead, we will summarize our results by un-derlining qualitative di ff erences whenever we observe themand show comparative corner plots only when necessary. Westart by discussing the structure of correlation between intrin-sic parameters, move to extrinsic parameters, then comparethe full Bayesian analysis with predictions from the Fishermatrix and, finally show the e ff ect on the PE of the LW ap-proximation to the response. A. Intrinsic parameters
One of the main features appearing in Fig. 2 is the strongcorrelation between intrinsic parameters, in particular the onebetween M c and η which is especially pronounced for 4 yearsof observation. The main reason for this degeneracy is thelimited evolution of the GW frequency: in 4 years of obser-vation the Fiducial system spans a very narrow range from f = . f = . f , χ + , χ − andall extrinsic parameters to the “true” values and investigate thecorrelation between M c and η for the Fiducial system. Keep-
FIG. 3. Analysis of the degeneracy between M c and η . The blue (orange) dots were obtained running a PE on the Fiducial system in the T obs = T obs = M c and η to vary. The injection point is indicated by the black dashed lines. The orange dottedand the green dashed curves are given by (4.5) using the full PhenomD phase and the 1.5 PN truncation of the phase respectively. The red solidline was obtained by minimising the phase di ff erence between the injected signal and templates over the whole frequency range spanned over4 years of observation.FIG. 4. Individual PN phase contributions ∆Φ n for the Fiducial sys-tem. The linestyle indicates the nature of the term, non-spinning(NS), spin-orbit (SO) or spin-spin (SS), while the color indicates thePN order. Note that these contributions are individually aligned at f , as explained in the text, and that interpreting the magnitude ofthese terms is not easy due to the alignment freedom. The verticalline shows f , and the greyed area shows the frequency range con-tributing less than 1 in SNR . ing these parameters fixed will collapse some of the degen-eracies seen in the full analysis, but this exercise will serve asan illustration of the di ff erences between a non-chirping andchirping system.First, it will be instructive to consider the magnitude of thedi ff erent post-Newtonian (PN) orders appearing in the phasing(see [64] for a review). We can write formally Ψ ( f ) = ηv (cid:88) i a i v i , (4.2) FIG. 5. Value of δ I Ψ along the curve in the ( M c , η ) plane that min-imises it for I = [ f , f y ears ] (blue) and I = [ f , f max , LIS A ] (orange).When LISA observes the system at low frequencies, the phase di ff er-ence can be kept small over an extended region far from the injection.When LISA observes the chirp of the system, the phase di ff erencebecomes very large immediately at the vicinity of the injection point,reducing the extent of the degeneracy between M c and η . where v = ( π M f ) / and where the a i are PN coe ffi cients (wescaled out the leading term, so that a =
1) that depend on themass ratio and on the spins and can be separated between non-spinning terms (NS), spin-orbit terms (SO) and spin-squareterms (SS). It was argued in [77] that most SBHBs would re-quire terms up to the 2PN order. In Fig. 4, we show the magni-tude of the know PN terms in the phasing for our
Fiducial sys-tem. In general, the magnitude of phase contributions is deli-cate to interpret because of the alignment freedom, as some ofthe phasing error can typically be absorbed in a time and phaseshift. In Fig. 4 we align contributions individually at f with azero phase and zero time according to (3.6). We see that, for T obs = FIG. 6. Samples of χ and χ obtained for di ff erent systems (defined in section III) in the T obs = − ≤ χ , ≤ χ PN = χ PN , lines. The samples follow the χ PN = χ PN , lines, showing thatthis is the specific combination of spins that can be measured. The orthogonal combination of spins is constrained only due to the boundariesof the physically allowed region. Due to the orientation of the χ PN = const lines, χ is better constrained than χ . High values of spins withsame (opposite) sign are the better (worse) constrained. the limited chirping in frequency, while more terms becomerelevant for T obs = f > , which wetake as a somewhat conventional limit to indicate that ignor-ing the signal beyond this point would not a ff ect the likeli-hood (2.4) and therefore the PE.Since for T obs = f to f , we can Taylor expand the phase around f : Ψ ( f ) (cid:39) Ψ ( f ) + d Ψ d f (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f ( f − f ) +
12 d Ψ d f (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f ( f − f ) . (4.3)We consider the inner product between the data d = A d ( f ) e − i Ψ d ( f ) and the template h = A h ( f ) e − i Ψ h ( f ) . From ourconvention, the initial phase at f is the same, Ψ d ( f ) =Ψ h ( f ). The initial time is zero at f , so the stationary phaseapproximation (3.6) gives: d Ψ d f | f =
0. The inner product be-comes:( d | h ) = (cid:90) f f d f A d ( f ) A h ( f ) e i ( Ψ d ( f ) − Ψ h ( f )) S n ( f ) d f (cid:39) A d ( f ) A h ( f ) 4Re (cid:34)(cid:90) f f d fS n ( f ) e i ( d2 Ψ d d f | f − d2 Ψ h d f | f ) ( f − f (cid:35) , (4.4) where we used the fact that the amplitude is a slowly varyingfunction of the frequency. The overlap is maximized whenwhen the template is in phase with the data, making the inte-grand non-oscillating. In our quadratic approximation to thedephasing, this defines a curve in the ( M c , η ) plane accordingto d Ψ d f (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f = d Ψ ( M c , , η )d f (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) f . (4.5)In Fig. 3 we display in blue (orange) dots points from thesampling in the ( M c , η ) plane in the T obs = T obs = T obs = ff er-ence between injection and template over the whole frequencyrange spanned by the injected signal. More specifically, defin-ing: δ I Ψ ( M c , η ) = max I | Ψ ( M c , , η )( f ) − Ψ ( M c , η )( f ) | , (4.6)1 FIG. 7. Distribution of η and χ + for the Later system ( t c = Fiducial system ( t c = T obs = T obs = Later system chirping, thedetermination of η and χ + is much better than for the Fiducial systemin the T obs = = . η = .
25, as an e ff ect of theprior. This is on contrast to the Fiducial system in the T obs = = .
1) which peaks at the injected value indicated byblack lines and squares. for each value of M c we find η such that δ I Ψ is minimised.Note that all parameters are kept fixed in the dephasing mea-sure we use here, in particular there is no optimization overa constant phase or time shift. The subscript I stands for thefrequency interval and we plot this curve for I = [ f , f y ears ]in Fig. 3. One can see that we almost perfectly reproduce theshape of the correlation between M c and η in the T obs = T obs = δ I Ψ for I = [ f , f y ears ]and for I = [ f , f LISAmax ] with f LISAmax = . ∼ . δ I Ψ to be quite small ( < . χ . Asthe bandwidth of the signal becomes broader, we cannot e ffi -ciently compensate for a change in the chirp mass by varying η which results in the significant reduction of the degeneracyand great improvement in measuring those two parameters asseen by the narrower region covered by the orange dots onFig. 3.We now come back to the full Bayesian analysis and con-sider the estimation of the black holes spins. Following [62, 78] we introduce: χ PN = (cid:32) χ + + q − q + χ − (cid:33) (4.7) = η (cid:32) (113 q + χ + ( 113 q + χ (cid:33) . (4.8)This term defines how the spins enter the GW phase at theleading (1.5 PN) order [64] and, therefore, should be deter-mined the best from observation of SBHBs. We found thisto be indeed well verified. As an illustration, we plot sam-ples obtained for the q3 , q8 , SpinUp , SpinDown , SpinOp12 and
SpinOp21 systems in the T obs = χ PN = χ PN , (fixing the mass ratio to its “true” (injec-tion) value) for all those systems in the ( χ , χ ) plane. In allthese cases χ PN is extremely well measured, within 10 − , butthe combination of spins orthogonal to χ PN is constrained onlyby the prior boundaries.For slowly evolving binaries, only terms up to 1.5 PN inthe GW phase are found to be relevant. At this order we ex-pect a strong correlation between χ PN and η : any change in χ PN can be e ffi ciently compensated by a change in η suchthat the 1.5 PN term ( − π + χ PN ) η − / is kept (al-most) constant. We have verified this by plotting the curve( − π + χ PN ) η − / = const on top of the samples obtainedfor the Fiducial system and reproducing the shape formedby the posterior samples. Thus, we obtain and explainedthe three-way correlation between chirp mass, mass ratio andspins for the mildly relativistic systems spanning a narrow fre-quency band during the observation time. The increase in theobservation time allows further chirping of the system, mak-ing the contribution of the 1PN and 1.5PN corrections in thephasing significant, thus breaking strong correlations betweenintrinsic parameters; however, the e ff ect of higher-order PNterms is weak, as consistent with [77] and 4, which leads toonly the spin combination χ PN to be measured. This studyalso suggests that χ PN , being the most relevant mass-weightedspin combination, should be used as sampling parameter. Thecomponent of χ PN along χ + is always much larger than theone along χ − (at least by a factor (cid:39) χ + isalso measured reasonably well and it is more relevant for theastrophysics of SBHBs. We will alternate between χ PN and χ + in our next discussions.In order to further quantify the dependence of PE on thefrequency bandwidth spanned by the signal during the obser-vation time, we consider the Earlier , Fiducial and
Later sys-tems which di ff er in the initial frequency chosen so that theSBHBs merge in 20, 8 and 2 years respectively. We com-pute the KL divergence between the marginalised posteriorand the marginalised prior for each intrinsic parameter, andreport our findings in Table III. The larger values of D KL in-dicate that knowledge has been gained from the GW observa-tions as compared to the prior. The results show a strong de-pendence on the observation time (therefore on the frequencybandwidth), especially for spins, for which the D KL variesby an order of magnitude. For the Earlier system we findthat only the chirp mass measurement is truly informative.2
Fiducial Earlier Later M c η χ + χ − χ PN M c η χ + χ − χ PN M c η χ + χ − χ PN Flatphys T obs = . . . . . . .
04 0 .
04 0 .
03 0 .
04 6 . . . . . T obs = . . . . . . . . . . / / / / / Flatmag T obs = . . .
07 0 .
04 0 . / / / / / / / / / / T obs = . . . . . / / / / / / / / / / Flatsampl T obs = . . . . . / / / / / / / / / / T obs = . . . . . / / / / / / / / / / TABLE III. Kullback-Leibler divergences between the marginalised posterior and prior distribution of the intrinsic parameters for di ff erentsystems and choices of prior. When observing the system at low frequencies, only M c shows a sensible deviation from the prior. The likelihoodis informative on η and χ + (and χ PN ) only for chirping systems. Di ff erent choices of prior give similar results. Note that the longer frequency evolution plays a bigger rolethan the SNR. For instance,
Later which leaves LISA after2 years with SNR = . Earlier with T obs = = .
2. We repeatedthis analysis using the
Flatmag and
Flatsampl priors for the
Fiducial system. For all choices of prior, the KL divergencesare similar, proving the η , χ + , χ − distributions are prior dom-inated when observing slowly evolving systems. Notice thatKL divergences for spins are slightly smaller when using the Flatmag prior, meaning that the posterior is even more dom-inated by the prior. This is because the
Flatmag prior peaksstrongly at χ + = χ − = D KL are always larger for χ PN than for the otherspin combinations, reflecting the fact that it is the best mea-sured spin combination. Still, for systems evolving througha narrow frequency interval, the χ PN distribution is also priordominated. The e ff ect of the prior is especially well seen forthe Fiducial system and T obs = η at 0.25 is what we expect due to prior (see section III B).The same peak is also observed for the Later system (predom-inantly due to low SNR) but η is much better constrained forthis system, the likelihood is informative enough to reduce thewidth of the distribution, but not large enough to supress theprior. Let us reiterate this important finding: for intrinsic pa-rameters beyond the chirp mass, the chirping (extend of thefrequency evolution) of the observed SBHB has stronger in-fluence on PE than the SNR or observation time per se.We note that, although the frequency is slowly evolving,the signal is far from monochromatic unlike most of galac-tic binaries (e.g. double withe dwarf binaries). As an ele-ment of comparison, using the quadrupole formula to computethe frequency derivative at f , for the Earlier system we find˙ f = . × − Hz which is 4 orders of magnitude higherthan the fastest evolving galactic binaries [18]. Thus, despitethe strong correlation between intrinsic parameters, the chirpmass is always well measured, with a relative error of order10 − for the Earlier system when observing for 4 years andbelow 10 − for the chirping systems. The tight constraint on M c leads to the banana-like shape correlation between m and m seen on the top right part of Fig. 2. As a result, we can de- termine individual masses (within 20 − Fiducial systemin Table IV. Whenever the marginalised distribution of a givenparameter is leaning against the upper (lower) boundary ofthe prior as for m ( m ) we define the 90% CI as the valuebetween the 0.1 and 1 quantile (0 and 0.9). Otherwise, in allother situations we define the 90% CI as the values betweenthe 0.05 and 0.95 quantiles. In all cases we report the medianas a point estimate.Systems with a higher mass ratio ( q3 and q8 , keeping thechirp mass the same as for Fiducial ) give an error on the chirpmass similar to the
Fiducial system, but the mass ratio is bet-ter determined. This is because, when keeping the chirp massfixed, the PN expansion of the GW phase features negativepowers of η , notably in the 1PN term. Moreover, what shouldmatter is the derivative of the phase with respect to η whichcontains only negative powers ( η − / , η − / ) which makes itmore sensitive to the small mass ratio as compared to theequal-mass systems. For an observation time of 4 years, theuncertainty on individual masses is still of order of 100%, butfor an observation time of 10 years, it reaches below 10% and1% for the q3 and q8 system respectively.We now discuss the e ff ect of priors on PE for high-spin sys-tems. Consider SpinUp system in the T obs = χ PN ) mass ratio η and the chirp mass. Inthe posterior we observe the interplay between the η and χ + priors which push samples towards η = .
25 and χ + = χ + (0.95). This,together with the correlation between parameters, leads to theresulting posterior distribution which has double peak in η andbroad distribution for χ + (the 2-D histogram is more informa-tive). The distribution (overall) is shifted away from the truevalues (well evident in the right panel of Fig. 8), though theyare still contained within 90% CI. In the case of T obs = FIG. 8. The left panel shows the inferred distribution on η and χ + for the SpinUp system. Because of a “competition” between the prior and thelikelihood the distributions of η and χ − peak away from the true value indicated by the black lines and and the square. The M c distribution, notshown, is marginally a ff ected. Because of the bias in η , the inferred distribution of masses is significantly biased. However with our definitions,the true value is within the 90% CI.FIG. 9. The left (right) panel shows the inferred distribution on η and χ + ( m , χ ) for the Fiducial system using the
Flatphys , Flatmag and
Flatsampl priors. Under the e ff ect of the prior, the posterior distribution can be significantly shifted away from the true value indicated byblack lines and squares. tems are badly constrained and closely resemble the priors.For chirping systems, the determination of spins can be un-derstood from Fig. 6. Because of the orientation of lines χ PN = const, χ is better constrained than χ . As the mass ra-tio increases the slope of these lines changes, accentuating thisdi ff erence. Spins of same (opposite) sign, are better (worse)determined as their magnitude increase because of the nar-rowing (broadening) of the allowed region. For the Fiducial system, the error on χ is quite large but we can infer that thespin is positive with 0 (and negative values) being outside the90% CI given in Table IV. The e ff ective spin ( χ + ) is measured within 0 . M c , s = M c / (1 + z ), we get: ∆ M c , s M c , s = ∆ z + z + ∆ M c M c . (4.9)As we dicuss in Sec. IV B 2, D L is typically measuredwithin 40 −
60% which implies a measurement of the redshift z within ∼ −
60% (at the low redshifts we are considering,4 T obs = T obs = M c / M c , + × − − × − + × − − × − M c , s / M c . s , . + . − . . + . − . q . + . − . . + . − . m / m , . + . − . . + . − . m / m , . + . − . . + . − . m , s / m , s , . + . − . . + . − . m , s / m , s , . + . − . . + . − . χ + . + . − . . + . − . χ − . + . − . . + . − . χ PN . + . − . . + . − . χ . + . − . . + . − . χ . + . − . . + . − . ∆ t c (s) 10 ∆Ω (deg ) 0 .
18 0 . D L / D L , . + . − . . + . − . z . + . − . . + . − . TABLE IV. 90% CI on the parameters of the
Fiducial system whoseparameters are given in Table I using the
Flatphys prior. For massesand distance, we give the relative errors. The redshifted chirp massis extremely well determined for both mission durations but indi-vidual masses can be measured only if the mission is long enoughand we can observe the system chirping. The measurement of thesource frame chirp mass is worse, being dominated by the error onthe distance measurement and therefore the redshift in (4.9). Theerror on individual masses is dominated by their intrinsic degener-acy. For chirping systems, we can also measure χ PN which translatesinto a good constraint on the e ff ective spin χ + . The error on individ-ual spins remains large for the chirping system, but we can start toconstrain the spin of the primary black hole (in our example, exclud-ing negative values). As a consequence of the overall improvementin the determination of the intrinsic parameters, the inference of thetime to coalescence improves drastically. The sky location (given byEq. 4.11) is very well determined for both mission durations, withinthe field of view of next generation electromagnetic instruments likeAthena and SKA [79, 80]. D L and z are linearly related). Thus, the second term on theright hand side of Eq. (4.9) is clearly subdominant and the er-ror on the source frame chirp mass is dominated by the erroron redshift, as the result we get ∆ M c , s M c , s (cid:39) . z + z . (4.10)This error is typically of the order of a few percent for systemsdetectable by LISA (up to z ∼ O (10 − )), which is better thancurrent LIGO / Virgo measurements [2]. This estimate is ingood agreement with the results presented in Table IV. Theerror for individual masses in the source frame are dominated by the error on the masses (like the second term on the lefthand side of Eq. (4.9)) due to poorly constrained mass ratio.The initial frequency is always extremely well determined,with relative errors below 10 − . Its determination improvesfor chirping systems due to reduction of correlation with otherintrinsic parameters. The frequency of the system at the begin-ning of the LISA observations is coincidental, as it is directlylinked to the time that is left for the system to coalesce t c . Weapply the stationary phase approximation (3.6) to the full GWphase to infer t c . This transformation involves all intrinsic pa-rameters, so the error on t c is typically smaller for chirpingsystems. We find an error of the order of 1 day for systems farfrom merger, while for more strongly chirping systems t c canbe determined to within 30 s.We find that increasing or decreasing the total mass of thesystem (while preserving the SNR) as in the Heavy and
Light systems has little consequence on the estimation of intrinsicparameters. The error on spins and symmetric mass ratio arethe same as in the
Fiducial case. The relative error on chirpmass and initial frequency is slightly smaller for lighter sys-tems (factor (cid:39) . Heavy and
Light systems) be-cause of the larger number of cycles. However, we do not finda simple scaling with the chirp mass of the system for a fixedlevel of SNR. In particular, we do not find the error on thechirp mass to scale with M / c as computed in [81, 82]. Thiswas to be expected since as discussed in this section the er-ror on intrinsic parameters depends crucially on the frequencyinterval through which we observe the binary.Finally, the choice of prior only marginally a ff ects the pos-terior distribution for chirping systems. On the other hand, itcan have a significant impact for non-chirping systems as canbe seen in Fig. 9. For example, the Flatmag prior completelydominates the posterior distribution of spins as the KL diver-gences suggested and shown in Fig. 9. Because of the notedcorrelations, the prior on spins propagates into determinationof mass ratio and individual masses.
B. Extrinsic parameters
1. Sky location
The sky location of the source is very well determined and,except for systems close to the equator, its posterior distribu-tion is very similar to a Gaussian unimodal distribution. Wedefine the solid angle as in [47]: ∆Ω = π (cid:112) ( Σ λ,λ )( Σ sin( β ) , sin( β ) ) − ( Σ λ, sin( β ) ) , (4.11)where Σ is the covariance matrix. This defines a 63% confi-dence region around the true location. The error for the Fidu-cial system, reported in Table IV, is below 0 . whichis within the field of view of most planned electromagneticinstruments such as Athena and SKA. [79, 80]. With theexception of the Equatorial system, the sky position is con-strained with a similar precision for all systems considered inthis work. The good localization comes from the complicatedmodulations imprinted on the signal by the orbital motion of5
FIG. 10. We plot f | ∂ i ˜ h | / S n (normalised to its maximum value) as afunction of the time to coalescence. This quantity is the integrand ofthe diagonal Fisher element (2.6) and indicates where in frequencythe information on the parameter θ i comes from. We show this quan-tity for λ (the behaviour for sin( β ) is very similar) and M c (the be-haviour for the other intrinisc parameters is very similar). We indi-cate the initial and end frequencies of the Later (red),
Fiducial (red)and
Early (black) systems for T obs = Polar , Fiducial and
Equatorial systems,with T obs = ff ects near the pole we do not compare the angles in the SSB frame( λ, β ) but transformed angles ( µ, γ ) defined by placing the injectionpoint at the equator in each case (note that the scale of the two axisis not the same). The injection corresponds to µ = γ = µ is equally well recovered inthe three cases. For the Equatorial system we find a tail extending tothe position β → − β . ∆Ω (deg ) T obs = T obs = Fiducial .
18 0 . Earlier .
70 0 . Later . / Polar .
14 0 . Equatorial .
74 0 . Polar sky position β (cid:39) ± π/
2, but much worsefor the
Equatorial sky position β (cid:39)
0. The sky localization is betterfor the
Later system than for the
Earlier system (despite a lower SNRin the T obs = LISA, according to (3.5). To understand how the sky localiza-tion evolves as a function of the frequency band we observea system, in Fig. 10 we plot f | ∂ λ ˜ h | / S n (normalised with re-spect to its maximum value) as a function of the frequency.The quantity f | ∂ i ˜ h | / S n is the integrand entering the compu-tation of the diagonal elements of the Fisher matrix (2.6) andindicates (for each parameter) the most informative frequencyrange. Using a logarithmic scale for frequencies, the factor f ensures that we can visualize the contributions to the integralas the area under the curve (up to a normalisation factor). Wealso indicate the corresponding values of the time to coales-cence t c on the upper x-axis. We indicate the initial (dashedline) and end (solid line) frequencies of the Later (red),
Fidu-cial (red) and
Early (black) systems for T obs = β ) is similar to λ . For comparison we showthe same quantity for M c , the behaviour for other intrinsic pa-rameters is similar. As discussed in the previous section, mostof the information on intrinsic parameters comes from highfrequencies. On the other hand, there is more information onthe sky location at low frequencies, where a given range offrequencies corresponds to more orbital cycles of the LISAconstellation. However, this is to be balanced with the nar-rower frequency range spanned by systems evolving at lowerfrequencies, for a fixed observation time. For this reason, the Later system gives a better localization than the
Earlier sys-tem even in the T obs = ).We can distinguish two main e ff ects in (3.5) informing usabout the sky localization: the time-dependency (through t f ,see (3.6)) of the response reflects the orbital cycles of LISA,and the Doppler modulation exp(2 i π f k · p ) of the phase.The Doppler modulation shows this time dependency, but alsoscales with f , so this term is larger for chirping signals reach-ing high frequencies. We find a better sky localization forlighter systems: ∆Ω Light < ∆Ω Fiducial < ∆Ω Heavy (ranging from0 . . in the T obs = .
02 to 0 . in the T obs = t c and the SNR (by adjusting the dis-6tance) for those systems. The GW signal from the lighter andheavier systems is displaced at higher and lower frequencies,since the evolution rate of the inspiral depends primarily onthe chirp mass. Namely, f = . , . , . Heavy , Fiducial , and
Light , respectively. Since we keepthe SNR fixed in this comparison, this means that the lightersystem has a stronger sky-dependent Doppler modulation ofthe phase, helping with the localization.When comparing the
Polar , Fiducial and
Equatorial sys-tems, a direct comparison of the sky localization could bequite misleading because the metric on a sphere depends onthe latitude, with a singularity at the pole. To evade this issuewe define a system of coordinates on the sphere ( µ, γ ) suchthat the injection point is always on the equator. The trans-formation from the ecliptic coordinates to this frame is sourcedependent. The spherical coordinates at the equator are lo-cally Cartesian and simplify the comparison of the results. Weshow the results of the sky localization in Fig. 11 for the
Po-lar , Equatorial and
Fiducial systems in the ( µ, γ ) frame andfor T obs = µ (azimuthal angle)similarly well but the determination of γ worsens as β → Equatorial system we find a tail extend-ing all the way to a secondary sky position corresponding to β → − β . This behaviour is due to the dominant Doppler phasein the frequency response which goes as cos β : although theamplitude of the e ff ect itself is maximized, its variation with β is minimal as cos β is flat for β =
0. For T obs = f -dependent in theresponse (3.5) are larger, and the total SNR itself is larger. Thesolid angle for the Equatorial system is larger as compared toother systems (as reported in Table V) but remains well belowthe current sky localization with ground-based observatories[2].
2. Other extrinsic parameters
We find strong correlations between inclination and dis-tance, and between the polarisation and the initial phase.These degeneracies are commonly seen in the analysis ofLIGO / Virgo sources when using only the dominant 2 , ± , ± h + ( f ) = ˜ A ( f ) 1 + cos ( ι )2 e i ϕ e − i Ψ ( f ) , (4.12a)˜ h × ( f ) = i ˜ A ( f ) cos ι e i ϕ e − i Ψ ( f ) , (4.12b)where ˜ h ( f ) = A ( f ) exp( − i Ψ ( f )) is the frequency domainamplitude and phase decomposition of the mode h , with˜ A ≡ √ / π A ( f ) absorbing conventional factors. We referto [51] for notations; in particular we exploit the symmetrybetween h and h , − for non-precessing systems to write thewaveform in terms of h only. Going to the SSB frame we FIG. 12. Comparison of the inferred distribution on ψ , φ , ι and D L for the Far and
Edge-on systems, both have similar SNR. We nor-malised the distance to the injection value. Black lines and squaresindicate the true values common for both systems and coloured linesand squares the value of ι for each system. For the Edge-on systemthe degeneracies between φ and ψ and between ι and D L are bro-ken giving a better estimation of each of these parameters. However,close to edge-on systems will usually have much lower SNR. Indeed,in order to keep a comparable SNR the distance to the Edge-on is lessthan half the distance to the
Far system. rotate by the polarisations angle ψ :˜ h SSB + = ˜ h + cos(2 ψ ) − ˜ h × sin(2 ψ ) , (4.13a)˜ h SSB × = ˜ h + sin(2 ψ ) + ˜ h × cos(2 ψ ) . (4.13b)For a face-on system, ι = h SSB + ( f ) = ˜ A ( f ) e i ( ϕ − ψ ) e − i Ψ ( f ) , (4.14a)˜ h SSB × ( f ) = i ˜ A ( f ) e i ( ϕ − ψ ) e − i Ψ ( f ) , (4.14b)Thus we see that ϕ and ψ appear only through the combination ψ − ϕ yielding a true degeneracy corresponding to ψ − ϕ = const. For systems close to face-on / face-o ff , like the Fiducial system, this gives the strong correlation between ψ and ϕ wellobserved in Fig. 2. For edge-on systems ( ι = π/
2) we haveinstead: ˜ h SSB + ( f ) = ˜ A ( f ) cos 2 ψ e i ϕ e − i Ψ ( f ) , (4.15a)˜ h SSB × ( f ) = ˜ A ( f ) sin 2 ψ e i ϕ e − i Ψ ( f ) , (4.15b)and the degeneracy between ψ and ϕ is then broken as alsoshown in Fig. 12. There we compare the distributions of ψ , ϕ , ι and D L for the Edge-on system to the
Far system (whichis almost face-on). Those systems have similar SNR. When7the degeneracy between ψ and φ is broken, we observe a cor-relation between ϕ and f . This is an artificial correlation dueto relating ϕ to the value of the phase at f for each template.Using a fixed reference frequency, such as the the initial fre-quency of the injected signal for example, eliminates this cor-relation. FIG. 13. Distribution of cos ι and D L using the Flatphys and
Flat-sampl priors for T obs = ff erent, the width of the 90% CI for D L is barely a ff ected and thetrue point, indicated by black lines and squares, is well within the CI. In Fig. 12 we also plot distance and inclination, which showa significant correlation for the
Far system, that is absent forthe
Edge-on system. Distance and inclination are purely ex-trinsic parameters, and the degeneracy features when subdom-inant (higher order) modes are negligible appear in the sameway for LIGO / Virgo and LISA. For LISA, see e.g. a dis-cussion in the context of galactic binaries in [83]. In short,in the limit of face-on / o ff systems the inclination acts as ascaling factor over a rather broad range of inclination values,thus changes in cos ι can be compensated by changes in D L .For close to edge-on systems, the × polarisation of the waveis suppressed (in the wave-frame, before transforming to theSSB frame as in (4.13)). The important point is that this sup-pression of h × depends itself quite rapidly on the inclination,so that reproducing the injected signal leads to a rather tightconstraint on ι , and as a consequence on D L . For MBHBsobservations with LISA, higher modes play an important roleand help breaking these degeneracies [51]; but SBHBs are ob-served by LISA far from coalescence and higher modes arenegligible for these signals.In Fig. 13 we show the e ff ect of the distance prior on theposterior distribution for cos ι and D L using the Flatphys and
Flatsampl priors for T obs = Flatsampl prior, the posterior distribution of cos ι is flat because the likelihood itself is very flat around ι = , π (cos ι is a slowly varying function around its extrema). Thus,the choice of prior shifts the peak of the posterior, but the90% CI still contains the true value and its width is largelyuna ff ected.Among all the cases we have considered, D L can be at bestdetermined within 40%, with the exception of the Edge-on system for which we can determine distance to within 20%.However, the edge-on systems will have lower SNR for a fixeddistance to the source, and, therefore there is a observationalselection e ff ect where the the face-on / o ff systems are preferred(that is what we observe with LIGO / Virgo). If we fix all otherparameters of the
Fiducial system and set ι = π/ − π/ T obs = / SNR.
C. Fisher matrix analysis
In this subsection we consider PE using a slightly improvedversion of Fisher information matrix analysis, inspired by[46]. We have introduced the Fisher matrix in subsection II Band discussed its augmented version, the e ff ective Fisher, insubsection II C for computing the covariance matrix. As wementioned in subsection II C and showed in subsection IV B 2,the likelihood is very flat around ι = , π leading Fisher-based-PE to overestimate the errors on cos ι and D L . To cor-rect for this, we add an additional term ( F t ) to the e ff ectiveFisher matrix: F e ff = F + F p + F t where F is the “original”Fisher matrix given by eq (2.6) and F p is introduced to ac-count for the prior on spins. Empirically, we found the choice F tcos ι, cos ι = . / SNR)) and 0 elsewehere to give good resultsfor cos ι and D L . The prior matrix F p does more than truncat-ing the error on spins: it mimicks the non-trivial prior on χ + and χ − . Indeed, recquiring the spins to be in the physicallyallowed range ( − ≤ χ , ≤
1) leads to a parabola-shapedprior on χ + , − as seen on 1. We approach this non-trivial priorby a Gaussian distribution centered at χ + , − = σ = .
5. We invert the e ff ective Fisher matrix toobtain the covariance matrix and use it to draw points froma multivariate Gaussian distribution. To fully account for thee ff ect of the prior on spins, the point at which the Gaussiandistribution is centred is shifted to: θ e ff = F − ff F θ . We onlykeep points within the boundaries given in Eq. (3.2). For ψ and φ we draw points in an interval of width π around the cen-tral value. In Figs. 14 and 15 we compare our Fisher analysisto the inferred distribution for the Fiducial system using the
Flatsampl prior.We find a very good agreement despite the rather low SNRof this system, especially in the T obs = ff ec-tive Fisher analysis. Naturally, this method cannot reproducethe secondary maximum we found for the Equatorial system8
FIG. 14. Comparison between the inferred distribution for the
Fiducial system using the
Flatsampl prior and our Fisher analysis with T obs = T obs = but it does predict a higher error as the system approaches theequatorial plane. The good agreement for χ + and χ − and for T obs = ff ective Fisher and posterior distri-bution are both prior dominated. For M c and η , Fisher agreeswith the full PE on a 2-sigma level but cannot reproduce thebanana-like correlation. In case of T obs = χ + , reducing the error predictedby the “original” Fisher while the χ − distribution is still priordominated. Without adding F t to the e ff ective Fisher, the di-rection for the correlation between cos ι and D L is predictedwell but the Fisher matrix severely overestimates the error fornearly face-on / face-o ff systems. For the Edgeon system, the9
FIG. 16. Evolution of the error as function of the time before merger we start observing the system in the T obs = T obs = M c , η and χ + correspond to the width of the 90% CIs, and ∆Ω is defined in (4.11). likelihood is not so flat so the error predicted by the “original”Fisher is already small and adding F t does not a ff ect the PE.Based on the rather good agreement we found withBayesian PE, we can exploit the simplicity of Fisher analy-sis to further explore how does the PE evolves with the time(left) to coalescence. In Fig. 16, assuming T obs = T obs = M c , η , χ + and ∆Ω as a func-tion of the time to coalescence t c , keeping all the parametersof the Fiducial system fixed but varying the initial frequencyin accordance with the chosen t c . We plot the correspondingevolution of the SNR in the top panel, with the lowest SNRof 8 being reached for t c (cid:39) t c = T obs in each case which corresponds to the maximum achievableSNR given the observation time and it also corresponds to thebest estimation of parameters. Note two di ff erent regimes onthe two sides of the dashed line: to the left, the PE is governedby the decrease in the signal duration in LISA band and reduc-tion in SNR, while to the right the PE is determined mainly bythe bandwidth of the signal spanned over the observation time.As discussed in Sec. IV B 1 the sky localization comes mainlyfrom modulations caused by the motion of LISA, therefore itworsens rapidly if the system spends too little time in band(below 1 year). D. Long-wavelength approximation
In Table VI we compare the SNR for the
Fiducial , Polar and
Equatorial systems using the
Full and LW responses fortwo observation times. We find that accounting for the degra- T obs = T obs = Full LW Full LWFiducial . . . . Polar . . . . Equatorial . . . . Fiducial , Polar and
Equatorial systems using the
Full and the LW response. dation at high frequencies (Eq. (3.8)), the LW approximationseems to barely a ff ects the PE as can be seen in Fig. 17. Wefind similar behaviour for the Polar and
Equatorial systems.Some care is needed in interpreting this result: this compari-son shows that the high frequency terms neglected in the LWapproximation have little impact on the posterior of the skyposition if the likelihood is computed self-consistently (signaland template are produced using same response, either LW orfull). However, when analysing real data these high frequencyterms cannot be neglected. In other words, these e ff ects in thefull response can indeed be subdominant in the parameter re-covery, if more information comes from other e ff ects like theLISA motion and the main Doppler modulation, while not be-ing negligible in the signal itself. To illustrate this, we simu-late data for the Fiducial , Polar and
Equatorial systems usingthe full response and perform a Bayesian analysis using the LW approximation to compute templates. In Table VII, we0 FIG. 17. Comparison of inferred distributions of intrinsic parmeters and sky location using the
Full and
Lowf responses for the
Fiducial systemin the T obs = T obs = T obs = L ( θ ) max(log L ) ˜ ρ log L ( θ ) max (log L ) ˜ ρ Fiducial − − . − −
38 0 . Polar − − . − −
30 0 . Equatorial − − . − −
34 0 . LW approximationin the Bayesian analysis for data generated with the Full response. give the log-likelihood evaluated at the true point, the maxi-mum likelihood and the maximum overlap:˜ ρ = max h (cid:32) ( d | h ) √ ( d | d )( h | h ) (cid:33) . (4.16)In practice, we compute max by optimizing over our samples.The quantity 1 − ˜ ρ indicates how much SNR would be lost ifwrong templates were used for the detection of signal. We findthat up to ∼
10% of the SNR could be lost, given the alreadylow SNR of SBHBs in LISA this would severely compromiseour chances of detecting such sources. The very small valueof the likelihood at the true point by itself shows that using the LW approximation will have an impact on the PE. In fig. 18,we compare posterior distributions obtained by using templategenerated with Full or LW response while analysing the Fidu-cial system, with T obs = T obs = ff erentresult, the LW template is e ff ectual enough to fit the signalrather well with the largest bias appearing only in ψ − ϕ distri-bution as a compensation for terms neglected in the responseand with a mild drop in the SNR. However, those signals arequite weak and we do not have the luxury to loose even a smallportion of SNR.Thus, our findings seem to validate the LW approximationfor prospective parameter estimation studies, if it is used con-sistently for injecting and recovering the signal, while it wouldbe inappropriate to analyze real data. However, we shouldremember that we did not explore the full parameter space,while (3.8) is valid as an average over orientations, so a di ff er-ent choice of parameters could yield worse results. We alsonote that the full response (3.5) is actually quite simple andnot more expensive computationnally, while being unambigu-ous and eliminating the need for the averaging entering (3.8).1 FIG. 18. Comparison of the infered distributions for the
Fiducial system in the T obs = Full and the LW approximation inthe Bayesian analysis. In both cases, data was generated with the Full response. Black lines and squares indicate the true values.
V. DISCUSSION
Merging binary stellar mass black holes are detected almostweekly during the third LIGO / Virgo observational run (O3).In this work we explored what LISA will be able to tell usabout those binaries. While ground-based detectors observethe last seconds before the merger, LISA will see the early in-spiral evolution of those systems. The results of O3 run arenot publicly available, so we used a GW150914-like systemas a fiducial system in our study. We varied parameters of thesystem in turn, investigating the corresponding changes in thePE. We worked on the simulated (noiseless) data and we usedthe full LISA response. We employed a Bayesian analysis for the parameter estimation and cross-checked our results usingtwo independent samplers. We have found that the PE resultsare most sensitive to the frequency span of the GW and itsposition within LISA sensitivity given the observation dura-tion, or in other words, how much the signal chirps during theobservation time.For weakly chirping systems that do not reach high frequen-cies during LISA’s observations, the GW phase dominated bythe leading PN order, with smaller contributions from higherPN terms. As a result, the best measured parameter is the chirpmass (entering at the leading order) with typical relative errorbelow 10 − . The weak contributions of sub-leading terms upto 1.5PN lead to a three-way correlation between spins, sym-2metric mass ratio and the chirp mass. The mass ratio is verypoorly constrained and the posterior for the spins is dominatedby the priors. We still can recover very well the sky positionof the source (typically < . ) which comes from themodulation of the amplitude and the phase of the GW signaldue to LISA’s motion. Such an area in the sky is within thefield of view of electromagnetic instruments, such as Athenaand / or SKA.For chirping systems, that coalesce during the observa-tions and therefore reach the high end of the LISA frequencyband, higher order PN terms become more important andhelp breaking the correlations in intrinsic parameters, lead-ing to a significant improvement in parameter estimation. Theindividual masses for chirping systems are measured within20 −
30% and even better for systems with higher mass ratio.The constraints on individual spins result from the combina-tion of the measurement of χ PN and the physical boundaries ofthe prior on spins: − ≤ χ , ≤
1. This suggests suggest using χ PN as a sampling parameter. We find that the measurement ofthe time to coalescence improves as we observe the systemscloser to merger, from O (1 day) (for mildly chirping binaries)to O (30 s). We note that the only way to increase our chancesto observe SBHBs chirping is by increasing LISA’s missionduration.The measurement of the luminosity distance is less im-pacted by whether the systems are chirping or not, and it isessentially a function of their SNR. Much like LIGO / Virgoobservations when higher modes can be neglected, the degen-eracy between distance and inclination is important. In ourexample, the distance is typically measured within 40 − / o ff and within 20% if thesystem is edge-on (when adjusting the distance to keep thesame SNR). The distance and therefore redshift uncertaintydominates the measurement of the source-frame chirp massat a percent level. The precision on individual masses in thesource frame is dominated by intrinsic degeneracies.We have suggested an augmentation of the usual Fisher ma-trix approach, that we called the e ff ective Fisher matrix, andwe have shown that it gives rather reliable results for the skyposition and intrinsic parameters of the system when com-pared to Bayesian PE. We also showed that combining the useof the long-wavelength approximation for LISA with the in-troduction of a degradation factor at high frequencies yieldsvery similar results as compared to using the full responsefor computing likelihoods self-consistently (using the sameresponse for injected data and templates). However, usingthe long-wavelength approximation to analyse real data coulddecrease the SNR by 10%, drastically reducing our chancesof detecting the signals and has a significant impact on thePE, in particular on the measurement of intrinisc parameters.The computational cost of the full response being essentiallythe same than the long-wavelength approximation, we recom-mend its use in future works.We can utilise the knowledge and understanding obtained inPE study for development of the search: (i) the PE for thosesystems is mainly monomodal, the secondary modes appeareither in special cases (like Equatorial system considered inthe main body of the paper) or under the e ff ect of prior for pa- rameters on which the likelihood is weakly informative (like η for non-chirping systems) (ii) the chirp mass and the sky co-ordinates are the best measured parameters, so we can makehierarchical search starting with those parameters taking intoaccount the correlations which are explored and understood(see section IV) (iii) the e ff ective Fisher is a good proposalfor Bayes-based search once we started to see sign of the GWsignal in the data. In addition we can perform data volume-increasing analysis starting with half year long data segmentprogressively increasing it. This works as a natural annealingscheme and should help in recovering (especially) chirpingsystems.Detection of even few SBHBs by LISA which merge some-what later with a very high SNR in the band of ground baseddetectors [37] will constitute “golden events”. Beyond all thebenefits of multiband detections per se, the information pro-vided by LISA itself will be very valuable. For example, [84]suggested that the measurement of the time to coalescencecould be used to inform ground based detectors and improveblack hole spectroscopy. The good estimate on the time to co-alescence and the sky location could be used for electromag-netic follow-up of the source as suggested in [27]. Finally,these measurements could be used to tighten the constraintson the Hubble constant ( H ) even if no electromagnetic coun-terparts are detected, using galaxy catalogues [85–87].Modifications to the GW phase induced by either modifiedtheories of gravity or environmental e ff ects will genericallyinvolve additional coe ffi cients parametrising the underlyingmechanisms and their correlation with the parameters of thesystem ( M c , η , . . . ). As a consequence, even for low fre-quency modifications the best constraints / measurements willcome from chirping systems as we found in [21] in the con-text of testing modified theories of gravity with LISA obser-vations.A major improvement to this work would be the inclusionof eccentricity. Astrophysical formation models predict thatbinaries formed dynamically should have large eccentricities[35, 36]. However, by the time these binaries reach the fre-quency band of ground based detectors they will have cir-cularised. Thus LISA could play an important role in thediscrimination between di ff erent formation channels [40–43].Furthermore, neglecting eccentricity could a ff ect PE and de-tection e ffi ciency. We are currently limited by the lack of fasteccentric waveforms but there is a good progress in this direc-tion [88–92]. Concerning spins, binaries formed dynamicallyare expected to have misaligned spins [37] but precession ef-fects should be weak in the early inspiral. Thus, althoughinteresting, the use of precessing waveforms should not ap-preciably change conclusions drawn in this work.We conclude with the claim that this work, together with[21, 27], has confirmed the scientific potential of the observa-tion of SBHBs with LISA and should be seen as a first steptowards an extensive study of the PE for the multiband obser-vations.3 ACKNOWLEDGEMENTS
A.T. is grateful to Nikolaos Karnesis for his help at the startof the project and his support all along it. This work hasbeen supported by the European Union’s Horizon 2020 re- search and innovation program under the Marie Skłodowska-Curie grant agreement No 690904. A.T. acknowledges finan-cial support provided by Paris-Diderot University (now part ofUniversit´e de Paris). The authors would also like to acknowl-edge networking support from the COST Action CA16104. [1] B. P. Abbott et al. Observation of Gravitational Waves froma Binary Black Hole Merger.
Phys. Rev. Lett. , 116(6):061102,2016.[2] B.P. Abbott et al. GWTC-1: A Gravitational-Wave TransientCatalog of Compact Binary Mergers Observed by LIGO andVirgo during the First and Second Observing Runs.
Phys. Rev.X , 9(3):031040, 2019.[3] B. P. Abbott et al. GW170814: A Three-Detector Observationof Gravitational Waves from a Binary Black Hole Coalescence.
Phys. Rev. Lett. , 119(14):141101, 2017.[4] B.. P.. Abbott et al. GW170608: Observation of a 19-solar-massBinary Black Hole Coalescence.
Astrophys. J. , 851(2):L35,2017.[5] B. P. Abbott et al. Binary Black Hole Mergers in the first Ad-vanced LIGO Observing Run.
Phys. Rev. , X6(4):041015, 2016.[erratum: Phys. Rev.X8,no.3,039903(2018)].[6] B.P. Abbott et al. GW170817: Observation of GravitationalWaves from a Binary Neutron Star Inspiral.
Phys. Rev. Lett. ,119(16):161101, 2017.[7] Benjamin P. Abbott et al. GW170104: Observation of a50-Solar-Mass Binary Black Hole Coalescence at Redshift0.2.
Phys. Rev. Lett. , 118(22):221101, 2017. [Erratum:Phys.Rev.Lett. 121, 129901 (2018)].[8] R. Abbott et al. GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses. 4 2020.[9] B.P. Abbott et al. GW190425: Observation of a Compact Bi-nary Coalescence with Total Mass ∼ . M (cid:12) . Astrophys. J. Lett. ,892:L3, 2020.[10] R. Abbott et al. GW190814: Gravitational waves from the co-alescence of a 23 solar mass black hole with a 2.6 solar masscompact object.
The Astrophysical Journal , 896(2):L44, jun2020.[11] B. P. Abbott et al. Binary Black Hole Population Properties In-ferred from the First and Second Observing Runs of AdvancedLIGO and Advanced Virgo.
Astrophys. J. , 882(2):L24, 2019.[12] Benjamin P. Abbott et al. Upper Limits on the Rates of BinaryNeutron Star and Neutron Star–black Hole Mergers From Ad-vanced Ligo’s First Observing run.
Astrophys. J. , 832(2):L21,2016.[13] B. P. Abbott et al. Tests of general relativity with the binaryblack hole signals from the ligo-virgo catalog gwtc-1. 2019.[14] B. P. Abbott et al. Tests of general relativity with GW150914.
Phys. Rev. Lett. , 116(22):221101, 2016. [Erratum: Phys. Rev.Lett.121,no.12,129902(2018)].[15] B. P. Abbott et al. Tests of General Relativity with GW170817.
Phys. Rev. Lett. , 123(1):011102, 2019.[16] Heather Audley et al. Laser Interferometer Space Antenna.2017.[17] Antoine Klein et al. Science with the space-based interfer-ometer eLISA: Supermassive black hole binaries.
Phys. Rev. ,D93(2):024003, 2016.[18] Valeriya Korol, Elena M. Rossi, Paul J. Groot, Gijs Nelemans,Silvia Toonen, and Anthony G.A. Brown. Prospects for detec-tion of detached double white dwarf binaries with Gaia, LSST and LISA.
Mon. Not. Roy. Astron. Soc. , 470(2):1894–1910,2017.[19] Alberto Sesana. Prospects for Multiband Gravitational-WaveAstronomy after GW150914.
Phys. Rev. Lett. , 116(23):231102,2016.[20] Pau Amaro-Seoane and Lucia Santamaria. Detection of IMBHswith ground-based gravitational wave observatories: A biogra-phy of a binary of black holes, from birth to death.
Astrophys.J. , 722:1197–1206, 2010.[21] Alexandre Toubiana, Sylvain Marsat, Enrico Barausse,Stanislav Babak, and John Baker. Tests of general relativitywith stellar-mass black hole binaries observed by LISA. 2020.[22] Giuseppe Gnocchi, Andrea Maselli, Tiziano Abdelsalhin,Nicola Giacobbo, and Michela Mapelli. Bounding AlternativeTheories of Gravity with Multi-Band GW Observations. 2019.[23] Zack Carson and Kent Yagi. Multi-band gravitational wavetests of general relativity. 2019.[24] Salvatore Vitale. Multiband Gravitational-Wave Astron-omy: Parameter Estimation and Tests of General Relativitywith Space- and Ground-Based Detectors.
Phys. Rev. Lett. ,117(5):051102, 2016.[25] Enrico Barausse, Nicol´as Yunes, and Katie Chamberlain.Theory-Agnostic Constraints on Black-Hole Dipole Radiationwith Multiband Gravitational-Wave Astrophysics.
Phys. Rev.Lett. , 116(24):241104, 2016.[26] Katie Chamberlain and Nicolas Yunes. Theoretical Physics Im-plications of Gravitational Wave Observation with Future De-tectors.
Phys. Rev. , D96(8):084039, 2017.[27] Andrea Caputo, Laura Sberna, Alexandre Toubiana, StanislavBabak, Enrico Barausse, Sylvain Marsat, and Paolo Pani.Gravitational-wave detection and parameter estimation for ac-creting black-hole binaries and their electromagnetic counter-part. 2020.[28] Enrico Barausse, Vitor Cardoso, and Paolo Pani. Can environ-mental e ff ects spoil precision gravitational-wave astrophysics? Phys. Rev. , D89(10):104059, 2014.[29] Enrico Barausse, Vitor Cardoso, and Paolo Pani. EnvironmentalE ff ects for Gravitational-wave Astrophysics. J. Phys. Conf. Ser. ,610(1):012044, 2015.[30] Nicola Tamanini, Antoine Klein, Camille Bonvin, Enrico Ba-rausse, and Chiara Caprini. The peculiar acceleration of stellar-origin black hole binaries: measurement and biases with LISA.2019.[31] Vitor Cardoso and Andrea Maselli. Constraints on the astro-physical environment of binaries with gravitational-wave ob-servations. 2019.[32] Konstantin A. Postnov and Lev R. Yungelson. The Evolutionof Compact Binary Star Systems.
Living Rev. Rel. , 17:3, 2014.[33] Matthew J. Benacquista and Jonathan M. B. Downing. Rela-tivistic Binaries in Globular Clusters.
Living Rev. Rel. , 16:4,2013.[34] Simeon Bird, Ilias Cholis, Julian B. Mu˜noz, Yacine Ali-Ha¨ımoud, Marc Kamionkowski, Ely D. Kovetz, Alvise Rac-canelli, and Adam G. Riess. Did LIGO detect dark matter? Phys. Rev. Lett. , 116(20):201301, 2016.[35] Fabio Antonini and Hagai B. Perets. Secular evolution ofcompact binaries near massive black holes: Gravitational wavesources and other exotica.
Astrophys. J. , 757:27, 2012.[36] Johan Samsing. Eccentric Black Hole Mergers Forming inGlobular Clusters.
Phys. Rev. , D97(10):103014, 2018.[37] Davide Gerosa, Emanuele Berti, Richard O’Shaughnessy,Krzysztof Belczynski, Michael Kesden, Daniel Wysocki, andWojciech Gladysz. Spin orientations of merging black holesformed from the evolution of stellar binaries.
Phys. Rev. ,D98(8):084036, 2018.[38] Vishal Baibhav, Davide Gerosa, Emanuele Berti, Kaze W. K.Wong, Thomas Helfer, and Matthew Mould. The mass gap, thespin gap, and the origin of merging binary black holes. 2020.[39] Davide Gerosa, Sizheng Ma, Kaze W. K. Wong, EmanueleBerti, Richard O’Shaughnessy, Yanbei Chen, and KrzysztofBelczynski. Multiband gravitational-wave event rates and stel-lar physics.
Phys. Rev. , D99(10):103004, 2019.[40] Johan Samsing and Daniel J. D’Orazio. Black Hole Merg-ers From Globular Clusters Observable by LISA I: EccentricSources Originating From Relativistic N -body Dynamics. Mon.Not. Roy. Astron. Soc. , 481(4):5445–5450, 2018.[41] Atsushi Nishizawa, Emanuele Berti, Antoine Klein, and Al-berto Sesana. eLISA eccentricity measurements as tracers ofbinary black hole formation.
Phys. Rev. , D94(6):064020, 2016.[42] Atsushi Nishizawa, Alberto Sesana, Emanuele Berti, and An-toine Klein. Constraining stellar binary black hole formationscenarios with eLISA eccentricity measurements.
Mon. Not.Roy. Astron. Soc. , 465(4):4375–4380, 2017.[43] Katelyn Breivik, Carl L. Rodriguez, Shane L. Larson, VassilikiKalogera, and Frederic A. Rasio. Distinguishing Between For-mation Channels for Binary Black Holes with LISA.
Astrophys.J. , 830(1):L18, 2016.[44] Seth E. Timpano, Louis J. Rubbo, and Neil J. Cornish. Charac-terizing the galactic gravitational wave background with LISA.
Phys. Rev. , D73:122001, 2006.[45] Christopher J. Moore, Davide Gerosa, and Antoine Klein. Arestellar-mass black-hole binaries too quiet for LISA?
Mon. Not.Roy. Astron. Soc. , 488(1):L94–L98, 2019.[46] Michele Vallisneri. Use and abuse of the Fisher informa-tion matrix in the assessment of gravitational-wave parameter-estimation prospects.
Phys. Rev. , D77:042001, 2008.[47] Curt Cutler. Angular resolution of the LISA gravitational wavedetector.
Phys. Rev. , D57:7089–7102, 1998.[48] Alberto Vecchio and Elizabeth D. L. Wickham. LISA re-sponse function and parameter estimation.
Class. Quant. Grav. ,21(5):S661–S664, 2004.[49] Alberto Vecchio and Elizabeth D. L. Wickham. The E ff ect ofthe LISA response function on observations of monochromaticsources. Phys. Rev. , D70:082002, 2004.[50] Sylvain Marsat and John G. Baker. Fourier-domain modula-tions and delays of gravitational-wave signals. 6 2018.[51] Sylvain Marsat, John G. Baker, and Tito Dal Canton. Exploringthe Bayesian parameter estimation of binary black holes withLISA. 2020.[52] S Hild, M Abernathy, F Acernese, P Amaro-Seoane, N Ander-sson, K Arun, F Barone, B Barr, M Barsuglia, M Beker, andet al. Sensitivity studies for third-generation gravitational waveobservatories.
Classical and Quantum Gravity , 28(9):094013,Apr 2011.[53] M. Punturo et al. The Einstein Telescope: A third-generation gravitational wave observatory.
Class. Quant. Grav. ,27:194002, 2010. [54] Stefan Ballmer and Vuk Mandic. New Technologies inGravitational-Wave Detection.
Ann. Rev. Nucl. Part. Sci. ,65:555–577, 2015.[55] B P Abbott, R Abbott, T D Abbott, M R Abernathy, K Ackley,C Adams, P Addesso, R X Adhikari, V B Adya, C A ff eldt, andet al. Exploring the sensitivity of next generation gravitationalwave detectors. Classical and Quantum Gravity , 34(4):044001,Jan 2017.[56] Johan Samsing and Daniel J. D’Orazio. How Post-NewtonianDynamics Shape the Distribution of Stationary Binary BlackHole LISA Sources in nearby Globular Clusters.
Phys. Rev. ,D99(6):063006, 2019.[57] Kyle Kremer et al. Post-Newtonian Dynamics in Dense StarClusters: Binary Black Holes in the LISA Band.
Phys. Rev. ,D99(6):063003, 2019.[58] Kaze W. K. Wong, Ely D. Kovetz, Curt Cutler, and EmanueleBerti. Expanding the LISA Horizon from the Ground.
Phys.Rev. Lett. , 121(25):251102, 2018.[59] Koutarou Kyutoku and Naoki Seto. Concise estimate of the ex-pected number of detections for stellar-mass binary black holesby eLISA.
Mon. Not. Roy. Astron. Soc. , 462(2):2177–2183,2016.[60] Bruce Allen, Warren G. Anderson, Patrick R. Brady, Duncan A.Brown, and Jolien D. E. Creighton. FINDCHIRP: An Algo-rithm for detection of gravitational waves from inspiraling com-pact binaries.
Phys. Rev. , D85:122006, 2012.[61] Sascha Husa, Sebastian Khan, Mark Hannam, Michael P¨urrer,Frank Ohme, Xisco Jim´enez Forteza, and Alejandro Boh´e.Frequency-domain gravitational waves from nonprecessingblack-hole binaries. I. New numerical waveforms and anatomyof the signal.
Phys. Rev. , D93(4):044006, 2016.[62] Sebastian Khan, Sascha Husa, Mark Hannam, Frank Ohme,Michael P¨urrer, Xisco Jim´enez Forteza, and Alejandro Boh´e.Frequency-domain gravitational waves from nonprecessingblack-hole binaries. II. A phenomenological model for the ad-vanced detector era.
Phys. Rev. , D93(4):044007, 2016.[63] Massimo Tinto and Sanjeev V. Dhurandhar. Time-delay inter-ferometry.
Living Rev. Rel. , 8:4, 2005.[64] Luc Blanchet. Gravitational radiation from postNewtoniansources and inspiraling compact binaries.
Living Rev. Rel. , 5:3,2002.[65] Alessandra Buonanno, Gregory B. Cook, and Frans Pretorius.Inspiral, merger and ring-down of equal-mass black-hole bina-ries.
Phys. Rev. , D75:124018, 2007.[66] Alessandra Buonanno, Bala Iyer, Evan Ochsner, Yi Pan, andB. S. Sathyaprakash. Comparison of post-Newtonian templatesfor compact binary inspiral signals in gravitational-wave detec-tors.
Phys. Rev. , D80:084043, 2009.[67] Rajeeva L. Karandikar. On the markov chain monte carlo(mcmc) method.
Sadhana , 31(2):81–104, Apr 2006.[68] Understanding the metropolis-hastings algorithm.
The Ameri-can Statistician , 49(4):327–335, 1995.[69] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H.Teller, and E. Teller. Equation of state calculations by fast com-puting machines.
J. Chem. Phys. , 21:1087–1092, 1953.[70] Andrew Gelman and Donald B. Rubin. Inference from IterativeSimulation Using Multiple Sequences.
Statist. Sci. , 7:457–472,1992.[71] Thomas D. Abbott et al. Improved analysis of GW150914using a fully spin-precessing waveform Model.
Phys. Rev. ,X6(4):041014, 2016.[72] N. Aghanim et al. Planck 2018 results. VI. Cosmological pa-rameters. 2018.[73] L. S. S. Team. LISA Science Requirements Document. [74] Travis Robson, Neil J. Cornish, and Chang Liu. The construc-tion and use of LISA sensitivity curves. Class. Quant. Grav. ,36(10):105011, 2019.[75] S. Kullback and R. A. Leibler. On information and su ffi ciency. Ann. Math. Statist. , 22(1):79–86, 03 1951.[76] Daniel Foreman-Mackey. corner.py: Scatterplot matrices inpython.
The Journal of Open Source Software , 24, 2016.[77] Alberto Mangiagli, Antoine Klein, Alberto Sesana, Enrico Ba-rausse, and Monica Colpi. Post-Newtonian phase accuracy re-quirements for stellar black hole binaries with LISA.
Phys. Rev. ,D99(6):064056, 2019.[78] Eric Poisson and Cli ff ord M. Will. Gravitational waves from in-spiraling compact binaries: Parameter estimation using secondpostNewtonian wave forms. Phys. Rev. , D52:848–855, 1995.[79] Ska white paper, 2014.[80] N. Meidinger. The Wide Field Imager instrument for Athena.
Contributions of the Astronomical Observatory Skalnate Pleso ,48(3):498–505, Jul 2018.[81] Lee Samuel Finn and David F. Cherno ff . Observing binary in-spiral in gravitational radiation: One interferometer. Phys. Rev.D , 47:2198–2219, 1993.[82] Curt Cutler and Eanna E. Flanagan. Gravitational waves frommerging compact binaries: How accurately can one extract thebinary’s parameters from the inspiral wave form?
Phys. Rev. ,D49:2658–2697, 1994.[83] Sweta Shah, Marc van der Sluys, and Gijs Nelemans. Usingelectromagnetic observations to aid gravitational-wave parame-ter estimation of compact binaries observed with LISA.
Astron.Astrophys. , 544:A153, 2012.[84] Rhondale Tso, Davide Gerosa, and Yanbei Chen. Optimiz-ing LIGO with LISA forewarnings to improve black-hole spec- troscopy.
Phys. Rev. , D99(12):124043, 2019.[85] Bernard F. Schutz. Determining the Hubble Constant fromGravitational Wave Observations.
Nature , 323:310–311, 1986.[86] Walter Del Pozzo, Alberto Sesana, and Antoine Klein. Stellarbinary black holes in the LISA band: a new class of standardsirens.
Mon. Not. Roy. Astron. Soc. , 475(3):3485–3492, 2018.[87] Koutarou Kyutoku and Naoki Seto. Gravitational-wave cos-mography with LISA and the Hubble tension.
Phys. Rev. ,D95(8):083525, 2017.[88] Zhoujian Cao and Wen-Biao Han. Waveform model for aneccentric binary black hole based on the e ff ective-one-body-numerical-relativity formalism. Phys. Rev. D , 96(4):044028,2017.[89] Brennan Ireland, Ofek Birnholtz, Hiroyuki Nakano, Eric West,and Manuela Campanelli. Eccentric Binary Black Holes withSpin via the Direct Integration of the Post-Newtonian Equationsof Motion.
Phys. Rev. D , 100(2):024015, 2019.[90] Ian Hinder, Lawrence E. Kidder, and Harald P. Pfei ff er. Eccen-tric binary black hole inspiral-merger-ringdown gravitationalwaveform model from numerical relativity and post-Newtoniantheory. Phys. Rev. D , 98(4):044015, 2018.[91] Tanja Hinderer and Stanislav Babak. Foundations of ane ff ective-one-body model for coalescing binaries on eccentricorbits. Phys. Rev. D , 96(10):104048, 2017.[92] E.A. Huerta et al. Eccentric, nonspinning, inspiral, Gaussian-process merger approximant for the detection and characteri-zation of eccentric binary black hole mergers.