Massive Black Hole Binaries Research Articles

On the formation of GW190521-like binary black hole merger systems

ABSTRACT GW190521 is the most massive merging binary black hole (BBH) system detected so far. At least one of the component BHs was measured to lie within the pair-instability supernova (PISN) mass gap (∼50–135 M⊙), making its formation a mystery. However, the transient observed signal allows alternative posterior distributions. There was suggestion that GW190521 could be an intermediate-mass ratio inspiral, with the component masses m1 ∼ 170 M⊙ and m2 ∼ 16 M⊙, happening to straddle the PISN mass gap. Under this framework, we perform binary population synthesis to explore the formation of GW190521-like systems via isolated binary evolution. We numerically calculate the binding energy parameter for massive stars at different metallicities, and employ them in our calculation for common envelope evolution. Our results prefer that the progenitor binaries formed in metal-poor environment with $\rm Z\le 0.0016$. The predicted merger rate density within redshift z = 1.1 is ${\sim} 4\times 10^{-5}{\text {--}}5\times 10^{-2} \, \rm Gpc^{-3}yr^{-1}$. We expect that such events are potentially observable by upcoming both space and ground-based gravitational wave detectors.

Can gravitational-wave memory help constrain binary black-hole parameters? A LISA case study

Besides the transient effect, the passage of a gravitational wave also causes a persistent displacement in the relative position of an interferometer's test masses through the nonlinear memory effect. This effect is generated by the gravitational backreaction of the waves themselves, and encodes additional information about the source. In this work, we explore the implications of using this information for the parameter estimation of massive binary black holes with LISA. Based on a Fisher analysis for nonprecessing black hole binaries, our results show that the memory can help to reduce the degeneracy between the luminosity distance and the inclination for binaries observed only for a short time ($\ensuremath{\sim}\mathrm{few}\text{ }\mathrm{hours}$) before merger. To assess how many such short signals will be detected, we utilized state-of-the-art predictions for the population of massive black hole binaries and models for the gaps expected in the LISA data. We forecast from tens to few hundreds of binaries with observable memory, but only $\ensuremath{\sim}\mathcal{O}(0.1)$ events in 4 years for which the memory helps to reduce the degeneracy between distance and inclination. Based on this, we conclude that the new information from the nonlinear memory, while promising for testing general relativity in the strong field regime, has probably a limited impact on further constraining the uncertainty on massive black hole binary parameters with LISA.

Post-Keplerian waveform model for binary compact object as sources of space-based gravitational wave detector and its implications

Binary compact objects will be among the important sources for the future space-based gravitational wave detectors. Such binary compact objects include stellar massive binary black hole, binary neutron star, binary white dwarf and mixture of these compact objects. Regarding to the relatively low frequency, the gravitational interaction between the two objects of the binary is weak. Post-Newtonian approximation of general relativity is valid. Previous works about the waveform model for such binaries in the literature consider the dynamics for specific situations which involve detailed complicated matter dynamics between the two objects. We here take a different idea. We adopt the trick used in pulsar timing detection. For any gravity theories and any detailed complicated matter dynamics, the motion of the binary can always be described as a post-Keplerian expansion. And a post-Keplerian gravitational waveform model will be reduced. Instead of object masses, spins, matter’s equation of state parameters and dynamical parameters beyond general relativity, the involved parameters in our post-Keplerian waveform model are the Keplerian orbit elements and their adiabatic variations. Respect to current planning space-based gravitational wave detectors including LISA, Taiji and Tianqin, we find that the involved waveform model parameters can be well determined. And consequently the detail matter dynamics of the binary can be studied then. For binary with purely gravitational interactions, gravity theory can be constrained well.

Probing z ≳ six massive black holes with gravitational waves

ABSTRACT We investigate the coalescence of massive black hole ($M_{\rm BH}\gtrsim 10^{6}~\rm {\rm M}_{\odot }$) binaries (MBHBs) at 6 < z < 10 by adopting a suite of cosmological hydrodynamical simulations of galaxy formation, zoomed-in on biased (>3σ) overdense regions (Mh ∼ 1012 M⊙ dark matter haloes at z = 6) of the Universe. We first analyse the impact of different resolutions and AGN feedback prescriptions on the merger rate, assuming instantaneous mergers. Then, we compute the halo bias correction factor due to the overdense simulated region. Our simulations predict merger rates that range between 3 and 15 $\rm yr^{-1}$ at z ∼6, depending on the run considered, and after correcting for a bias factor of ∼20−30. For our fiducial model, we further consider the effect of delay in the MBHB coalescence due to dynamical friction. We find that 83 per cent of MBHBs will merge within the Hubble time, and 21 per cent within 1 Gyr, namely the age of the Universe at z > 6. We finally compute the expected properties of the gravitational wave (GW) signals and find the fraction of LISA detectable events with high signal-to-noise ratio (SNR > 5) to range between 66 per cent and 69 per cent. However, identifying the electro-magnetic counterpart of these events remains challenging due to the poor LISA sky localization that, for the loudest signals ($\mathcal {M}_c\sim 10^6~{{\rm M}_{\odot }}$ at z = 6), is around 10 $\rm deg^2$.

Orbital evolution of binaries in circumbinary discs

ABSTRACT We present the to-date largest parameter space exploration of binaries in circumbinary discs (CBDs), deriving orbital evolution prescriptions for eccentric, unequal mass binaries from our suite of hydrodynamic simulations. In all cases, binary eccentricities evolve towards steady state values that increase with mass ratio, and saturate at an equilibrium eccentricity eb,eq ∼ 0.5 in the large mass ratio regime, in line with resonant theory. For binaries accreting at their combined Eddington limit, a steady state eccentricity can be achieved within a few megayears. Once at their steady state eccentricities, binaries with qb ≳ 0.3 evolve towards coalescence, while lower mass ratio systems expand due to CBD torques. We discuss implications for population studies of massive black hole binaries, protostars in binary systems, and post-common envelope binaries observed by ground-based gravitational wave detectors.

Hypercritical accretion during common envelopes in triples leading to binary black holes in the pair-instability-supernova mass gap

ABSTRACT Hydrodynamic studies of stellar-mass compact objects (COs) in a common envelope (CE) have shown that the accretion rate onto the CO is a few orders of magnitude below the Bondi–Hoyle–Lyttleton (BHL) estimate. This is several orders of magnitude above the Eddington limit and above the limit for neutrino-cooled accretion (i.e. hypercritical accretion or HCA). Considering that a binary system inside the CE of a third star accretes material at nearly the same rate as a single object of the same total mass, we propose stellar-evolution channels which form binary black hole (BBH) systems with its component masses within the pair-instability supernova (PISN) mass gap. Our model is based on HCA onto the BBH system engulfed into the CE of a massive tertiary star. Furthermore, we propose a mass transfer mode which allows to store mass lost by the binary onto a third star. Through the use of population synthesis simulations for the evolution of BBHs and standard binary-evolution principles for the interaction with a tertiary star, we are able to produce BBHs masses consistent with those estimated for GW190521. We also discuss the massive binary system Mk34 as a possible progenitor of BBHs in the PISN gap.

Rapid search for massive black hole binary coalescences using deep learning

The coalescences of massive black hole binaries are one of the main targets of space-based gravitational wave observatories. Such gravitational wave sources are expected to be accompanied by electromagnetic emissions. Low latency detection of the massive black hole mergers provides a start point for a global-fit analysis to explore the large parameter space of signals simultaneously being present in the data but at great computational cost. To alleviate this issue, we present a deep learning method for rapidly searching for signals of massive black hole binaries in gravitational wave data. Our model is capable of processing a year of data, simulated from the LISA data challenge, in only several seconds, while identifying all coalescences of massive black hole binaries with no false alarms. We further demonstrate that the model shows robust resistance to a wide range of generalization cases, including various waveform families and updated instrumental configurations. This method offers an effective approach that combines advances in artificial intelligence to open a new pathway for space-based gravitational wave observations.

Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3

We report on the population properties of compact binary mergers inferred from gravitational-wave observations of these systems during the first three LIGO-Virgo observing runs. The Gravitational-Wave Transient Catalog 3 (GWTC-3) contains signals consistent with three classes of binary mergers: binary black hole, binary neutron star, and neutron star–black hole mergers. We infer the binary neutron star merger rate to be between 10 and 1700 Gpc−3 yr−1 and the neutron star–black hole merger rate to be between 7.8 and 140 Gpc−3 yr−1, assuming a constant rate density in the comoving frame and taking the union of 90% credible intervals for methods used in this work. We infer the binary black hole merger rate, allowing for evolution with redshift, to be between 17.9 and 44 Gpc−3 yr−1 at a fiducial redshift (z=0.2). The rate of binary black hole mergers is observed to increase with redshift at a rate proportional to (1+z)κ with κ=2.9−1.8+1.7 for z≲1. Using both binary neutron star and neutron star–black hole binaries, we obtain a broad, relatively flat neutron star mass distribution extending from 1.2−0.2+0.1 to 2.0−0.3+0.3M⊙. We confidently determine that the merger rate as a function of mass sharply declines after the expected maximum neutron star mass, but cannot yet confirm or rule out the existence of a lower mass gap between neutron stars and black holes. We also find the binary black hole mass distribution has localized over- and underdensities relative to a power-law distribution, with peaks emerging at chirp masses of 8.3−0.5+0.3 and 27.9−1.8+1.9M⊙. While we continue to find that the mass distribution of a binary’s more massive component strongly decreases as a function of primary mass, we observe no evidence of a strongly suppressed merger rate above approximately 60M⊙, which would indicate the presence of a upper mass gap. Observed black hole spins are small, with half of spin magnitudes below χi≈0.25. While the majority of spins are preferentially aligned with the orbital angular momentum, we infer evidence of antialigned spins among the binary population. We observe an increase in spin magnitude for systems with more unequal-mass ratio. We also observe evidence of misalignment of spins relative to the orbital angular momentum.22 MoreReceived 4 February 2022Revised 28 October 2022Accepted 19 December 2022DOI:https://doi.org/10.1103/PhysRevX.13.011048Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.Published by the American Physical SocietyPhysics Subject Headings (PhySH)Research AreasGravitational wave sourcesGravitational wavesTechniquesGravitational wave detectionGravitation, Cosmology & Astrophysics

Numerical relativity higher order gravitational waveforms of eccentric, spinning, nonprecessing binary black hole mergers

We use the open source, community-driven, numerical relativity software, the Einstein Toolkit to study the physics of eccentric, spinning, nonprecessing binary black hole mergers with mass-ratios $q=\{2, 4, 6\}$, individual dimensionless spin parameters $\chi_{1z}=\pm0.6$, $\chi_{2z}=\pm0.3$, that include higher order gravitational wave modes $\ell\leq4$, except for memory modes. Assuming stellar mass binary black hole mergers that may be detectable by the advanced LIGO detectors, we find that including modes up to $\ell=4$ increases the signal-to-noise of compact binaries between $3.5\%$ to $35\%$, compared to signals that only include the $\ell=|m|=2$ mode. We use two waveform models, TEOBResumS and SEOBNRE, which incorporate spin and eccentricity corrections in the waveform dynamics, to quantify the orbital eccentricity of our numerical relativity catalog in a gauge-invariant manner through fitting factor calculations. Our findings indicate that the inclusion of higher order wave modes has a measurable effect in the recovery of moderately and highly eccentric black hole mergers, and thus it is essential to develop waveform models and signal processing tools that accurately describe the physics of these astrophysical sources.

Cover Your Basis: Comprehensive Data-driven Characterization of the Binary Black Hole Population

We introduce the first complete nonparametric model for the astrophysical distribution of the binary black hole (BBH) population. Constructed from basis splines, we use these models to conduct the most comprehensive data-driven investigation of the BBH population to date, simultaneously fitting nonparametric models for the BBH mass ratio, spin magnitude and misalignment, and redshift distributions. With GWTC-3, we report the same features previously recovered with similarly flexible models of the mass distribution, most notably the peaks in merger rates at primary masses of ∼10M ⊙ and ∼35M ⊙. Our model reports a suppressed merger rate at low primary masses and a mass-ratio distribution consistent with a power law. We infer a distribution for primary spin misalignments that peaks away from alignment, supporting conclusions of recent work. We find broad agreement with the previous inferences of the spin magnitude distribution: the majority of BBH spins are small (a < 0.5), the distribution peaks at a ∼ 0.2, and there is mild support for a nonspinning subpopulation, which may be resolved with larger catalogs. With a modulated power law describing the BBH merger rate’s evolution in redshift, we see hints of the rate evolution either flattening or decreasing at z ∼ 0.2–0.5, but the full distribution remains entirely consistent with a monotonically increasing power law. We conclude with a discussion of the astrophysical context of our new findings and how nonparametric methods in gravitational-wave population inference are uniquely poised to complement to the parametric approach as we enter the data-rich era of gravitational-wave astronomy.

Observability of lensing of gravitational waves from massive black hole binaries with LISA

The gravitational waves emitted by massive black hole binaries in the Laser Interferometer Space Antenna (LISA) band can be lensed. Wave-optics effects in the lensed signal are crucial when the Schwarzschild radius of the lens is smaller than the wavelength of the radiation. These frequency-dependent effects can enable us to infer the lens parameters, possibly with a single detection alone. In this work, we assess the observability of wave-optics effects with LISA by performing an information-matrix analysis using analytical solutions for both point-mass and singular isothermal sphere lenses. We use gravitational-waveform models that include the merger, ringdown, higher harmonics, and aligned spins to study how waveform models and source parameters affect the measurement errors in the lens parameters. We find that previous work underestimated the observability of wave-optics effects and that LISA can detect lensed signals with higher impact parameters and lower lens masses.

Dynamical evolution of massive black hole pairs in the presence of spin-dependent radiative feedback

ABSTRACTThe putative ubiquity of massive black holes (MBHs) at the centre of galaxies, and the hierarchical progress of structure formation along the cosmic history, together necessarily imply the existence of a large population of cosmic MBH binaries. Such systems are understood to be the loudest sources of gravitational waves at MHz frequencies, the regime that will be probed by the next Laser Interferometer Space Antenna. It has been proposed that the rate at which MBHs pair and then bind to form binaries is critically dependent upon the feedback exerted by the MBHs on the surrounding gaseous environment. Using the publicly available code gizmo, we perform a suite of simulations aimed at studying the dynamics of an MBH pair embedded in a gaseous disc on ∼100-pc scale. By means of dedicated modules, we follow the dynamics of MBHs in the presence of different spin-dependent radiative feedback models, and compare the results to a benchmark case with no feedback at all. Our main finding is that feedback causes the secondary MBH to shrink its orbit at a reduced pace, when compared with models where feedback is absent. Moreover, such slower inspiral occurs on eccentric orbits, as feedback has the net effect of hampering the circularization process. Though idealized in many aspects, our study highlights and quantities the importance of including spin-dependent feedback recipes in hydrodynamic simulations of MBH pairs, and ultimately in assessing the cosmological coalescence rate of such systems in view of their detection through gravitational waves.

Forecasting constraints on deviations from general relativity in f(Q) gravity with standard sirens

In this work, we explore how modified gravity theories based on the nonmetricity scalar, known as $f(Q)$ gravity, affect the propagation of gravitational waves from inspiraling of binary systems. We discuss forecast constraints on $f(Q)$ gravity by considering standard siren events in two contexts: (i) simulated sources of gravitational waves as black hole--neutron star binary systems, emitting in the frequency band of the third-generation detector represented by the Einstein Telescope (ET); (ii) three standard siren mock catalogs based on the merger of massive black hole binaries that are expected to be observed in the operating frequency band of the Laser Interferometer Space Antenna (LISA). We find that, within the ET sensitivity, in combination with supernova and cosmic chronometer data, it will be possible to test deviations from general relativity at $&lt;3%$ accuracy in the redshift range $0&lt;z&lt;5$, while the main free parameter of the theory is globally constrained at 1.6% accuracy within the same range. In light of LISA's forecasts, combined to supernova and cosmic chronometer data, in the best scenario, we find that the main free parameter of the theory will be constrained at 1.6% accuracy up to high redshifts. Therefore, we conclude that future gravitational wave observations by ET and LISA will provide a unique way to test, with good accuracy, the nature of gravity up to very large cosmic distances.

Eccentric binary black holes: Comparing numerical relativity and small mass-ratio perturbation theory

The modeling of unequal mass binary black hole systems is of high importance to detect and estimate parameters from these systems. Numerical relativity (NR) is well suited to study systems with comparable component masses, ${m}_{1}\ensuremath{\sim}{m}_{2}$, whereas small mass ratio (SMR) perturbation theory applies to binaries where $q={m}_{2}/{m}_{1}\ensuremath{\ll}1$. This work investigates the applicability for NR and SMR as a function of mass ratio for eccentric nonspinning binary black holes. We produce 52 NR simulations with mass ratios between $1:10$ and $1:1$ and initial eccentricities up to 0.8. From these we extract quantities like gravitational wave energy and angular momentum fluxes and periastron advance, and assess their accuracy. To facilitate comparison, we develop tools to map between NR and SMR inspiral evolutions of eccentric binary black holes. We derive post-Newtonian accurate relations between different definitions of eccentricity. Based on these analyses, we introduce a new definition of eccentricity based on the (2,2)-mode of the gravitational radiation, which reduces to the Newtonian definition of eccentricity in the Newtonian limit. From the comparison between NR simulations and SMR results, we quantify the unknown next-to-leading order SMR contributions to the gravitational energy and angular momentum fluxes, and periastron advance. We show that in the comparable mass regime these contributions are subdominant and higher order SMR contributions are negligible.

Unveiling the hosts of parsec-scale massive black hole binaries: morphology and electromagnetic signatures

ABSTRACT Parsec-scale massive black hole binaries (MBHBs) are expected to form in hierarchical models of structure formation. Even though different observational strategies have been designed to detect these systems, a theoretical study is a further guide for their search and identification. In this work, we investigate the hosts properties and the electromagnetic signatures of massive black holes gravitationally bound on parsec-scales with primary mass $\rm {\gt }\, 10^7\, M_{\odot }$. For that, we construct a full-sky light-cone by the use of the semi-analytical model L-Galaxies in which physically motivated prescriptions for the formation and evolution of MBHBs have been included. Our predictions show that the large majority of the MBHBs are placed either in spiral galaxies with a classical bulge structure or in elliptical galaxies. Besides, the scaling relations followed by MBHBs are indistinguishable from the ones of single massive black holes. We find that the occupation fraction of parsec-scale MBHBs reaches up to ${\sim }\, 50{{\ \rm per\ cent}}$ in galaxies with $\rm M_{stellar}\, {\gt }\, 10^{11}\, M_{\odot }$ and drops below 10 per cent for $\rm M_{stellar}\, {\lt }\, 10^{11}\, M_{\odot }$. Our model anticipates that the majority of parsec-scale MBHBs are unequal mass systems and lie at $z\, {\sim }\, 0.5$, with ${\sim }\, 20$ objects per $\rm deg^2$ in the sky. However, most of these systems are inactive, and only $\rm {1-0.1}$ objects per $\rm deg^2$ have an electromagnetic counterpart with a bolometric luminosity in excess of 1043 erg s−1. Very luminous phases of parsec-scale MBHBs are more common at $z\, {\gt }\, 1$, but the number of binaries per $\rm deg^2$ is ${\lesssim }\, 0.01$ at $\rm L_{\rm bol}\, {\gt }\, 10^{45}\,\rm erg\,s^{-1}$.

Evolution of massive black hole binaries in collisionally relaxed nuclear star clusters – Impact of mass segregation

ABSTRACT Massive Black Hole (MBH) binaries are considered to be one of the most important sources of Gravitational Waves (GW) that can be detected by GW detectors like LISA. However, there are a lot of uncertainties in the dynamics of MBH binaries in the stages leading up to the GW-emission phase. It has been recently suggested that Nuclear Star Clusters (NSCs) could provide a viable route to overcome the final parsec problem for MBH binaries at the centre of galaxies. NSCs are collisional systems where the dynamics would be altered by the presence of a mass spectrum. In this study, we use a suite of high-resolution N-body simulations with over 1 million particles to understand how collisional relaxation under the presence of a mass spectrum of NSC particles affects the dynamics of the MBH binary under the merger of two NSCs. We consider MBH binaries with different mass ratios and additional non-relaxed models. We find that mass-segregation driven by collisional relaxation can lead to accelerated hardening in lower mass ratio binaries but has the opposite effect in higher mass ratio binaries. Crucially, the relaxed models also demonstrate much lower eccentricities at binary formation and negligible growth during hardening stages leading to longer merger time-scales. The results are robust and highlight the importance of collisional relaxation on changing the dynamics of the binary. Our models are state-of-the-art, use zero softening, and high enough particle numbers to model NSCs realistically.

Optical follow-up of the tick-tock massive black hole binary candidate

ABSTRACT The observation of a population of massive black hole binaries (MBHBs) is key for our complete understanding of galaxy mergers and for the characterization of the expected gravitational waves (GWs) signal. However, MBHBs still remain elusive with only a few candidates proposed to date. Among these, SDSSJ143016.05 + 230344.4 (’tick-tock’ hereafter) is the only candidate with a remarkably well sampled light curve showing a clear reduction of the modulation period and amplitude over three years of observations. This particular feature has been recently claimed to be the signature of a MBHB that is about to merge. In this paper, we provide an optical follow-up of the tick-tock source using the Rapid Eye Mount (REM) telescope. The decreasing luminosity observed in our follow up is hardly explained within the binary scenario. We speculate about an alternative scenario that might explain the observed light curve through relativistic Lense-Thirring precession of an accretion disc around a single massive black hole.

Massive black hole binaries in LISA: Multimessenger prospects and electromagnetic counterparts

In the next decade, the Laser Interferometer Space Antenna (LISA) will detect the coalescence of massive black hole binaries (MBHBs) in the range $[10^4, 10^8] \, \rm M_{\odot}$, up to $z\sim10$. Their gravitational wave (GW) signal is expected to be accompanied by an electromagnetic counterpart (EMcp), generated by the gas accreting on the binary or on the remnant BH. In this work, we present the number and characteristics (such as redshift and mass distribution, apparent magnitudes or fluxes) of EMcps detectable jointly by LISA and some representative EM telescopes. We combine state-of-the-art astrophysical models for the galaxies formation and evolution to build the MBHBs catalogues, with Bayesian tools to estimate the binary sky position uncertainty from the GW signal. Exploiting additional information from the astrophysical models, such as the amount of accreted gas and the BH spins, we evaluate the expected EM emission in the soft X-ray, optical and radio bands. Overall, we predict between 7 and 21 EMcps in 4 yrs of joint observations by LISA and the considered EM facilities, depending on the astrophysical model. We also explore the impact of the hydrogen and dust obscuration of the optical and X-ray emissions, as well as of the collimation of the radio emission: these effects reduce the number to EMcps to 2 or 3, depending on the astrophysical model, again in 4 yrs of observations. Most of the EMcps are characterised by faint EM emission, challenging the observational capabilities of future telescopes. Finally, we also find that systems with multi-modal sky position posterior distributions represent only a minority of cases and do not affect significantly the number of EMcps.

Preferential accretion and circumbinary disc precession in eccentric binary systems

ABSTRACT We present a suite of high-resolution hydrodynamic simulations of binaries immersed in circumbinary accretion discs (CBDs). For the first time, we investigate the preferential accretion rate as a function of both eccentricity eb and mass ratio qb in a densely sampled parameter space, finding that when compared with circular binaries, (i) mass ratios grow more efficiently in binaries on moderately eccentric orbits (0.0 ≲ eb ≲ 0.4), and (ii) high eccentricities (eb ≳ 0.6) suppress mass ratio growth. We suggest that this non-monotonic preferential accretion behaviour may produce an observable shift in the mass ratio distributions of stellar binaries and massive black hole binaries. We further find that the response of a CBD can be divided into three regimes, depending on eccentricity and mass ratio: (i) CBDs around circular binaries always precess freely, whereas CBDs around eccentric binaries either (ii) undergo forced precession or (iii) remain locked at an angle with respect to the binary periapsis. Forced precession in eccentric binaries is associated with strong modulation of individual accretion rates on the precession time-scale, a potentially observable signature in accreting binaries with short orbital periods. We provide CBD locking angles and precession rates as a function of eb and qb for our simulation suite.

Constraint on the mass of graviton with gravitational waves

We consider the effects of the mass of graviton on both the waveform of gravitational waves and the antenna response to gravitational waves. We find the effect on the response function is negligible for small mass. By using the Fisher matrix method, we make parameter estimations with space-based gravitational wave detectors for massive black hole binaries in massive gravity theory. The wavelength of massive graviton is constrained to be $\lambda_g\gtrsim 10^{17}$ m and the mass is constrained to be $m_g\lesssim 10^{-60}$ kg by one year's observation of massive black hole binaries with space-based gravitational wave detectors.