Systematic biases in parameter estimation on LISA binaries. II. The effect of excluding higher harmonics for spin-aligned, high-mass binaries

TL;DR

Massive black hole binaries observed by LISA can suffer sizable systematic biases when higher-order waveform modes are neglected, especially for high-mass systems with aligned spins. The authors extend prior work to cover total detector-frame masses up to $M=10^8\,M_\odot$ and reveal that higher-order modes can dominate the signal in parts of parameter space, with biases strongly depending on $q$, $\iota$, and spins; sky localization can be severely biased for the heaviest, shortest signals. They develop an improved likelihood-optimization workflow, combining dual annealing, reparameterization, and Fisher-informed priors, to robustly predict these biases in a computationally efficient manner. They also analyze sky-position degeneracies (octants) and show how multimodal likelihoods can be navigated or mitigated, while clearly indicating regimes where full Bayesian PE remains necessary. The study underscores the necessity of accurate higher-mode modeling for LISA MBHB science and provides practical strategies for rapid, robust bias estimation in planning analyses.

Abstract

The Laser Interferometer Space Antenna (LISA) will observe massive black hole binaries (MBHBs) with astoundingly high signal-to-noise ratio, leaving parameter estimation with these signals susceptible to seemingly small waveform errors. Of particular concern for MBHBs are errors due to neglected higher-order modes. We extend Yi et al. [arXiv:2502.12237] to examine errors due to neglected higher-order modes for MBHBs with nonzero (aligned) progenitor spins and total mass up to $10^8\,M_\odot$. For these very massive systems, there can be regions of parameter space in which the $(\ell, |m|)=(2,\,2)$ modes are no longer dominant with respect to higher-order ones. We find that the extent of systematic bias can change significantly when varying the progenitor spins of the binary. We also find that for the heaviest, and therefore shortest, MBHB signals, slight systematic errors can cause severe mis-inference of the sky localization parameters. We propose an improved likelihood optimization scheme with respect to previous work as a way to predict these effects in a computationally efficient manner.

Figures (14)

• Figure 1: SNR vs. detector-frame total mass for a few values of mass ratio, $q$, and inclination angle, $\iota$, at redshift $z=1$ ($D_L=6791.8$ Mpc). The lines are shown for nonspinning binaries ($\chi_1=\chi_2=0$). The remaining extrinsic parameters are the same as in Yi:2025pxe (namely, [$\delta t,\,\phi,\,\lambda_L,\,\beta_L,\,\Psi_L$] = [0, 0.2, 1.8, $\pi/6$, 1.2]).
• Figure 2: Hierarchy of SNR contribution by different angular harmonics as a function of total mass ($y$-axes), mass ratio ($x$-axes), inclination angle (column), and progenitor spins (row). As in Yi:2025pxe, we use log spacing between $q=1$ and $q=3$, as the importance of higher-order modes tends to change more rapidly between these mass ratios. The spacing is linear above $q=3$ (marked by the dashed gray line). Compared to Yi:2025pxe, the individual mode SNRs are ordered in many more different ways across this much broader region of parameter space including significantly higher masses and nonzero spins. Most notably, in the purple, blue, and darker green regions, the (2, 2) mode contributes less SNR than the (3, 3) and/or (4, 4) modes. Finally, the presence of a galactic white dwarf background causes nontrivial changes in the relative mode contribution around $M\sim10^7\,M_\odot$ ($10^{-4}\lesssim f_{\rm ISCO}/\rm{Hz}\lesssim 10^{-3}$). A comparison with Fig. \ref{['fig:mode_by_mode_noWD']} of Appendix \ref{['app:frac_snr']} clarifies that the features observed here are indeed due to the galactic background.
• Figure 3: Hierarchy of SNR contribution by different angular harmonics as a function of mass ratio. At the top, we show the hierarchy for a reference system of total mass $10^7\,M_\odot$, inclination $\pi/3$, and nonspinning progenitors. In the lower panels, we examine how this hierarchy changes as we modify one set of parameters at a time compared to the reference system (i.e., varying the inclination, total mass, and spins). In each panel, we show the total SNR in a thick light gray line.
• Figure 4: Same as Fig. \ref{['fig:spin_dominance_condensed']}, but for the cross-term contributions to the signal. We again show the total SNR in a thick light gray line for reference. To improve readability, in this figure we denote $\left(h_{\ell m}|h_{\ell ' m'}\right)$ as simply $\left(\ell m|\ell'm'\right)$.
• Figure 5: Comparison of the systematic biases on two MBHB parameters with 4 modes $=\{(2,\,2),\,(2,\,1),\,(3,\,3),\,(4,\,4)\}$ (excluding the $(3,\,2)$ mode), with biases determined by a full PE run (green violin plots and quantiles), the Cutler-Vallisneri approach (magenta squares), and the original likelihood optimization procedure used in Yi:2025pxe (blue circles labeled "NM" for the Nelder-Mead algorithm). The systematic biases are normalized by the respective statistical errors on each parameter. In these cases, simple likelihood optimization via the Nelder-Mead algorithm continues to work well in estimating the systematic biases due to excluding higher-order modes for heavier events than were considered in Ref. Yi:2025pxe. Note that one of the systems considered here has nonzero aligned progenitor spins ($\chi_1=0.5,\chi_2=-0.5$).