IMRPhenomXPHM_NSBH Waveform Model

• IMRPhenomXPHM_NSBH is a precessing NSBH waveform model that integrates higher harmonic modes, spin precession via twisting-up, and matter effects to capture tidal disruption physics.
• The model employs a frequency-domain modal decomposition and analytic Euler angle prescriptions to efficiently incorporate mode asymmetry and calibrated merger–ringdown amplitude suppression.
• Extending previous quadrupole-only approaches, XPNSBH combines NRTidalv3 corrections with NR-tuned amplitude fits, enhancing parameter estimation over a broad calibration domain.

Searching arXiv for the specified papers to ground the article in the cited literature. IMRPhenomXPHM_NSBH, denoted XPNSBH, is a fast, frequency-domain inspiral–merger–ringdown waveform model for quasi-circular neutron star–black hole binaries that combines higher-order harmonic modes beyond the dominant quadrupole, spin precession via twisting-up of an aligned-spin baseline, matter effects in both phase and amplitude, and a merger–ringdown amplitude suppression calibrated to NSBH simulations to capture the physics of tidal disruption (Vidal et al., 3 Jun 2026). It extends IMRPhenomXHM_NSBH by importing the IMRPhenomXPHM precession machinery into an NSBH model whose phases include NRTidalv3 matter effects and whose amplitudes inherit NSBH-specific tidal and disruption features.

1. Lineage and modeling scope

XPNSBH was introduced together with IMRPhenomXHM_NSBH and SEOBNRv5HM_ROM_NRTidalv3_NSBH as part of a framework for gravitational-wave signals from quasi-circular, aligned-spin NSBH binaries including higher-order modes beyond the dominant quadrupole; XPNSBH is the spin-precessing extension of the Phenom branch (Vidal et al., 3 Jun 2026). Its default mode content is

$(\ell,|m|)=\{(2,2),(2,1),(3,3),(3,2),(4,4)\},$

and it omits $m=0$ modes in the aligned-spin baseline. The model is intended for systems in which asymmetry-enhanced higher modes, precession, and matter effects can all be relevant.

The immediate precursor in the Phenom NSBH line is PhenomNSBH, a frequency-domain, quadrupole-only model for aligned-spin NSBH systems. PhenomNSBH combines an IMRPhenomD phase baseline with NRTidalv2 tidal phase corrections and a PhenomC-based aligned-spin amplitude extended by disruption-aware NSBH phenomenology, but it does not include precession or higher harmonics and models only $(\ell,|m|)=(2,2)$ (Thompson et al., 2020). Relative to that earlier construction, XPNSBH propagates tidal corrections and NSBH-specific amplitude features into a precessing higher-mode setting. This suggests a shift from a quadrupole-centered NSBH description toward a model class in which mode asymmetry and spin geometry are treated as first-order ingredients rather than perturbative add-ons.

2. Frequency-domain construction and precession prescription

The waveform is written as a modal decomposition in the frequency domain,

$$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$

with Fourier convention

$$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$

For the non-precessing baseline, equatorial symmetry is imposed as

$$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$

These relations define the aligned-spin structure that XPNSBH subsequently rotates into the inertial frame (Vidal et al., 3 Jun 2026).

Precession is introduced through twisting-up. In a co-precessing frame, the signal is approximated by the aligned-spin IMRPhenomXHM_NSBH waveform, and the inertial-frame strain is obtained by a time-dependent rotation with Euler angles $(\alpha,\beta,\gamma)$,

$$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$

or equivalently,

$$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$

The Euler angles are computed in closed form using multiple-scale analysis of the orbit-averaged post-Newtonian precession equations, exploiting the hierarchy between radiation-reaction and precession timescales; a more costly alternative numerically integrates SpinTaylorT4 precession. The default is the analytic-angle prescription for speed. The effective aligned spin is

$$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$

and the effective precession parameter used in validation and sampling is

$m=0$0

Near merger and ringdown, the aligned-spin baseline adjusts post-merger frequencies and damping times through calibrated remnant-property fits, mitigating unphysical phase behavior and improving the ringdown description. No explicit spherical–spheroidal mixing is modeled. Instead, mode-specific amplitude suppression and a tidal phase prescription are used to represent NSBH physics around merger. A common misconception is that twisting-up alone provides a complete precessing NSBH merger model; the actual construction is more specific, because the precessing sector inherits an aligned-spin NSBH baseline already modified by higher-mode tidal phasing and disruption-aware amplitude suppression.

3. Tidal phasing, higher-mode scaling, and equation-of-state dependence

XPNSBH adds NRTidalv3 matter effects to the black-hole-binary phase mode by mode,

$m=0$1

For the dominant mode,

$m=0$2

where $m=0$3 is the leading PN constant, $m=0$4 is the dynamical tidal parameter, and $m=0$5 is a Padé-resummed rational function constructed to reproduce the 7.5PN tidal phase at low frequency and fitted to binary-neutron-star numerical-relativity waveforms. Higher-mode tidal phases are then included through the scaling

$m=0$6

EOS-dependent spin terms are included separately because the NRTidal fits were obtained from non-spinning NR: spin-squared terms up to 3.5PN and leading spin-cubed terms at 3.5PN are scaled mode by mode in the same way (Vidal et al., 3 Jun 2026).

The model contains an explicit robustness prescription near merger. If $m=0$7 develops poles or $m=0$8 becomes non-concave through a spurious inflection, the affected region is replaced by a second-order Taylor extrapolation beginning at

$m=0$9

with continuity enforced in $(\ell,|m|)=(2,2)$0 and its first two derivatives. This is operationally important because tidal-resummed phase prescriptions can otherwise develop unphysical structure precisely where waveform systematics are most consequential.

The neutron-star tidal deformability is defined as

$(\ell,|m|)=(2,2)$1

with $(\ell,|m|)=(2,2)$2 the quadrupolar gravito-electric Love number. Additional EOS-dependent quantities are mapped from $(\ell,|m|)=(2,2)$3 using quasi-universal relations: the magnetic Love number for an irrotational fluid, the octupolar Love number, and the spin-induced quadrupole and octupole moments. In this formulation the EOS enters through a compact set of deformability-dependent effective quantities rather than through an explicit microphysical EOS model. A plausible implication is that XPNSBH is designed to remain practical for parameter estimation while still carrying a richer matter sector than purely black-hole-binary baselines.

4. Amplitude construction and merger–ringdown suppression

The aligned-spin baseline inherited by XPNSBH uses a piecewise inspiral–merger–ringdown amplitude ansatz with smooth blending,

$(\ell,|m|)=(2,2)$4

where $(\ell,|m|)=(2,2)$5. In the inspiral region,

$(\ell,|m|)=(2,2)$6

and $(\ell,|m|)=(2,2)$7 is obtained by stationary-phase approximation applied to time-domain adiabatic tidal amplitude corrections up to high PN orders. Spin-sector PN terms already present in the BBH baseline are subtracted with their BBH values and reinstated in general form using spin-induced multipole quasi-universal relations. In the merger–ringdown region,

$(\ell,|m|)=(2,2)$8

with inverse-power-law suppression

$(\ell,|m|)=(2,2)$9

The coefficients $$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$0 are determined mode by mode by weighted least-squares fits to collocation points, while $$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$1 is fixed per mode by a Nelder–Mead optimization against NR (Vidal et al., 3 Jun 2026).

The smooth transition window is

$$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$2

with

$$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$3

The transition frequencies are

$$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$4

The choice of $$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$5 is coupled to a safeguard requiring

$$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$6

at the transition, so that PN tidal amplitude corrections do not dominate where they begin to diverge.

Calibration is organized around collocation frequencies derived from MECO, tidal disruption, and ringdown scales: $$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$7

$$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$8

and, when $$h(f,\theta,\phi;\sigma)=\sum_{\ell\ge 2}\sum_{m=-\ell}^{\ell} {}_{-2}Y_{\ell m}(\theta,\phi)\,h_{\ell m}(f;\sigma), \qquad h_{\ell m}(f)=A_{\ell m}(f;\sigma)e^{-i\Psi_{\ell m}(f;\sigma)},$$9, an additional

$$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$0

The tidal MECO frequency is modeled as

$$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$1

while the disruption frequency $$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$2 is scaled to higher modes by $$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$3. Ringdown frequencies $$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$4 and $$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$5 come from fits to the remnant black-hole mass and spin, using NSBH-calibrated remnant-property fits but the same functional form as IMRPhenomXHM to preserve the BBH limit.

The amplitude calibration fits ratios rather than raw amplitudes,

$$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$6

with $$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$7 enforced through the fit strategy. The phenomenological ansatz is

$$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$8

with coefficients $$h(f)=\int_{-\infty}^{\infty} h(t)e^{-i2\pi f t}\,dt.$$9 fitted per mode and collocation frequency, monotonicity enforced in $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$0, and residuals weighted to avoid underweighting small $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$1 values. The overall design is to inherit BBH accuracy where matter effects are negligible and suppress the merger–ringdown amplitude where disruption removes black-hole-binary-like post-merger power.

5. Numerical-relativity calibration, validation, and computational behavior

The higher-mode amplitude baseline was tuned to several NR data sets: 162 SACRA simulations with $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$2, $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$3, and piecewise-polytrope EOSs, but with $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$4 only; 25 BAM higher-mode simulations with subdominant modes up to $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$5, hybridized with DALI for calibration and with NRSurTidal for validation; 12 near-equal/equal-mass BAM simulations with $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$6, modes up to $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$7, and eccentricity $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$8, corrected for center-of-mass drift; and SXS BHNS higher-mode simulations including $$h_{\ell m}(f;\sigma)=(-1)^\ell h^{*}_{\ell,-m}(-f;\sigma).$$9 cases (Vidal et al., 3 Jun 2026). Misaligned black-hole or neutron-star spins and neutron-star spins were excluded from calibration but used for validation. For higher modes, SXS waveforms were weighted by a factor of 4 relative to BAM because the SXS data were cleaner and longer.

Against NR hybrids built from NRSurTidal inspiral plus NR merger, the aligned-spin XPNSBH baseline yields significantly lower mismatches and narrower spreads than $(\alpha,\beta,\gamma)$0-only models. The improvements are largest for asymmetric mass ratios, $(\alpha,\beta,\gamma)$1, and diminish near equal mass where higher modes are symmetry-suppressed. For precession, validation against the SXS:BHNS:0010 configuration, a $(\alpha,\beta,\gamma)$2 precessing case with $(\alpha,\beta,\gamma)$3, showed good agreement through merger and ringdown in the time domain; over a grid in inclination and polarization, matches over $(\alpha,\beta,\gamma)$4 Hz had median mismatch $(\alpha,\beta,\gamma)$5 and maximum $(\alpha,\beta,\gamma)$6. Across the broader parameter space, XPNSBH’s aligned-spin baseline and SEOBNRv5HM_ROM_NRTidalv3_NSBH show good mutual agreement, although mismatches increase toward high $(\alpha,\beta,\gamma)$7 and high aligned $(\alpha,\beta,\gamma)$8.

The recommended coverage is tied to the calibration and validation domain: mass ratio $(\alpha,\beta,\gamma)$9 in higher-mode calibration, extensive $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$0 calibration up to $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$1, black-hole spin $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$2 to $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$3 for aligned configurations, and $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$4 as the sampling prior. Accuracy may degrade for very extreme mass ratios $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$5, strong prograde black-hole spins approaching $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$6 without direct calibration, or very strong precession with large neutron-star in-plane spins, because XPNSBH has no neutron-star-spin calibration. Another delicate regime is the vicinity of the tidal disruption boundary, where amplitude transitions are sharp and the model is tuned to suppress merger–ringdown power but requires care in highly disruptive cases.

The model is designed to remain fast. For evaluation over $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$7 Hz with $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$8 Hz and no OpenMP, XNSBH is about $$h(t;\hat n)=\sum_{\ell m m'} D^{\ell}_{m m'}[\alpha(t),\beta(t),\gamma(t)]\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t)\,{}_{-2}Y_{\ell m}(\hat n),$$9 slower than the BBH baseline XHM, while XPNSBH is about $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$0 slower than XPHM because precession routines dominate. Restricting to $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$1, XNSBH becomes about $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$2 and $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$3 faster than IMRPhenomNSBH and SEOBNRv4_ROM_NRTidalv2_NSBH, respectively. This combination of higher modes, precession, and modest overhead is one of the model’s defining practical features.

6. Inference use, empirical comparisons, and open limitations

For practical analyses, the default XPNSBH higher-mode set is $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$4, with starting frequencies around $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$5 Hz, ending frequencies $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$6 Hz, and $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$7 Hz for long signals (Vidal et al., 3 Jun 2026). The model internally handles NRTidalv3 poles and inflections through extrapolation, so no special user intervention is required. In parameter estimation, standard LVK priors are used, with $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$8 fixed for the black hole and $$h_{\ell m}(t)=\sum_{m'}D^\ell_{m m'}(\alpha,\beta,\gamma)\,h^{\mathrm{co\text{-}prec}}_{\ell m'}(t).$$9 sampled for the neutron star. Edge-on systems are a notable use case because higher modes can dominate, although caution is required near equal masses where odd-$$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$0 modes are symmetry-suppressed and template-bank mode content should be checked.

Comparisons with other models clarify XPNSBH’s position. Relative to IMRPhenomXHM_NSBH, XPNSBH adds precession via twisting-up and propagates higher-mode and tidal corrections into the precessing case; it is recommended whenever precession cannot be excluded. Relative to SEOBNRv5HM_ROM_NRTidalv3_NSBH, both models include higher modes and NRTidalv3 phasing and calibrate amplitudes to the same NR data sets in a similar spirit; pairwise mismatches are low across broad parameter ranges, with differences increasing in highly asymmetric, high-spin configurations where the underlying BBH baselines diverge most. Relative to time-domain GIOTTO and DALI, GIOTTO includes higher modes and precession but does not individually fit tidal higher-mode corrections, which can produce occasional amplitude spikes and larger mismatches in higher-mode-rich regions, while DALI partially mitigates spikes but exhibited pathological behavior, including NaNs and excess power, in some configurations across wide parameter-estimation samples.

Application to real NSBH candidates from GWTC-3 and GWTC-4 yielded posteriors that were largely consistent across XNSBH, XPNSBH, v5HMROM_NSBH, and baseline models at signal-to-noise ratios around $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$1. Higher modes and tides nonetheless produced noticeable shifts in some cases. For GW230518, Bayes factors favored higher-mode NSBH models over IMRPhenomNSBH with $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$2, and the inferred chirp mass and distance shifted when higher modes were included. For GW200105, tides broadened and shifted the posterior on $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$3 and $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$4, and XNSBH/v5HMROM_NSBH favored somewhat larger $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$5 medians than IMRPhenomNSBH. For GW230529, small chirp-mass shifts appeared between BBH and NSBH models, while the NSBH models remained mutually consistent. The overall result was no conclusive evidence for precession or strong tidal imprints at current SNRs, but higher-mode inclusion improved consistency and reduced potential biases.

Several limitations are explicit. XPNSBH does not model orbital eccentricity, non-quasi-circular transients, large neutron-star spin magnitudes, finite-temperature EOS effects, or explicit spherical–spheroidal mixing of quasinormal modes. Planned or needed improvements include extending calibration to more extreme $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$6 and $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$7, incorporating eccentricity, refining higher-mode tidal modeling and disruption-boundary handling, and exploring double-spin precession refinements together with additional higher-mode content. A recurrent misconception is that support for precession implies broad precessing NSBH calibration; the current model supports precession through the XPHM twisting-up construction, but its direct precessing validation is limited to a single-spin SXS case at $$\chi_{\rm eff}=\frac{M_1\chi_{1z}+M_2\chi_{2z}}{M},$$8, and misaligned spins were excluded from calibration.