Black Hole Mass Function of Coalescing Neutron Star-Black Hole Binary Systems: The Prospect of Reconstruction with the Gravitational Wave Observations
Shao-Peng Tang, Hao Wang, Yuan-Zhu Wang, Ming-Zhe Han, Yi-Zhong Fan, Da-Ming Wei
DDraft version March 31, 2020
Typeset using L A TEX twocolumn style in AASTeX63
Black Hole Mass Function of Coalescing Neutron Star-Black Hole Binary Systems: The Prospect ofReconstruction with the Gravitational Wave Observations
Shao-Peng Tang,
1, 2
Hao Wang, Yuan-Zhu Wang, Ming-Zhe Han,
1, 2
Yi-Zhong Fan,
1, 2 and Da-Ming Wei
1, 2 Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing, 210033,People’s Republic of China School of Astronomy and Space Science, University of Science and Technology of China, Hefei, Anhui 230026, People’s Republic ofChina Department of Physics and Astronomy, Purdue University, 525 Northwestern Avenue, West Lafayette, IN 47906, USA
ABSTRACTThe discovery of gravitational waves from compact objects coalescence opens a brand-new windowto observe the universe. With more events being detected in the future, statistical examinations wouldbe essential to better understand the underlying astrophysical processes. In this work we investigatethe prospect of measuring the mass function of black holes that are merging with the neutron stars.Applying Bayesian parameter estimation for hundreds of simulated neutron star–black hole (NSBH)mergers, we find that the parameters for most of the injected events can be well recovered. We alsotake a Bayesian hierarchical model to reconstruct the population properties of the masses of blackholes, in the presence of a low mass gap, both the mass gap and power-law index ( α ) of black holemass function can be well measured, thus we can reveal where the α is different for binary black hole(BBH) and NSBH systems. In the absence of a low mass gap, the gravitational wave data as well asthe electromagnetic data can be used to pin down the nature of the merger event and then measure themass of these very light black holes. However, as a result of the misclassification of BBH into NSBH,the measurement of α is more challenging and further dedicated efforts are needed. Keywords:
Gravitational waves; Black holes; Compact objects INTRODUCTIONThe successful detection of a gravitational wave (GW)signal from the merger of a binary black hole (BBH) byAdvanced LIGO (aLIGO; Abbott et al. 2016a)) on 2015September 14 marks the onset of the era of GW as-tronomy, which opens a new window into observing theuniverse. Since then, dozens of GW events have beendetected (Abbott et al. 2019c), including 10 confidentdetections of BBH mergers, a binary neutron star (BNS)merger event GW170817 (Abbott et al. 2017a) with as-sociated gamma-ray burst (Goldstein et al. 2017) andmacronova/kilonova (Abbott et al. 2017b; Pian et al.2017), and candidates with low false alarm rates (FAR)claimed at the LIGO/Virgo O3 public alerts. In a fewyears, aLIGO and Advanced Virgo (AdV) are antici-pated to reach their design sensitivities, therefore manymore GW signals will be detected (Abbott et al. 2018b).Coalescing BNS and neutron star–black hole (NSBH)
Corresponding author: [email protected] (YZF) binaries attract wide attention, because, in addition togiving rise to GWs, these mergers can also produce elec-tromagnetic transients such as short/long-short GRBsand macronovae/kilonovae, as widely investigated in theliterature (e.g., Eichler et al. 1989; Li & Paczy´nski 1998).In the absence of GW observations, the identificationof macronova/kilonova signals in the afterglow of a fewshort/long-short GRBs provides the strongest supportto their compact object merger origin (see Jin et al. 2016,and the references therein). The GW/GRB/macronovaassociation provides a wealth of physical informationabout the source(s) and allows novel tests of fundamen-tal physics (e.g., Sivaram 1999; Del Pozzo et al. 2013;Li et al. 2016; Miller 2016; Wu et al. 2016; Paschalidis2017), as demonstrated in the case of GW170817/GRB170817A/AT2017gfo (e.g., Abbott et al. 2017c; Wanget al. 2017). Moreover, with the increasing sensitivi-ties of the LIGO/Virgo/KAGRA detectors, the numberof events will accumulate considerably in the next fewyears, reliable statistical studies will become possible. a r X i v : . [ a s t r o - ph . H E ] M a r Figure 1.
Dynamical ejecta masses (left panel) and disk masses (right panel) produced by different chirp masses M c and thespins of BHs χ BH . In this work, we focus on the black hole mass func-tion (hereafter BHMF) of the merging NSBH binaries.Though NSBH binary systems have not been identifiedin the Galaxy yet, they are widely believed to exist inthe universe (Abadie et al. 2010) and the NSBH mergermodel for long-short GRB 060614 has been adopted toreproduce the luminous macronova/kilonova signal (Jinet al. 2015; Yang et al. 2015). Due to the current lim-ited samples of stellar mass BHs, the BHMFs are notwell determined, yet. However, previous studies have al-ready identified some possible characteristics of BHMFfrom the observations of Galactic BHs. For example, thelightest black hole measured in X-ray binaries is ∼ M (cid:12) ( ¨Ozel et al. 2010), much heavier than the upper limit ofneutron stars. Such a result leads us to suspect the exis-tence of a mass gap between the lightest black holes andthe heaviest neutron stars. Population synthesis expectsa high mass cutoff on the power-law mass distributionof black holes (Dominik et al. 2015), because massivestars will lose their masses by stellar wind. Thanks toa high merger rate, such characteristics are expected tobe identified in merging BBH systems via gravitationalwave detection, as demonstrated in Kovetz et al. (2017).For merging NSBH binary systems, the BHMF may bemore challenging because of the expected smaller num-ber of events. Nevertheless, an advantage of construct-ing BHMF of merging NSBH binaries is that the smallchirp mass leads to better mass measurement (Cutler &Flanagan 1994) for the same signal-to-noise ratio (S/N).The other advantage is a good prior knowledge of neu-tron star distribution (Kiziltan et al. 2013), which willcompensate the large measurement error of mass ra-tio. One interesting question is whether the BHMFsare different between the merging NSBH and BBH bi- naries. This consideration is mainly motivated by thefact that neutron star distributions are slightly differentin binary neutron stars and neutron star–white dwarfbinaries ( ¨Ozel et al. 2012; Kiziltan et al. 2013). On theother hand, BBHs may have multiple formation chan-nels, such as binary stellar evolution and dynamical cap-ture. While NSBH binaries have more difficulty formingthrough a dynamical process because of the small massof neutron stars. The BHMF for different binary sys-tems could then be different. However, it is beyond thescope of this work to quantify the prospect of identifyingsuch a difference.As for Bayesian parameter estimation with strain dataof the NSBH merger events, the degeneracy betweenthe mass and spin of BH may produce asymmetric er-rors or biases in mass measurements. Recently, Barbi-eri et al. (2019) showed that the electromagnetic (EM)counterparts information of the NSBH merger can helpto break the degeneracies in the GW parameter space,leading to an unbiased estimation of BH mass comparedto the sole GW data analysis (Veitch et al. 2015, Ta-ble IV). With the works of Kawaguchi et al. (2016)and Foucart et al. (2018), it is straightforward to use( M BH , M NS , χ BH , Λ NS ) to deduce the dynamical ejectamasses M dyn and disk masses M disk which are responsi-ble for powering the electromagnetic emission. As shownin Fig.1, the chance of observing an NSBH merger withEM counterparts may be low for M BH > M (cid:12) , dueto the high ejecta mass requiring low mass and highspin for BH. Thus, for generality, we only consider thesole GW data injected in the Advanced LIGO/Virgo de-tectors with design sensitivities (Abbott et al. 2018b).Therefore, we generate the simulated events and make afull Bayesian parameter estimation for the injected data,then apply a Bayesian hierarchical model to evaluate theprospect of characterizing BHMF.Our work is organized as follows: in Section 2 we in-troduce our BHMF model and the process of generatingsimulated events, analysis of single event using Bayesianparameter estimation, and Bayesian hierarchical modelfor constructing population properties of the BH masses.We present the results and discuss the implications inSection 3. Section 4 contains our discussion and sum-mary. METHODSFor a long time, the mass function of stellar massBHs has remained a topic of interest and a few mod-els have been proposed/investigated ( ¨Ozel et al. 2010;Dominik et al. 2015; Kovetz et al. 2017; Abbott etal. 2019b). BBH merger events detected by AdvancedLIGO/Virgo provide us with a powerful tool to mea-sure the mass and spin of the source, which may tracethe formation channels of BBH systems. However, themass function of BH in NSBH systems still remains un-known because no such event has been reliably identi-fied before. As reported in the LIGO/Virgo O3 publicalerts (GraceDB ), there were four NSBH candidates de-tected in the first six month run (S190814bv, S190910d,S190923y, and S190930t). This indicates a reasonablyhigh merger rate of NSBH systems (note that the suc-cessful detection of NSBH events in late O2 or early O3runs of Advanced LIGO/Virgo has been predicted by Liet al. (2017) based on the macronova/kilonova observa-tions/modeling), thus it is possible to statistically revealthe BHMF with an accumulation of merger events in thenext decade. This work aims to investigate the feasibil-ity of reconstructing the BHMF with dozens of NSBHevents. 2.1. Injection Configurations
We use a phenomenological model to characterize theblack hole and neutron star mass distributions and as-sume that they do not evolve with the redshift (given thelimited distance range of the NSBH events detectablefor Advanced LIGO/Virgo, this approximation is likelyreasonable). The black hole population is assumed toobey a power-law distribution as adopted in Abbott etal. (2016b). In addition, current observations in X-raybinaries suggest a cutoff at ∼ M (cid:12) ( ¨Ozel et al. 2010),while the population synthesis predicts cutoffs in bothlow and high mass bands (Dominik et al. 2015). Veryrecently, a massive unseen companion with a mass of3 . +2 . − . M (cid:12) was identified in the binary system 2MASS https://gracedb.ligo.org/superevents/public/O3/ J05215658+4359220 (Thompson et al. 2019), and a fewMassGap candidate events (S190924h, S190930s, andS191216ap) were claimed in GraceDB. Thus, it is worthinvestigating both scenarios, i.e., the absence and pres-ence of the low mass gap. In this work, we take theBHMF as, P ( M BH ) ∝ M − α BH exp( − M BH /M cut ) , for M gap (cid:54) M BH (cid:54) M (cid:12) , (1)where we set the fiducial values to α = 2 . M gap = 5 M (cid:12) (3 M (cid:12) ; i.e., without the low mass gap),and M cut = 60 M (cid:12) following Abbott et al. (2016b)and Kovetz et al. (2017). These parameters Λ = { α, M gap , M cut } are called hyperparameters that we tryto reconstruct in Sec.2.3.With the data of the first and second observing runsof Advanced LIGO/Virgo, Abbott et al. (2019b) havefurther examined the BBH population properties (e.g.,mass and spin distributions) with different phenomeno-logical models. It is found that components of BBHswith large spins aligned with the orbital angular mo-mentum are unlikely, while a low and restricted (LR)distribution of spin is favored. As show in Fig.1, the lackof GRB and kilonova observations for the four NSBHcandidates also indicates that BHs may have a low spinor alternatively a too “large” BH mass. Therefore, weadopt a low ( L ) spin magnitude distribution with prob-ability density function (PDF) p ( a BH ) = 2 · (1 − a BH ),and a restricted ( R ) distribution of spin’s tilt angle withPDF p (cos θ BH ) = 1 (0 < cos θ BH < F ) spin magnitude distribution (uniformlyspanning in range [0 , . • Case A: With MassGap ( M gap = 5 M (cid:12) ), Low ( L )spin magnitude distribution; • Case B: With MassGap ( M gap = 5 M (cid:12) ), Flat ( F )spin magnitude distribution; • Case C: Without MassGap ( M gap = 3 M (cid:12) ), Low( L ) spin magnitude distribution; • Case D: Without MassGap ( M gap = 3 M (cid:12) ), Flat( F ) spin magnitude distribution.As for the neutron star mass function (NSMF), a trun-cated gaussian distribution is adopted, i.e., P ( M NS ) ∝ √ πσ exp (cid:20) − ( M NS − µ ) σ (cid:21) , for M min (cid:54) M NS (cid:54) M max , (2) Table 1.
Distributions for Injecting Signals and Priors of Bayesian Inference Adopted for GW Parameters (cid:126)θ GW Names Parameters Injection Configurations Priors of Parameter InferenceSource frame mass of BH m /M (cid:12) (1+ z ) a BHMF (Eq.(2.1)) Bounded in (3, 100)Source frame mass of NS m /M (cid:12) (1+ z ) NSMF (Eq.(2.1)) Bounded in (1.1, 2.1)Detector frame chirp mass M c /M (cid:12) ( m m ) / ( m + m ) / Uniform [1.5 × (1+z), 9.8 × (1+z)]Mass ratio q m /m Uniform (0.011, 0.7)Spin magnitude of BH a Low ( L )/Flat ( F ) Uniform (0, 0.99)Spin magnitude of NS a (cid:126) L b cos θ Restricted ( R ) Restricted ( R )Tilt angle between the NS’s spin and (cid:126) L θ φ Uniform(0,2 π ) 0Azimuthal position of (cid:126) L φ JL Uniform(0,2 π ) Uniform (0,2 π )Luminosity distance d L / Mpc Uniform comoving-volume MarginalizedInclination angle θ JN Uniform Sine Uniform SineRight ascension R . A . Uniform(0,2 π ) Uniform (0,2 π )Declination Decl . Uniform Cosine Uniform CosineCoalescence phase φ π ) Uniform (0, π )Geocentric GPS time of the merger t c / s 60 MarginalizedTidal deformability of BH Λ Λ NS az is the cosmic redshift calculated with luminosity distance assuming ΛCDM cosmology b (cid:126) L means the orbital angular momentum where µ , σ , M min = 1 . M (cid:12) , and M max = 2 . M (cid:12) are the mean value, standard deviation, lower and up-per bounds of NS masses, respectively. Based on cur-rent observation data (Kiziltan et al. 2013), we choose µ = 1 . M (cid:12) and σ = 0 . M (cid:12) assuming that NSmass distribution in NSBH systems is similar to thatin BNS. Though the minimum/maximum mass of NS isstill uncertain, our choice of M NS ∈ [1 . , . M (cid:12) is rea-sonable (e.g., Suwa et al. 2018; Cromartie et al. 2019;Tang et al. 2020), and has little influence on our sim-ulations due to the narrow distribution of NSMF (theprobability of injecting very low/high NS mass is prettylow). Additionally, neutron stars would spin down dueto the magnetic dipole radiation and lose their angu-lar momentum during the long merging time scale. Forsimplicity, we only consider the nonrotating NS case,which is in agreement with expectations from GalacticBNS spin measurements (Tauris et al. 2017; Zhu et al.2018), and the approximation of fixing the spin of NS tozero when we inject signals has negligible effect on ourstudy. De et al. (2018) showed that the relation betweentidal deformability and mass of NS approximately obeyΛ( M ) ∝ M − in a relevant mass range. In this workwe take Λ . = 330 (e.g., Abbott et al. 2018a; Jiang et al. 2019) and the tidal deformability is injected asΛ NS ( M NS ) = 330 × ( M NS / . M (cid:12) ) − .All parameters ( (cid:126)θ GW ) used to generate GW wave-forms and their corresponding distributions are summa-rized in Table.1, where we take a uniform comoving-volume distribution up to 500Mpc for luminosity dis-tance ( d L ), and a uniform sky distribution for the loca-tion parameters, i.e., right ascension (R . A . ) and dec-lination (decl . ). To inject the simulated signals, aninspiral-only post-Newtonian waveform template namedSpinTaylorT4Fourier is adopted, which is competent forcomponents with arbitrary, precessing spins (Klein etal. 2014; Veitch et al. 2015). Besides, we take the powerspectral density (PSD ) of design sensitivities into ac-count, which is appropriate for NSBH merger in theaLIGO/AdV era (Abbott et al. 2018b). We set a typ-ical condition that the S/N of a single interferometersatisfying S / N > .
0, as the definition of a GW eventbeing “detected.” This approximately translates into anetwork S / N >
12, which is conventionally used as thethreshold for a network GW detector to identify the GW https://dcc.ligo.org/LIGO-P1200087-v42/public signals (Abadie et al. 2010; Kovetz et al. 2017; Thrane& Talbot 2019).2.2. Single Event Analysis
To examine how well the parameters of BHMF can beconstrained, we first perform a Bayesian parameter in-ference for each simulated event, using the
Bilby package(Ashton et al. 2019) and
P yM ultinest sampler (Buch-ner 2016). Based on the Bayes’ Theorem, the posteriorPDF is proportional to the product of the prior PDF p ( (cid:126)θ GW ) and the likelihood L ( d inj | (cid:126)θ GW ) of the injectedsignal d inj given the waveform model described by (cid:126)θ GW ,i.e., p ( (cid:126)θ GW | d inj ) ∝ L ( d inj | (cid:126)θ GW ) p ( (cid:126)θ GW ). If we assumestationary Gaussian noise, then the log-likelihood of sin-gle detector usually takes the function form,log L ( (cid:126)θ GW ) = − (cid:90) f max f min | d inj ( f ) − h ( (cid:126)θ GW , f ) | S n ( f ) df + C, (3)where S n ( f ), d inj ( f ), and h ( (cid:126)θ GW , f ) represent the one-sided PSD of the noise, the injected signal, and thefrequency domain waveform generated using parameter (cid:126)θ GW , respectively. Due to the lack of reliable numeri-cal simulation based waveform template including tidaleffect for NSBH merger, we only consider the frequen-cies bounded in the range of (23 Hz , f ISCO ) to previewthe situation of analyzing the future real data with aninspiral-only template, e.g., SpinTaylorT4Fourier. Fur-thermore, f ISCO is calculated with the following formu-lae (Bardeen et al. 1972; Apte & Hughes 2019), Z =1+(1 − χ z ) / [(1+ χ z ) / +(1 − χ z ) / ] , Z =(3 χ z + Z ) / ,f ISCO = 6 / [3+ Z − sign( χ z ) (cid:112) (3 − Z )(3+ Z +2 Z )] / + χ z m src1 + m src2 , (4) where χ z = a cos θ is the projection of BH spin alongthe direction of orbital angular momentum, and m src1 , m src2 are source frame masses of the components in unitof M (cid:12) . This procedure has little influence on extractingthe mass and spin information from GW signal, becausethe properties, e.g., chirp mass M c , and mass ratio q , arepredominantly determined by inspiral stage (Damour etal. 2012).Though the real data recorded by AdvancedLIGO/Virgo may suffer from glitch and nonstationarynoise, which may produce biased PSD estimation. Somepowerful tools, e.g., BayesLine and BayesWave, havebeen developed to solve this problem (Cornish & Lit-tenberg 2015; Littenberg & Cornish 2015; The LIGOScientific Collaboration et al. 2020). Here, we only con-sider the ideal case by assuming the PSD can be wellestimated, and use the same PSD and waveform tem-plate as injecting signals to infer the GW parameters () / t o t Figure 2.
Dashed line represents the relation between theratio of visible volume to total spacetime volume V (Λ) / V tot and hyperparameter α , while the shaded area marks the 68%confidence region. of each simulated event. The priors of (cid:126)θ GW are listedin Table.1, where we marginalize the likelihood over thephase φ , geocentric time t c , and luminosity distance d L to accelerate Nest sampling (Allen et al. 2012; Lange etal. 2018; Thrane & Talbot 2019). Due to the compo-nent masses m and m being partially degenerate, wesample the chirp mass M c and mass ratio q instead ofthese parameters to improve the convergence rate of thestochastic sampler (Abbott et al. 2019a). Additionally,we request that the component masses are constrained inreasonable ranges (i.e., m src1 ∈ [3 , m src2 ∈ [1 . , . M c and q .2.3. Bayesian Hierarchical Model
Bayesian hierarchical inference (Adams et al. 2012;Thrane & Talbot 2019) allows us to probe the popula-tion properties of an ensemble of events, and has beenwidely used in various fields, e.g., studying the evo-lutionary scenarios of binary stellar (Taylor & Gerosa2018), revealing the origin of BHs from effective spinmeasurements (Fernandez & Profumo 2019), construct-ing mass distribution of galactic BNS (Farrow et al.2019), constraining of the equation of state (EoS) of NS(Hernandez Vivanco et al. 2019), and investigating thejet properties of short gamma-ray bursts (sGRB; Bis-coveanu et al. 2019).Based on the method introduced by Thrane & Tal-bot (2019) and Galaudage et al. (2019), we apply thistechnique to infer hyperparameters Λ = { α, M gap , M cut } with the likelihood L tot ( (cid:126)d, N | Λ) = N (cid:89) i Z ø ( d i ) n i n i (cid:88) k π ( θ ki | Λ) π ( θ ki | ø) , (5) Table 2.
The 5%-95% Confidence Intervals of Some Recovered Parameters M srcc ( M (cid:12) ) Mass ratio q M srcBH ( M (cid:12) ) χ z R . A . (rad) Decl . (rad) Network S/NInferred 2 . +0 . − . . +0 . − . . +1 . − . . +0 . − . − . +0 . − . . +0 . − . . . .Injected 2 . . . . − . . . . +0 . − . . +0 . − . . +1 . − . . +0 . − . . +0 . − . . +0 . − . . . .Injected 2 . . . . . . . srcc / M q % % RA D E C % % m src1 / M z % % q z % % srcc / M q % % RA D E C % % m src1 / M z % % q z % % Figure 3.
The results of single event analysis. Red lines are injected values, and blue dashed lines represent (5%, 50%, and95%) percentiles. where N , n i , Z ø ( d i ), π ( θ ki | Λ), and π ( θ ki | ø) represent thetotal number of events, the size of downsampled poste-rior samples, the Bayesian evidence of each event, thenormalized BHMF, and the prior of the BH’s sourceframe mass, respectively.Through single event analysis described in Sec.2.2, theBayesian evidences are directly obtained by Nest sam-pling, and the samples of source frame masses of BHs( m src1 ) can be transformed from the posterior samples of M c , q , and the reconstructed d L via m src1 = q − / (1 + q ) / M c z ( d L ) , (6)where z is the cosmic redshift calculated with luminositydistance d L assuming ΛCDM cosmology (Planck Collab-oration et al. 2016).However, due to the fact that higher mass mergersare relatively easier to detect than lower mass merg-ers (Fishbach & Holz 2017), we must take the selectioneffects into account. With the works of Abbott et al.(2019b) and Thrane & Talbot (2019), and assuming auniform-in-log prior of rate, we can marginalize over the Poisson-distributed rate to produce the detection prob-ability p det (Λ | N ) as p det (Λ | N ) ∝ (cid:18) V (Λ) V tot (cid:19) N , (7)where V (Λ) means the “visible volume” which can benumerically calculated with injected signals, and V tot refers to the total spacetime volume. The ratio of detec-tion V (Λ) / V tot is mainly determined by hyperparameter α , because the power-law index describes the profile ofthe population properties. We simulate thousands ofevents with the Monte Carlo method, and collect theevents above the threshold (network S / N >
12) of “de-tecting” GW to approximately evaluate the detectionratio (the relation between this ratio and α is presentedin Fig.2). Thus the likelihood Eq.(5) is modified to L tot ( (cid:126)d, N | Λ , det ) = 1 p det (Λ | N ) L tot ( (cid:126)d, N | Λ) . (8)Finally, we take priors of the hyperparameters of BHMFas α ∼ U (0 . , M cut ∼ U (50 , M (cid:12) , and M gap ∼ U (2 . , . M (cid:12) . P D F =2.21 +0.430.40 M c u t / M % % M cut =66.18 +12.5114.38 M M g a p / M % %
60 80 M cut / M % % M gap / M M gap =5.01 +0.080.31 M P D F =2.09 +0.420.42 M c u t / M % % M cut =65.88 +12.5914.35 M M g a p / M % %
60 80 M cut / M % % M gap / M M gap =4.51 +0.200.49 M P D F =2.27 +0.450.39 M c u t / M % % M cut =65.96 +12.5414.38 M M g a p / M % %
60 80 M cut / M % % M gap / M M gap =3.18 +0.160.25 M P D F =2.42 +0.480.43 M c u t / M % % M cut =65.18 +13.4113.46 M M g a p / M % %
60 80 M cut / M % % M gap / M M gap =3.27 +0.190.26 M Figure 4.
Distributions of hyperparameters reconstructed using 50 simulated events. The top left, top right, bottom left, andbottom right panels show the results of Cases A, B, C, and D, respectively. Red lines are our fiducial values, while the dashedblue lines and the confidence intervals represent (5%, 50%, and 95%) percentiles.3.
RESULTSApplying the Bayesian parameter estimation to eachsimulated event, we can obtain the posterior distribu-tions of the intrinsic parameters, e.g., mass and spin ofBH. Some inference results are shown in Fig.3, and the90% confidence intervals are summarized in Table.2. Asfound in Veitch et al. (2015), the slightly biased me-dian values are owing to the degeneracy between themass ratio q and spin of BH χ z , whose posterior dis-tributions present a strong correlation. The degrees ofbiases are dependent on the S/N and the magnitude ofthe BH’s spin χ z , usually low S/N and high χ z can leadto larger mass measurement error. Because mass ratioand spin are high order post-Newtonian (PN) parame-ters that have minor contributions to the gravitationalwaves (Damour et al. 2012; Baiotti 2019), the inferenceof such parameters will heavily rely on the qualities ofGW data. Besides, high S/N can also reduce the un-certainties of chirp mass (∆ M c / M c ∝ M c / S / N; Cut-ler & Flanagan 1994), while the deviation of M srcc is usually caused by biases of estimating luminosity dis-tance d L and inclination angle θ JN . With the upgradeof Advanced LIGO/Virgo detectors, we expect to detectmore and more high S/N events. If the spins of BHs inNSBH systems share the similar properties with BBHs(i.e., Low and Restricted cases), which is beneficial forparameter estimation, then we can reduce the errors orbiases of mass measurements to a certainly low level.Therefore, it is feasible to perform a Bayesian hier-archical inference to investigate the population proper-ties of BH masses. Though the NSBH merger rate isquite uncertain (Abadie et al. 2010; Li et al. 2017), fourNSBH merger candidate events (S190814bv, S190910d,S190923y, and S190930t) have been claimed in publicalerts of the first half-year LIGO/Virgo O3 run. The im-proved sensitivity in the O4 and full sensitivity runs willfurther enhance the detection rate significantly. Then itis reasonable to assume a sample of ∼ −
100 eventsin the next decade. Fig.4 shows the results of the hy-perparameters reconstructed using 50 events (randomlytaken from 200 simulated events). We note that the M c u t / M
50 60 70 80 90 100N3.54.04.55.05.5 M g a p / M M c u t / M
50 60 70 80 90 100N3.03.54.04.55.0 M g a p / M M c u t / M
50 60 70 80 90 100N2.83.03.23.43.63.8 M g a p / M M c u t / M
50 60 70 80 90 100N3.03.23.43.63.8 M g a p / M Figure 5.
Results of population properties of BH masses obtained with different numbers of simulated events. The top left,top right, bottom left, and bottom right panels show the confidence regions of hyperparameters reconstructed using the inferredposteriors of N events for Cases A, B, C, and D, respectively. The edges of the deep and light blue areas are smoothed with thescatter points (best-fit values and the conservative uncertainties), representing the 1 σ and 2 σ regions, respectively. The dashedblack lines are true values of our BHMF models. high mass cutoff cannot be well identified in NSBH bi-naries, due to the very low expected number of eventswith M BH > M cut . While the gap between NS and BHis mainly determined by the event with the smallest BHmass and can be well constrained. The power-law index α also lies in a relatively narrow region compared to thatin BBH systems (Abbott et al. 2019b).To make a robust evaluation of the uncertain-ties, we use the bootstrap method that ran- domly takes N events repeating 200 times to get { M N , M N , . . . , M iN , . . . , M N } . For each subset M iN in-cluding N groups of inferred posteriors of the N simu-lated events, we fit them with the Bayesian hierarchicalmodel using Nest sampling, and obtain the best-fit val-ues of the hyperparameters together with their fit uncer-tainties σ i fit and the statistical uncertainties σ stat (amongthe best-fit values). Then, we choose σ = (cid:112) σ + σ as the conservative uncertainties, where σ fit is the mean M c / M q M NS = 1.4 M NSBHGAPBNSBBH 5 10 15 20 25 30 m inj1 / M M i s i d e n t i f i c a t i o n P r o b a b ili t y m inj2 = 3 M Figure 6.
Left panel: classification of four mass classes (BNS, NSBH, BBH, and MassGap). Right panel: the probability ofidentifying BBH events as NSBH, calculated using the recovered light component masses of a series of simulated BBH mergerevents. The histogram and errorbars represent the ratios of simulated events with inferred masses m , m , m (medianand 1 σ percentiles of posteriors of m ) that are less than 3 M (cid:12) . value of σ i fit . As shown in Fig.5, the statistical uncer-tainties caused by fluctuation have been greatly reduced,ensuring a better robustness on our results. By increas-ing the number of events, the uncertainties ( σ ) gradu-ally reduced. Though, both α and M gap still suffer fromslight biases, we conclude that the mass gap would beidentified in high significance if dozens of events are de-tected.Note that in Fig.4 and Fig.5, the possible contamina-tion of the NSBH sample caused by the “misclassifica-tion” of the BBH events (see Fig.6(a)), due to the uncer-tainty of the measurement of q , has not been taken intoaccount (some NSBH mergers in principle could also bemisidentified as the BBH mergers. However, as long asthe NS masses do follow a narrow distribution shown ineq.(2.1), such a chance is very low and can be ignored).We have carried out some simulations and found outthat if there is a significant mass gap between neutronstars and black holes (i.e., M gap ∼ M (cid:12) ), such a con-tamination can be ignored. However, in the absence ofthe mass gap, the contamination could be serious (seealso Yang et al. 2018) and the inferred α will be biased.In Fig.6(b) we present the misclassification probability P mis , which is the chance to identify the BBH mergerevents with the light component mass m inj2 = 3 M (cid:12) im-properly as the NSBH ones. Therefore, if in the futurethe absence of the low mass gap has been establishedin the BBH merger events, dedicated simulations withreal PSDs are necessary to reliably infer P mis as a func-tion of m . Together with the well measured BHMF ofthe merging BBH systems, the contamination to the ob-served NSBH merger events can be effectively removed. With the “cleaned” sample, the BHMF of the mergingNSBH systems can be reasonably reconstructed. Sucha detailed approach is of course beyond the scope of thecurrent work. Though the measurement of α is morechallenging, the absence of the low mass gap can bestraightforwardly established because the mergers of the ∼ M (cid:12) BHs with the neutron stars usually are able toproduce short GRBs and bright macronovae/kilonovae.Without the energetic neutrino emission from the cen-tral remnant and due to the higher amount of dynamicalejecta, the macronovae/kilonovae of NSBH mergers areexpected to be different from that from BNS mergers(e.g., Hotokezaka et al. 2013; Jin et al. 2015; Kawaguchiet al. 2020). Moreover, the accurately measured chirpmass can help to distinguish between the NSBH andBNS mergers and the corresponding uncertainty is bet-ter constrained for a relatively high q . SUMMARY AND DISCUSSIONIn this work, we carry out simulations of NSBH merg-ers under four configurations of the population proper-ties of BHs’ spins and low mass breaks. In each case,we perform full Bayesian parameter estimations for allof the simulated events, and apply a Bayesian hierar-chical model to reconstruct the parameters of popula-tion properties of BHs’ masses, i.e., the hyperparame-ters Λ = { α, M gap , M cut } . Though there are still bi-ases of the recovered GW parameters in the analysis ofsome simulated events, the BHMF of all the cases arereconstructed with relatively small uncertainties. In thepresence of a low mass gap (i.e., M gap ≈ M (cid:12) ), our re-sults show a promising prospect of well measuring such0a gap and studying the behavior of BHMF in differentbinary systems. Thus, characterizing BHMF in coalesc-ing NSBH systems from GW measurements is feasible,which may shed new light on the formation or evolu-tionary paths of BHs. So far, it is unclear whether theBHMFs are different for the merging NSBH and BBHsystems. Although the qualification of the prospect ofidentifying such a difference is beyond this work, ourmeasurement errors of distribution parameters are rela-tively small (with fewer event numbers) compared withsimilar works on BBHs (Abbott et al. 2016b; Kovetzet al. 2017). If the BHMFs are considerably different,e.g., ∆ α > . α , it would be plausible to characterizesuch a difference. In the absence of a low mass gap (i.e., M gap ∼ M (cid:12) ), the reconstruction of the BHMF of merg-ing NSBH systems is more challenging because of the(substantial) contamination of the BBH merger events.In this case, we need both the well-reconstructed BHMFfor the merging BBH systems and the misclassificationpossibility of the BBH merger events into NSBH, whichis obtainable via Monte Carlo numerical simulations, toreliably measure α . The determination of the lightestBH mass ( ∼ M (cid:12) ), however, is very straightforward.This is because for such light BHs, the mergers withneutron stars will give rise to bright GRBs and in par-ticular macronovae/kilonovae. Together with the gravi-tational wave data and the benefit of a relatively high q ,this electromagnetic information can help us accuratelyinfer the masses of the BHs. The improvements madeby adding more detectors, such as KAGRA and LIGO-India (Abbott et al. 2017e, 2018b), will be investigatedin the further work. Even if the measurement preci-sion of the parameters of a single event may not greatlyincrease, the increase of sensitive volume will lead tomore detection events and then reduce the statistic er-rors, ensuring a more robust construction of BHMF inthe future.Finally, we would like to note some caveats of our re-sults due to some model dependencies and uncertaintiesin the investigation. In the source parameter estimation,we have ignored some measurement errors that wouldappear in the real data analysis. One of them is the de-tector calibration error which creates uncertainties re-garding strain’s scale and phase. Abbott et al. (2016c)reported that such an error would greatly influence thesky localization but has little effect on mass measure-ment. So the exclusion of such an error does not influ-ence our results. We also fix the PSD as a certain curvein the likelihood function, Eq.(3). In reality, PSD willslowly change with time. One needs to obtain the PSDfrom a piece of data that does not contain signals (at the time period near the event), and parameterize the PSDestimation uncertainty in likelihood function (Veitch etal. 2015). In our simulations, such detailed considera-tion is not possible. We do not consider the systematicerror caused by a waveform template, either. The tem-plate adopted in our work (i.e., SpinTaylorT4Fourier) isan inspiral-only waveform without the merger and ring-down phases, but Abbott et al. (2017d) showed thatcompared with numerical simulation waveforms contain-ing the full inspiral-merger-ringdown phases, this wave-form works well on parameter estimation. We thus ex-pect that such template approximation is fairly good.The term of gravitational wave selection effect (i.e., p det (Λ | N )) is an approximate expression but has beenproved effective in real data analysis (Abbott et al.2016c). If we consider the effect of false alarm rate andcalibration error, the only method of evaluating this is toperform Monte Carlo simulation. Another uncertaintyis the detection rate that relies on the binary mass dis-tribution and LIGO/Virgo/KAGRA’s final sensitivity.Considering the high merger rate of BBHs (e.g., Ab-bott et al. 2016d) and binary neutron stars (e.g., Ab-bott et al. 2018a; Jin et al. 2018), there is no motiva-tion to assume a very low merger rate of NSBH. More-over, there are already four NSBH candidates claimed inthe first half-year O3 run of aLIGO/AdV network. Re-cently, GW190425 (The LIGO Scientific Collaborationet al. 2020) is also shown to be consistent with being anNSBH merger (Han et al. 2020). Therefore, a moder-ately large sample of NSBH mergers is expected to beaccumulated in the near future, with which the blackhole mass function can be reasonably reconstructed.ACKNOWLEDGMENTSThis work was supported in part by NSFC un-der grants of No. 11525313 (i.e., Funds for Distin-guished Young Scholars) and No. 11921003, the ChineseAcademy of Sciences via the Strategic Priority ResearchProgram (Grant No. XDB23040000), and the Key Re-search Program of Frontier Sciences (No. QYZDJ-SSW-SYS024). Software:
Bilby (Ashton et al. 2019, version0.5.5, ascl:1901.011, https://git.ligo.org/lscsoft/bilby/),PyCBC (The PyCBC Team 2018, version 1.13.6,ascl:1805.030, https://github.com/gwastro/pycbc), Py-MultiNest (Buchner 2016, version 2.6, ascl:1606.005,https://github.com/JohannesBuchner/PyMultiNest)1REFERENCES
Abadie, J., Abbott, B. P., Abbott, R., et al. 2010, Classicaland Quantum Gravity, 27, 173001Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016,PhRvL, 116, 061102Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016,Physical Review X, 6, 041015Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016,PhRvL, 116, 241102Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016,ApJL, 833, L1Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017,PhRvL, 119, 161101Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017,ApJL, 848, L12Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017,ApJL, 848, L13Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017,PhRvL, 119, 141101Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017,Classical and Quantum Gravity, 34, 044001Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018,PhRvL, 121, 161101Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018,Living Reviews in Relativity, 21, 3Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019,PhRvX, 9, 011001Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019,ApJL, 882, L24Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019,Physical Review X, 9, 031040Adams, M. R., Cornish, N. J. & Littenberg, T. B. 2012,PhRvD, 86, 124032Allen, B., Anderson, W. G., Brady, P. R., Brown, D. A. &Creighton, J. D. E. 2012, PhRvD, 85, 122006Apte, A. & Hughes, S. A. 2019, PhRvD, 100, 084031Ashton, G., H¨ubner, M., Lasky, P. D., et al. 2019, Bilby:Bayesian inference library, ascl:1901.011Baiotti, L. 2019, Progress in Particle and Nuclear Physics,109, 103714Barbieri, C., Salafia, O. S., Perego, A., et al. 2019, A&A,625, A152Bardeen, J. M., Press, W. H. & Teukolsky, S. A. 1972, ApJ,178, 347Biscoveanu, A., Thrane, E. & Vitale, S. 2019, APS AprilMeeting Abstracts 2019, 64, 3, H16.006Buchner, J. 2016, PyMultiNest: Python interface forMultiNest, ascl:1606.005Cornish, N. J. & Littenberg, T. B. 2015, Classical andQuantum Gravity, 32, 135012 Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2019,Nature Astronomy, 3, 439Cutler, C. & Flanagan, ´E. E. 1994, PhRvD, 49, 2658Damour, T., Nagar, A. & Villain, L. 2012, PhRvD, 85,123007De, S., Finstad, D., Lattimer, J. M., et al. 2018, PhRvL,121, 091102Del Pozzo, W., Li, T. G. F., Agathos, M., et al. 2013,PhRvL, 111, 071101Dominik, M., Berti, E., O’Shaughnessy, R., et al. 2015,ApJ, 806, 263Eichler, D., Livio, M., Piran, T., et al. 1989, Nature, 340,126Farrow, N., Zhu, X.-J. & Thrane, E. 2019, ApJ, 876, 18Fernandez, N. & Profumo, S. 2019, JCAP, 2019, 022Fishbach, M. & Holz, D. E. 2017, ApJL, 851, L25Foucart, F., Hinderer, T. & Nissanke, S. 2018, PhRvD, 98,081501Galaudage, S., Talbot, C. & Thrane, E. 2019, arXive-prints, arXiv:1912.09708Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJL, 848,L14Han, M.-Z., Tang, S.-P., Hu, Y.-M., et al. 2020, ApJL, 891,L5Hernandez Vivanco, F., Smith, R., Thrane, E., et al. 2019,PhRvD, 100, 103009Hotokezaka, K., Kyutoku, K., Tanaka, M., et al. 2013,ApJL, 778, L16Jiang, J.-L., Tang, S.-P., Shao, D.-S., et al. 2019, ApJ, 885,39Jin, Z.-P., Li, X., Cano, Z., et al. 2015, ApJL, 811, L22Jin, Z.-P., Hotokezaka, K., Li, X., et al. 2016, NatureCommunications, 7, 12898Jin, Z.-P., Li, X., Wang, H. et al. 2018, ApJ, 857, 128Kawaguchi, K., Kyutoku, K., Shibata, M., et al. 2016, ApJ,825, 52Kawaguchi, K., Shibata, M., & Tanaka, M. 2020, ApJ, 889,171Kiziltan, B., Kottas, A., De Yoreo, M., et al. 2013, ApJ,778, 66Klein, A., Cornish, N. & Yunes, N. 2014, PhRvD, 90,124029Kovetz, E. D., Cholis, I., Breysse, P. C., et al. 2017,PhRvD, 95, 103010Lange, J., O’Shaughnessy, R. & Rizzo, M. 2018, arXive-prints, arXiv:1805.10457Li, L.-X. & Paczy´nski, B. 1998, ApJL, 507, L59Li, X., Hu, Y.-M., Fan, Y.-Z., et al. 2016, ApJ, 827, 75Li, X., Hu, Y.-M., Jin, Z.-P., et al. 2017, ApJL, 844, L222