An aligned-spin neutron-star--black-hole waveform model based on the effective-one-body approach and numerical-relativity simulations

After the discovery of gravitational waves from binary black holes (BBHs) and binary neutron stars (BNSs) with the LIGO and Virgo detectors, neutron-star--black-holes (NSBHs) are the natural next class of binary systems to be observed. In this work, we develop a waveform model for aligned-spin neutron-star--black-holes (NSBHs) combining a binary black-hole baseline waveform (available in the effective-one-body approach) with a phenomenological description of tidal effects (extracted from numerical-relativity simulations), and correcting the amplitude during the late inspiral, merger and ringdown to account for the NS's tidal disruption. In particular, we calibrate the amplitude corrections using NSBH waveforms obtained with the numerical-relativity spectral Einstein code (SpEC) and the SACRA code. Based on the simulations used, and on checking that sensible waveforms are produced, we recommend our model to be employed with NS's mass in the range $1-3 M_\odot$, tidal deformability $0\mbox{-}5000$, and (dimensionless) BH's spin magnitude up to $0.9$. We also validate our model against two new, highly accurate NSBH waveforms with BH's spin 0.9 and mass ratios 3 and 4, characterized by tidal disruption, produced with SpEC, and find very good agreement. Furthermore, we compute the unfaithfulness between waveforms from NSBH, BBH, and BNS systems, finding that it will be challenging for the advanced LIGO-Virgo--detector network at design sensitivity to distinguish different source classes. We perform a Bayesian parameter-estimation analysis on a synthetic numerical-relativity signal in zero noise to study parameter biases. Finally, we reanalyze GW170817, with the hypothesis that it is a NSBH. We do not find evidence to distinguish the BNS and NSBH hypotheses, however the posterior for the mass ratio is shifted to less equal masses under the NSBH hypothesis.


I. INTRODUCTION
In their first two observing runs (O1 and O2), Advanced LIGO [1] and Advanced Virgo [2] have observed gravitational waves (GWs) from ten binary black holes (BBHs) and one binary neutron star (BNS), GW170817 [3]. Recently, in the third observing run (O3), a second BNS, GW190425, was discovered [4]. Other groups have reported additional GW observations analyzing the public data from the first two runs [5][6][7]. Neutron-starblack-holes (NSBHs) may be the next source class to be discovered. Given the lack of a detection in O1 and O2, the rate of NSBHs is uncertain. However, based on estimates from Ref. [8], the expected number of NSBH detections ranges from 0 to 19 in O3 and from 1 to 92 in O4 [9]. As of this writing, in O3, the LIGO and Virgo Collabora-tions have published seven circulars via the Gamma-ray Coordinates Network (GCN) describing detection candidates for which the probability of the system being a NSBH is larger than 1%, and for which the candidate has not been retracted [10][11][12][13][14][15][16]. Furthermore, GW data alone does not exclude the possibility that GW170817 is a NSBH [17][18][19], and it has also been suggested that GW190425 could be a NSBH [20,21]. Therefore it is timely to develop methods that can be used to study NSBHs in GW data.
NSBH binaries exhibit a rich phenomenology that is imprinted on the gravitational waveform (for a review see Ref. [22]). First, as is the case for BNS systems, finite size effects cause a dephasing of the waveform relative to a BBH with the same masses and spins [23][24][25]. Additionally, the amplitude of NSBH waveforms can be  I: Dictionary relating the names we use in this paper for several waveforms from the SEOBNR family and the corresponding names of the waveforms implemented in LAL. The second two waveforms use tidal effects within the NR-Tidal approach.
affected by tidal forces. For unequal mass ratios and slowly spinning BHs, the amplitude of the waveform is well-described by a BBH [26]. On the other hand, for near-equal mass ratios or for highly spinning BHs, depending on the NS equation of state (EOS), the NS can undergo tidal disruption, in which the star is ripped apart as it approaches the BH [27][28][29][30][31]. If the disruption takes place before the NS crosses the innermost stable circular orbit, then the material ejected from the NS can form a disk around the BH [32,33]. If so, starting at a characteristic (cutoff) frequency [34,35], the amplitude of the waveform is strongly suppressed, and the ringdown stage is reduced or even effaced. The details of this process contain information about the NS EOS. Additionally, the disk around the remnant BH and dynamical ejecta can provide the engine for the kilonova signal [36,37], like the one observed for GW170817 [3,38].
In order to take advantage of this potentially rich source of information, it is crucial to have a fast and accurate waveform model capturing effects due to relativistic matter, which can be used in analyzing GW data. Several approaches exist for describing finite-size effects in BNS systems. Tidal corrections [23-25, 39, 40] have been incorporated in the effective one-body (EOB) formalism [41][42][43] in Refs. [44][45][46][47][48][49]. References [50,51] developed a flexible technique that starts from a point-mass BBH baseline waveform, and applies tidal-phase modifications by fitting a Padé-resummed post-Newtonian (PN)-based ansatz to the phasing extracted from numerical-relativity (NR) simulations (henceforth, we refer to this as the NR-Tidal approach). These corrections have been applied to BBH baselines produced within the EOBNR framework [52], and within the inspiral-merger-ringdown phenomenological (IMRPhenom) approach [53,54].
There have been several previous works constructing NSBH waveforms. An aligned-spin NSBH waveform model was developed in Refs. [55,56], but it covered a limited range of mass ratios. In Ref. [57], this waveform model was used in parameter and population studies in conjunction with a former version of the EOBNR BBH baseline [58]. A NSBH model called PhenomNSBH, which was constructed using a similar approach to modeling NSBHs as the one discussed in this paper but developed within the IMRPhenom approach, was recently put forward in Ref. [59]. This model uses the method of Ref. [60] to describe tidal disruption of the amplitude, and uses the tidal phase corrections from Ref. [51].
In this work we develop a frequency-domain model for the dominant, quadrupolar multipole of GWs emitted by aligned-spin NSBH systems. Together with the recent waveform model of Ref. [59], these are the first NSBH models covering a wide range of mass ratios and spin that can be used to analyze GW data. In this paper, we refer to our model as SEOBNR NSBH, which has already been implemented in the LIGO Algorithms Library (LAL) [61]. In Table I we provide a dictionary between the names we use in this work, and the name as implemented in LAL. The amplitude is based on an EOBNR BBH baseline model that we refer to as SEOBNR BBH [52]. We apply corrections inspired by Pannarale et al. [60] to account for tidal disruption. We have adapted the corrections of Ref. [60], originally developed for a former version of the IMRPhenom BBH model [62], for use with EOBNR waveforms [52], augmented with reduced-order modeling (ROM) [63,64] to enhance the speed. Differently from Ref. [59], which uses Pannarale et al. [60]'s fit, here we have performed a fit incorporating results from the new NSBH simulations at our disposal. The phase is computed by applying tidal corrections to the EOBNR BBH baseline [52] using the NRTidal approach, as in SEOBNR BNS [51]. As shown in Ref. [65], even though the tidal corrections from SEOBNR BNS were derived from BNS simulations, they give good agreement with NSBH simulations.
The rest of this paper is organized as follows. In Sec. II, we describe the construction of the waveform model. We summarize the NR waveforms that we use in Sec. II A, review properties of NSBH systems in Sec. II B, give an outline of the waveform model in Sec. II C, summarize the procedure used to calibrate the amplitude correction, assess their accuracy by computing the unfaithfulness, and compare the NSBH waveforms to NR simulations in Sec. II D. Then, we discuss the regime of validity in Sec. II E. Next, we apply the waveform model to several data analysis problems. First, in Sec. III A we estimate when the advanced LIGO-Virgo detector network at design sensitivity can distinguish NSBH and BBH, and NSBH and BNS, systems. Then, in Sec. III B we perform parameter-estimation Bayesian analysis on an NRwaveform, hybridized to an analytical waveform at low frequency, and show the differences between recovering this waveform with a BBH and NSBH model. Finally, we reanalyze GW170817 under the hypothesis that it is a NSBH in Sec. III C. We conclude in Sec. IV by summarizing the main points and lay out directions for future improvements. Finally, in Appendix A, we give explicit expressions defining the waveform model.
We work in units with G = c = 1. Unless otherwise stated, we also assume that the total mass of the binary, M tot = M BH + M NS , where M BH is the mass of the BH and M N S is the mass of the NS, is unity, M tot = 1. The symbol M refers to the mass of the sun.

II. CONSTRUCTING THE NSBH WAVEFORM MODEL
A. Numerical-relativity waveforms In this section we briefly describe the NR data used to construct and validate the model. The Simulating eXtreme Spacetimes (SXS) collaboration has publicly released data from seven simulations described in Refs. [65,66], which were produced using the Spectral Einstein Code (SpEC) [67]. The hyrodynamical part of the code is described in [68,69]. These configurations do not contain spinning BHs, but do include mergers with and without tidal disruption. These simulations use an ideal gas EOS with polytropic index Γ = 2, except for the mass ratio 3 simulation SXS:BHNS:0003, which uses a piecewise polytropic ansatz calibrated to the H1 EOS [70]. We refer the reader to Ref. [65] for further explanation. For five of these simulations the NS spin is zero, and we use these simulations to fit the model as described in Sec. II D. We use the other two simulations for verification. Additionally, SpEC has simulated nine systems with large BH spin in Ref. [71], using the more advanced temperature and composition dependent LS220 EOS [72], which we also use to fit our waveform model. Finally, we validate our NSBH model also against two new SXS waveforms, Q3S9 and Q4S9, which are highly accurate simulations describing disruptive mergers with large BH spin. These simulations were also performed using the Γ = 2 EOS. We give the parameters of all SXS waveforms used here in Table II B. In fitting the model, we also use 134 simulations performed with the SACRA code [73], which were presented in Refs. [28,29].
These simulations span the mass ratios Q = {2, 3, 4, 5}, BH spins χ BH = {−0.5, 0, 0.25, 0.5, 0.75}, and a range of piecewise polytropic EOS. The parameters for all of the waveforms and the EOS used are given Table II of Ref. [56]. Whereas the large number of SACRA waveforms lets us probe a wide parameter range, these waveforms are shorter and of lower accuracy than the publicly available SpEC waveforms as well as Q3S9 and Q4S9, due to finite numerical resolution and non-negligible eccentricity in the initial data.

B. NSBH binary properties
We now turn to a brief description of the final stages of a NSBH coalescence, identifying the main features of the process and the physical properties of the remnant BH. Here, we mainly follow the discussion in Refs. [22,32,60].
The two bodies spiral in due to the loss of energy from the emission of GWs. If the NS approaches close enough to the BH, tidal forces exerted by the BH on the NS can overcome the self gravity of the NS, causing the star to loose mass. This process is called mass shedding. This in turn often leads to tidal disruption, in which the NS is completely torn apart by the strong gravitational field of the BH. Let us denote with r tide the binary's separation at which mass shedding begins. To understand the fate of the NS and the characteristics of the GW signal emitted during the last stages of inspiral, plunge and merger, we compare r tide to the location of the binary's innermost-stable circular orbit (ISCO) (which marks the beginning of the plunge). If r tide < r ISCO , the NS is swallowed by the BH, without loss of material. By contrast, if r tide > r ISCO , mass is ejected from the NS before it plunges. If the NS is far away from the ISCO when it is disrupted, matter may form an accretion disk (torus) around the BH after merger. It has been shown (e.g., see Refs [28,29,32,33] and also below) that the disruption affects the GW signal for NSBH binaries with either nearly equal masses or large BH spins aligned with the orbital angular momentum, because for those systems the condition r tide < r ISCO is satisfied. In Fig. 1, we show an illustration of the effect of tidal disruption on the GW waveform for an example NR hybrid with mass ratio 1.5.   Table II of Ref [56]. We also report the number of GW cycles, NGW, computed up to merger, that up to the peak of the dominant GW mode.
Let us now estimate the radial separation at which mass shedding occurs, r tide , by imposing that the tidal force from the BH balances the self-gravity of the NS. As described in Ref. [32], in the Newtonian limit, r tide can be estimated as r tide ≈ ξ Newt R NS , where ξ Newt = (3Q) 1/3 , R NS is the NS radius, and Q ≡ M BH /M NS is the mass ratio. The factor of 3 is an estimate obtained by matching with NR simulations. This estimate can be improved by accounting for relativistic effects due to the large compactness of the NS, C NS ≡ M NS /R NS [32]. First, r tide is reduced by a factor (1 − 2C NS ) relative to the Newtonian waveforms. In the left panel, we show the waveforms during the long inspiral and mark 20 Hz, which is the lower frequency typically used for LIGO-Virgo parameter-estimation analyses. We also indicate the region toward merger where the NR data are available. In the right panel, we zoom into the last stages of the inspiral and merger. Due to tides, the hybrid waveform and SEOBNR NSBH have a faster inspiral, and end earlier than the corresponding point-mass SEOBNR BBH waveform. Furthermore, because of tidal disruption, the hybrid and SEOBNR NSBH waveforms have no ringdown phase, but end abruptly when the NS gets disrupted. Overall, we find very good agreement between SEOBNR NSBH and NR-hybrid waveforms throughout the entire coalescence. estimate; this factor enforces the absence of tidal disruption in the BH limit C NS → 1/2. Second, point-mass motion in the Kerr metric leads to a correction factor ξ which differs from the Newtonian estimate [32,74]. Combining these effects, we have The relativistic correction parameter ξ is determined by solving the algebraic equation (we take the largest positive root of this equation) where χ BH is the BH's spin. We can associate to the tidal-disruption separation a frequency, which is more useful in the context of modeling the gravitational waveform, as follows which is obtained from the (circular orbit) relation between radial separation and (angular) orbital frequency in the Kerr geometry. 1 1 Note that in [60], the formula for f tide is written in terms of the final, rather than initial, BH mass. In LAL, SEOBNR NSBH is implemented with the final mass. We became aware of this point during the LSC review. The fits in this work were done selfconsistently using the final mass. We have checked that when we replace M f by M BH in the expression for f tide , mismatches with SEOBNR NSBH are O(10 −4 ) or less across parameter space.
The NS compactness C NS , which depends on the NS EOS, enters the expression for r tide . In order to avoid making an assumption about the EOS, it is more convenient to work in terms of the dimensionless tidal-deformability parameter Λ NS , which relates the quadrupole moment of the NS to the tidal field of the companion. The tidal parameter is determined by the compactness of the NS and the tidal Love number k 2 as follows We can relate Λ NS and C NS in an equation-of-state independent way with the Λ NS -C relation [75] with a 0 = 0.360, a 1 = −0.0355, a 2 = 0.000705. In order to achieve continuity with BBH waveforms in the limit Λ NS → 0, for Λ NS ≤ 1, we replace the Λ NS − C relation with a cubic polynomial which interpolates from Λ NS = 1 to C NS = 1/2 at Λ NS = 0, and it is continuous and once differentiable at Λ NS = 1. The matter ejected from the NS, during tidal disruption, can remain bound, forming a disk (torus) around the remnant BH. The mass of this remnant torus, M b,torus , can be determined in terms of the baryonic mass of the NS using fits from Ref. [32] (see also more recent simulations performed in Ref. [76]) where the ISCO radius (r ISCO ) in the Kerr spacetime is given by where the ∓ sign holds for prograde (retrograde) orbits. As mentioned above, the onset of mass shedding occurs when the objects approach within a distance r tide before the NS cross the ISCO. However, the ISCO does not introduce a definite feature in the gravitational waveform. In order to identify the onset of tidal disruption with a definite feature in an NR waveform, in our model we compare f tide to the ringdown frequency of the final BH, f RD , which is the frequency of least-damped quasinormal mode of the final BH. The ringdown frequency can be computed from the final mass and spin using fitting formulas from [77]. To obtain the final mass and spin from the initial parameters of the binary, we use the fits performed by Ref. [78], which account for the ejected mass.

C. Parameterization of the NSBH waveform model
We limit the waveform modeling to the dominant quadrupolar multipole, notably the modes = 2, m = ±2 in the -2 spin-weighted spherical harmonic decomposition of the gravitational polarizations h +,× , and to alignedspin NSBHs. In the frequency domain, we can write the waveform as Henceforth, we focus on the dependence of the amplitude A(f ) and phase φ(f ) on the binary's parameters where we indicate with χ BH and χ NS the (dimensionless) components of the spin aligned with the orbital angular momentum, for the BH and NS, respectively. To compute the GW phase φ(f ), we use the point-mass baseline SEOBNR BBH model, and apply tidal corrections from the NRTidal framework [51]. As shown in Ref. [65] (and as we verify in Figs. 5 and 6), applying NRTidal corrections gives a reasonable approximation of the phase, until the last few cycles.
In order to model the amplitude A(f ), we start with the BNS model SEOBNR BNS as a baseline. Since this model includes tapering beyond the BNS merger frequency [79], we first remove this tapering. This is necessary since the tapering depends on the tidal parameters of both objects, Λ 1 and Λ 2 , and does not vanish as Λ 1 → 0. We note that this means that the Λ 1 → 0 limit of SEOBNR BNS does not correctly describe the amplitude of a NSBH system. We then apply a correction to the amplitude that describes the tidal disruption effects discussed in the previous section. More precisely, we relate the amplitude of SEOBNR NSBH, A(f ), to the amplitude of SEOBNR BNS with no tapering applied, where the correction function w corr is given by and w ± (f ; f 0 , σ) are the hyperbolic-tangent window functions We illustrate the behavior of w corr in Fig. 3. When = 0, w corr (f ) cuts off the amplitude before the end expected for a BBH system with the same masses and spins of the NSBH, and therefore describes tidal disruption. When > 0, the final part of the inspiral and the post-merger signal are still present, but are suppressed relative to the BBH case. The parameters f 0 , σ, , which determine the precise nature of these corrections, are determined by comparing with NR simulations.
Following Ref. [60], we classify the waveforms into four cases: non-disruptive, disruptive, mildly disruptive without torus remnant, and mildly disruptive with torus remnant, depending on the intrinsic parameters of the system. To determine the three parameters {f 0 , σ 0 , } in Eq. (9), we adapt the amplitude model of Ref. [60], which was developed for a different BBH baseline, to the SEOBNR BBH model. We then calibrate the parameters of this model, using the method described in Sec. II D.
Case 1 (Non-disruptive mergers): f RD < f tide , M b,torus = 0. When the tidal frequency is larger than the ringdown frequency of the final BH, the NS reaches the ISCO before crossing r tide . In this case the NS remains intact as it plunges, but with a slightly suppressed amplitude of the ringdown. The waveform is very similar to a BBH. To describe this, we use f 0 = f ND , σ 0 = σ ND , = ND , where ND stands for non-disruptive, and The explicit expressions relating f ND , σ ND , and ND to the intrinsic parameters of the binary are given in Appendix A.
Case 2 (Disruptive mergers): In this case, tidal disruption occurs and a remnant torus of matter forms. For such systems, the typical merger and ringdown stages present for BBHs are exponentially suppressed. To model this case, we set = 0, so that the waveform decays above the frequency f D . This leads to the expression The precise definition is given in Appendix A.
Case 3 (Mildly disruptive mergers with no torus remnant): f RD > f tide , M b,torus = 0. In this case, the NS undergoes mass shedding, but no torus forms around the remnant BH. We combine the information from Case 1 and 2 to determine the cutoff frequency and the width of the tapering. We set Case 4 (Mildly disruptive mergers with torus remnant): f RD < f tide , M b,torus > 0. In this scenario the tidal frequency is above the ringdown frequency, but there is a remnant disk of matter around the BH. As discussed, for example, in Ref. [22], this scenario occurs at large BH spins, and represents the case in which the NS is disrupted before crossing the ISCO, but the size of the tidally disrupted material in the vicinity of the BH is smaller than the BH surface area. Thus, in this case, although a remnant disk eventually forms, the mat-ter does not distribute uniformly around the BH quickly enough to cancel coherently or suppress the BH oscillations. As a consequence, the ending part of the NSBH waveform contains a ringdown signal. In this case, we again combine information from Cases 1 and 2, and fix In Fig. 2, we show the regions of these different parameter spaces, along with relevant NR simulations from the SACRA and SpEC codes.

D. Fitting procedure
The amplitude correction described in the previous section has 20 free parameters, which we denote with the vector λ. The definition of these parameters is given in Appendix A. We fix the coefficients in λ by requiring that the SEOBNR NSBH waveforms agree, as much as pos- sible, with the SpEC and SACRA waveforms described in Sec. II A. For a given NR waveform indexed by I, let us denote the Fourier-domain amplitude of the dominant mode by A NR I (f ; θ). Given the intrinsic binary's parameters θ, and a set of fit parameters λ, we compute the following quantity ∆ I ( λ), to estimate the difference between the frequency-domain amplitude of the model A(f ) and of the NR simulation A NR I (f ). We choose the lower bound of the integral f min to be the frequency at which A NR I (f ) falls to 90% of its initial (lowest-frequency) value; this is a low enough frequency to ensure w corr (f min ) ≈ 1 while avoiding possible contamination from eccentricity in the initial data. For the upper frequency, we take a definition inspired by the cutoff frequency given in [35]. First, we define f max to be the frequency at which f 2 A(f ) takes its maximum value. Then we define f cut to be the frequency (larger than f max ) which satisfies This frequency is larger the ringdown frequency for nondisruptive mergers. For disruptive mergers, f cut gives a characteristic frequency at which the frequency domain amplitude has been suppressed. For the error function in Eq. (14), we consider a constant relative error at each frequency given by σ I (f ) = k I A I (f ). We use k I = 1 for the SACRA waveforms, and k I = 0.1 for the SpEC waveforms, to account for the difference in length and accuracy in We minimize ∆ 2 ( λ) with respect to λ using the Nelder-Mead algorithm [80]. We first use the parameter values from Ref. [60] as an initial guess, and minimize the error over the parameters of the non-disruptive (Case 1) and disruptive (Case 2) window functions separately. We then use the results of this fit as an initial guess for a global fit, including all of the available waveforms in S.
The final results of this global fit are used to define the model, and the numerical values are given in Table VII in Appendix A.
We now turn to a quantitative assessment of the model's performance by comparing against NR simulations. We additionally compare with the recently developed PhenomNSBH model of Ref. [59], in order to understand the performance of the two approximants relative to NR and to each other. To this end, we employ the faithfulness function [81], which is commonly used in LIGO and Virgo data analysis to assess the agreement of two waveforms, e.g., the template τ and the signal s. Let us first introduce the inner product (overlap) between   two waveforms a and b [82,83] a|b where a star denotes the complex conjugate and S n (f ) is the one-sided, power spectral density (PSD) of the detector noise. Here, we use the Advanced LIGO design sensitivity PSD [84]. We compute the faithfulness F by maximizing the overlap over the coalescence time t c , the initial phase φ 0τ of the template τ , and setting the phase of the signal φ 0s to zero at merger, while fixing the same binary's parameters θ for the template and the signal, that is We find it convenient to discuss results also in terms of the unfaithfulness, that isF = 1 − F. Henceforth, we consider the NR waveform as the signal, and the SEOBNR NSBH or PhenomNSBH as the template. In Table III we list the unfaithfulness obtained against all the SXS NSBH waveforms at our disposal, for both SEOBNR NSBH and PhenomNSBH. We also specify the lower frequency f low used to compute the match. We see that both models have broadly similar performance. Since the NR waveforms do not cover the entire bandwidth of the detector, we compute the faithfulness also between both NSBH waveform models, and NR hybrids. We construct hybrids with both SEOBNR NSBH and PhenomNSBH, and compare both waveform models to the two hybrids. The four comparisons have two distinct purposes. First, the low frequency part of SEOBNR NSBH and the SEOBNR NSBH hybrid, and the PhenomNSBH and PhenomNSBH hybrid, are identical up to a shift in the time and phase of the waveform, so that the unfaithfulness of SEOBNR NSBH with an SEOBNR NSBH hybrid quantifies the error of the waveform model failing to capture the NR; the same is true of the unfaithfulness between PhenomNSBH and a Phe-nomNSBH hybrid. Second, comparing SEOBNR NSBH with a PhenomNSBH hybrid, and vice versa, includes the error from the NR part of the waveform, and additionally the error of waveform modeling uncertainty.
To construct the hybrids, we follow the hybridization procedure given in Refs. [51,79]. We first align the waveforms by adjusting the time and phase of SEOBNR NSBH to maximize the overlap with the NR waveform, then we apply a Hann window to smoothly transition from the model to the NR waveform. We refer to the initial and final times of the alignment window as t min and t max . These are chosen for each waveform to produce good agreement in the early part of the waveform. We provide the windows used in Table IV. In Fig. 4, we compare the frequency-domain amplitude of the SEOBNR NSBH model and PhenomNSBH against two publicly available non-spinning SXS waveforms which were used to calibrate SEOBNR NSBH. For context, we additionally show the BBH baseline model, SEOBNR BBH. For SXS:BHNS:0001, which is a non-disruptive merger, the amplitudes of the NR data, SEOBNR BBH, SEOBNR NSBH, and PhenomNSBH, agree well. For the disruptive merger SXS:BHNS:0002, SEOBNR NSBH and PhenomNSBH capture the tapering of the amplitude due to tidal disruption. In Fig. 5  , the NSBH waveform models SEOBNR NSBH and PhenomNSBH, and the BBH model SEOBNR BBH, that is used as a baseline for SEOBNR NSBH. SEOBNR NSBH is able to capture the effects of tidal disruption on the amplitude, while also reducing to BBH-like waveform for large mass ratios when tidal disruption does not occur. same two NR simulations in the time domain. We include the NR error for those waveforms for which it is available, estimated using the methods described in Refs. [65,85]. In Fig 6, we compare SEOBNR NSBH and PhenomNSBH to the accurate spinning simulations simulations Q3S9 and Q4S9, which we use for validation. We align the waveforms using the same procedure to construct the hybrids. We perform these comparisons using the N = 3 extrapolation order.

E. Regime of validity
In Table V, we provide the regime of validity of our SEOBNR NSBH waveform model, which we justify as follows.
• Mass ratio Q. We take the lower limit for the mass ratio to be 1, given that in our fit we in-  clude NR simulations with these mass ratios. For large enough mass ratios, for any spin and Λ NS , Case 1 becomes active and the model reduces to the SEOBNR BBH waveform model. We have checked that there is always a range of parameter space at large mass ratios where this transition occurs, within the regime of validity of the model. Therefore we inherit the upper limit on Q coming from SEOBNR BBH, which is of 100.
• NS mass M NS . Based on EOS's expectations of the maximum NS's mass, we restrict the NS mass to be less than 3M . We also suggest restricting the NS mass to be larger than 1M , which is consistent with the range that we choose for the tidal parameter Λ NS .
• NS tidal-deformability Λ NS . We have verified that sensible waveforms are generated with Λ NS varying from 0 up to 5000, and on this basis suggest the waveform model can be used in this range.
We have also performed a calibration and comparison against available NR simulations to verify the model accurately describes simulations with tidal disruption, as we have described. However the available NR simulations have a more limited range of Λ NS , depending on the NS mass and equation of state, as seen in Fig. 2. Thus we caution that tidal disruption effects are uncertain for large Λ NS , in particular Λ NS 1000 for 1.4M NSs. Even this restricted range includes the bound Λ NS < 800 for a 1.4M NS, obtained from measurements of GW170817 [86].
• BH spin χ BH . In the fit we include simulations with positive spins as large as 0.9, and negative spins as low as -0.5. Since negative spins tend to make the merger less disruptive (i.e., more BBHlike), in order to obtain a symmetric range we suggest [−0.9, 0.9] as a range for the spin.
• NS spin χ NS . While we do not include simulations with NS's spin in the fit, from PN theory we expect that the main effect of the spin enters via the effective aligned spin parameter χ eff Except for mass ratios close to 1 and small spins, χ eff is dominated by the BH spin. We also see rea- We show two of the publicly available SXS simulations with zero spin which were used for calibration of SEOBNR NSBH: the non-disruptive merger SXS:BHNS:0001, with mass ratio 6, and the disruptive merger SXS:BHNS:0002 with mass ratio 2. Also plotted is phase difference for both NSBH models against the relevant NR simulation. The gray band shows the region used to align the model waveforms and NR. We also show the NR phase error for SXS:BHNS:0002 as a horizontal gray band; for SXS:BHNS:0001, which is an older waveform, the NR error is not available. The NR waveform has been shifted in time so the peak amplitude occurs at t = 0, and that the phase is zero there. We see that, for both waveform families, the agreement and NR is very good at the beginning of the NR waveform, but there is dephasing toward the end.
sonable agreement with simulations when the NS's spin is nonzero, as shown in Table III. We therefore recommend that the NS spin is bounded by the low-spin prior that has been used in the literature [3,4], |χ NS | < 0.05.
Through a thorough study, we have verified that the SEOBNR NSBH waveforms look sensible in the region in which we suggest to use this model.

III. APPLICATIONS
Having constructed the SEOBNR NSBH waveform model and checked that it agrees well with existing NR waveforms, we now apply the model to a few dataanalysis problems. In particular, in Sec. III A, we compute the unfaithfulness of the SEOBNR NSBH model against SEOBNR BBH and SEOBNR BNS models in order to obtain an estimate of the regions of parameter space where the advanced-detector network may be able to distinguish different source classes. In Sec. III B, we perform a Bayesian parameter-estimation analysis in which we inject a synthetic NSBH signal (notably a disruptive NSBH merger) and infer the source's proper-ties and parameter's biases when recovering it with the SEOBNR NSBH model and the SEOBNR BBH model. Finally, in Sec. III C, we reanalyze the LIGO/Virgo event GW170817 under the hypothesis that it is a NSBH binary, instead of a BNS.

A. Distinguishing different source classes
When is it possible to determine whether a given binary system is a BBH, BNS, or NSBH based only on the gravitational waveform? We can address this question with our waveform model by considering how similar a SEOBNR NSBH waveform is to a waveform from another source class.
First, we consider the case of distinguishing the hypotheses that a given signal is a BBH or a NSBH. Suppose the signal is a NSBH with a given set of parameters, θ NSBH . We compute the unfaithfulness between the SEOBNR NSBH and SEOBNR BBH models, with the same masses and spins. In the left panel of Fig.,  context, following Refs. [87,88], we estimate that two waveforms are distinguishable at the 1-σ level when the signal-to-noise ratio (SNR) ρ satisfiesF = (D − 1)/2ρ 2 , where D is the number of intrinsic parameters. In our case, considering the intrinsic parameters of the NSBH {m BH , m NS , χ BH , χ NS , Λ NS }, we have D = 5. Then, an unfaithfulness ofF = 10 −3 , corresponds to an SNR of ρ ≈ 45. However, we emphasize that this criterion is sufficient, but it is not necessary, and also it does not say which parameters are biased and by how much. A more detailed discussion of the use of this criterion can be found in Ref. [89].
We also compute the unfaithfulness between the SEOBNR NSBH and SEOBNR BNS models, with the same masses, spins, and tidal parameters. We find that, for zero spin, the unfaithfulness between NSBH and BNS are always less than 10 −3 when the NS mass is less than 3M . This suggests it will be very difficult to distinguish NSBH and BNS systems on the basis of tidal effects on the waveform alone. However, inference on the component masses provides additional useful information that can help distinguish different source classes.
As said above, computing the unfaithfulness does not allow us to quantify its impact on the inference of the binary's parameters and quantify possible biases. Therefore, in the next section, at least for one particular case, we perform a Bayesian parameter-estimation study and extract those biases, and compare with the distinguisha-bility criterion of Refs. [87,88].

B. Parameter-estimation case study
In this section, because of computational costs, we perform a Bayesian parameter-estimation analysis for one specific NSBH system, and postpone to the future a more comprehensive analysis.
We first create a synthetic NSBH signal consisting of an NR hybrid built by stitching together the SEOBNR NSBH waveform to the SXS:BHNS:0006 waveform, which has Q = 1.5 and both spins equal to zero. We do not add a noise realization (i.e., we work in zero noise) which is equivalent to averaging over different noise realizations [90]. We perform four injections, with SNRs of 25, 50, 75, and 100 in the advanced LIGO-Virgo network. While the masses are not astrophysically motivated, this system is interesting to study because it is disruptive. Further, SXS:BHNS:0006 is the simulation with the largest number of cycles of the publicly available SXS waveforms.
We apply the Markov Chain Monte Carlo (MCMC) sampling algorithm implemented in LALInference [91] to these four signals, and recover the signal with both the SEOBNR BBH and SEOBNR NSBH waveform models. Due to limited computational resources, we run the parameter estimation with a lower cutoff frequency of 30 As discussed in the main text, the unfaithfulness can be used to provide an estimate of the SNR at which data can be used to distinguish between two waveforms. Note that it is easier to distinguish BBH and NSBH systems for smaller Q and larger χBH.
Hz. We use a uniform prior on the detector frame component masses. For SEOBNR NSBH, we impose a constraint that M NS < 3M , |χ NS | < 0.05, and |χ BH | < 0.9 consistent with the range of validity of the model. We take a prior on Λ NS that is uniform between 0 and 5000. For SEOBNR BBH, we do not impose a constraint on the maximum mass but do require that the spins of both objects were less than 0.9. Since the two approximants make different assumptions about the nature of the component objects, in describing the results of the Bayesian analysis, we refer to the masses as m 1 and m 2 rather than M BH and M NS . In Fig. 8, we show posteriors in the m 1 − m 2 plane, as well as the q − χ eff plane, for the SNR=25 and SNR=75 injections. For ease of comparison with other PE results by LIGO-Virgo analyses, we show the posterior in terms of the mass ratio q ≡ Q −1 = m 2 /m 1 . For the SNR=25 injection, we see the posteriors from the two waveforms agree very well and are consistent with the injected value within the 90% credible interval. For larger SNRs, posteriors derived using the two waveforms are in tension, and at large enough SNR, the injected value lies outside of the 90% credible interval of the posterior for each model. For the SNR=50 injection, and for larger SNRs, there is a bias in the masses and χ eff recovered using SEOBNR BBH. In particular, SEOBNR BBH recovers a larger total mass. We interpret the results as follows. In order to match the data with a BBH, one raises the total mass of the system with the chirp mass approximately fixed. This brings the ringdown frequency to lower values, in effect mimicking tidal disruption. The masses and spins recovered by SEOBNR NSBH are consistent at the 90% level with the injected values for the SNR=50 case, but are only marginally consistent for the SNR=75 injection, and for SNR=100, the injected values of the masses and χ eff lie outside 90% credible interval. This bias is due to differences with the NR-hybrid waveform. We show the recovery of the SNR=75 injection, for which the SEOBNR NSBH recovery is marginally consistent with the true parameters, in the right two panels of Fig. 8.
In Fig. 9, we show recovery of the tidal parameter Λ NS obtained using SEOBNR NSBH for the 4 different cases. In all four scenarios, the injected tidal parameter is consistent with the 90% credible interval of the Λ NS posterior, although this is only marginally true for the SNR=100 injection. It is interesting to compare the difference between the recovered and injected values, with what is expected from the indistinguishability criterion discussed in the previous section. The unfaithfulness from 30 Hz between the NR hybrid used and SEOBNR NSBH is 2 × 10 −3 , using the advanced LIGO design sensitivity PSD. From the indistinguishability criterion discussed in the previous section, we would expect to see deviations at the 1 − σ level between the posterior recovered with SEOBNR NSBH and the injected value an SNR of 32. A full Bayesian analysis reveals that this level of bias for the recovery ofΛ only arises at a larger value of the SNR. However as we have emphasized, the criterion is only sufficient, it does not specify which parameters are biased, and it has been shown to be quite conservative [89].
This case study illustrates the importance of having accurate NSBH models that can account for tidal disruption in order to derive correct conclusions about astrophysical parameters. However, we emphasize that these injections are only meant as an example. Larger mass ratios may be less tidally disruptive and have tidal effects on the phase suppressed. Conversely, systems with large BH spin will tend to be more disruptive, which will enhance the differences between the BBH and NSBH waveforms.
C. Inference of GW170817 as a NSBH As a final application, we reanalyze GW170817 [92] under the hypothesis that it is a NSBH (see also Ref. [17,18] for related studies). Indeed, it is interesting to ask whether GW data alone can be used to distinguish the hypotheses that this event is a BNS or a NSBH.
We run the Bayesian inference study with the MCMC code implemented in LALInference, using publicly available data of GW170817 from the GW open science center (GWOSC) [93] (discussion of these data for O1 and O2 are contained in Ref. [94]). We run with both the SEOBNR NSBH model, as well as the SEOBNR BNS model in order to be able to do a fair comparison. As far as we know, this is the first time that the new version of the SEOBNR BNS model has been used to analyze GW170817. We compare our results to those from the In the bottom panel, we show the q − χ eff plane. We show the injected value as a black dot. For SNR=25, both SEOBNR NSBH and SEOBNR BBH recover the injected value within the 90% credible interval. For larger SNRs, the recovery with SEOBNR BBH is biased. The posterior with the BBH waveform is peaked around a larger total mass, than the injected one. This effectively shortens the waveform, mimicking tidal disruption. We show this explicitly for SNR=75. Additionally, we use this software injection test to explore how the difference between SEOBNR NSBH and the NR hybrid affects parameter estimation. We see that, at SNR=75, the injected values of the masses and χ eff are marginally consistent with the SEOBNR NSBH posterior at the 90% level; at larger SNRs we find the 90% credible interval does not include the true value.
runs obtained in the GWTC-1 catalog [3], which used a former version of the SEOBNR BNS model. We use the same priors as the GWTC-1 analysis [3], except where otherwise stated. 2 For SEOBNR BNS we assume a flat prior on Λ 1 and Λ 2 , while for SEOBNR NSBH we assume a flat prior on Λ 2 and fix Λ 1 to zero. The posteriors contain support only in the interior of the prior domain for both waveform models with these priors. The prior on the component mass ranges from 0.5 − 7.7M , and therefore the prior does not require that both objects 2 Our analysis of SEOBNR BNS used a prior on Λ 1 ranging from 0 to 3000 and Λ 2 from 0 to 5000. To make a consistent comparison, we remove posterior samples from the GWTC-1 results with Λ 1 > 3000 and reweigh the resulting posterior with the same prior we used for SEOBNR BNS. Including or not including samples with Λ 1 > 3000 does not have a visible effect on the posteriors shown in Fig. 10 or on the results in Table VI. have masses below the maximum mass of a NS. First, we obtain that the median-recovered matchedfilter SNR for each waveform model is 32.7. Since the SEOBNR NSBH and SEOBNR BNS models recover the signal with a similar SNR, we do not find a clear preference either for a NSBH or BNS signal, when we only consider the GW data. Moreover, in Fig. 10, we show the recovery of the mass ratio and tidal deformabilityΛ which is given bỹ In order to more easily compare with results in Ref. [3], we show the mass ratio q ≡ Q −1 = M NS /M BH . There is a preference for unequal mass ratios in the SEOBNR NSBH case, due to tidal disruption that occurs for higher mass ratios. SinceΛ depends non-trivially on the mass ratio and individual tidal parameters, and since Λ 1 is fixed to zero in the prior for SEOBNR NSBH but not for the BNS models, the priors onΛ for the BNS and NSBH models are quite different. In order to make a fair comparison between the posteriors onΛ, we divide each posterior by the prior onΛ, effectively obtaining a flat prior onΛ, as was done in Ref. [3]. We see that under the NSBH hypothesis, smaller values ofΛ are preferred, compared to SEOBNR BNS. Nevertheless, this smaller value ofΛ is consistent with a stiffer equation of state (larger Λ NS ) than in the BNS case, since only one star contributes toΛ in Eq. 20. We give the median and 90% credible intervals for the masses, χ eff ,Λ, and matched filter network SNR, in Table VI.

IV. CONCLUSION
In this work we have built an aligned-spin NSBH waveform model based on the EOB framework, the NRTidal approach and NR simulations: SEOBNR NSBH. This model incorporates a suitable tapering of the frequencydomain waveform's amplitude in regions where tidal disruption occurs, as evinced from physical considerations of tidal effects as the NS plunges into the BH, and from NR simulations of those sources. Tidal corrections to the frequency-domain waveform's phase have been computed using the NRTidal framework. We have shown that SEOBNR NSBH gives good agreement with NR simulations by comparing the waveforms in the frequency domain (Fig. 4) and time domain (Fig. 5), as well as by computing the unfaithfulnesses shown in Tables III and IV. In Fig. 6, we compare the model with two new, highly accurate, simulations from the SXS Collaboration of disruptive NSBHs with highly spinning BHs, which we used for validation. We find very good agreement across a large number of cycles. We also performed the same comparisons with the recently published model PhenomNSBH, and find similar levels of agreement with NR. In Figs. 8 and 9, we have demonstrated that the model can be used to infer properties of NSBH sys- FIG. 10: Reanalysis of GW170817. We show posteriors for the effective tidal deformabilityΛ and the mass ratio q ≡ m2/m1, for three waveform models: in green we show the recovery with an older version of SEOBNR BNS which performed in GWTC-1. In red we show a recovery with the current version of SEOBNR BNS, which has been recalibrated in [51]. The results are broadly consistent, though there are small differences consistent with changes to the waveform model. Finally, in blue we show the recovery with SEOBNR NSBH. The NSBH model has a slight preference for less equal mass ratio and smaller effective tidal parameter, compared to the two BNS models. However, this is consistent with a stiffer equation of state in the NSBH case (larger ΛNS), because there is only one star contributing toΛ. Note that we have reweighed the samples by dividing by the prior onΛ, as done in Ref [3].
tems using software injections, and that at large enough SNR assuming the wrong source class can lead to biased astrophysical inferences. Finally, we have reanalyzed GW170817 with the hypothesis that it is a NSBH instead of a BNS. In Table VI, we see the results are broadly consistent, although there seems to be a slight preference for smaller tidal deformability and unequal masses when recovering with SEOBNR NSBH.
In the future, we plan to extend and improve the SEOBNR NSBH waveform model in various ways. A relatively simple, but important extension is to incorporate information from modes beyond the quadrupolar one using SEOBNRv4HM [95] as a baseline. This is particularly relevant, since the NS can be tidally disrupted also in cases in which the mass ratio is larger than one and the BH's spin is large. Another crucial improvement  is to extend the model to precessing NSBH binaries, since some astrophysical scenarios predict that the BH's spin may be misaligned with the orbital angular momentum.
As more high quality NR simulations of NSBHs become available, it will also be possible to develop a more accurate model for the transition from disruptive to nondisruptive mergers. It will also be interesting to study the effect of using different tidal models in order to quantify uncertainty in the tidal part of the waveform.  VII: Parameters for the amplitude correction wcorr(f ). The parameters above the line appear in the correction for tidally disruptive mergers, the parameters below the line appear in the correction for disruptive mergers.