New and Robust Gravitational-Waveform Model for High-Mass-Ratio Binary Neutron Star Systems with Dynamical Tidal Effects

For the analysis of gravitational-wave signals, fast and accurate gravitational-waveform models are required. These enable us to obtain information on the system properties from compact binary mergers. In this article, we introduce the NRTidalv3 model, which contains a closed-form expression that describes tidal effects, focusing on the description of binary neutron star systems. The model improves upon previous versions by employing a larger set of numerical-relativity data for its calibration, by including high-mass ratio systems covering also a wider range of equations of state. It also takes into account dynamical tidal effects and the known post-Newtonian mass-ratio dependence of individual calibration parameters. We implemented the model in the publicly available LALSuite software library by augmenting different binary black hole waveform models (IMRPhenomD, IMRPhenomX, and SEOBNRv5_ROM). We test the validity of NRTidalv3 by comparing it with numerical-relativity waveforms, as well as other tidal models. Finally, we perform parameter estimation for GW170817 and GW190425 with the new tidal approximant and find overall consistent results with respect to previous studies.


I. INTRODUCTION
Two years after the discovery of a binary black hole (BBH) merger [1], which initiated the era of gravitational-wave (GW) astronomy, the discovery of the binary neutron star (BNS) merger GW170817 inaugurated a new era in multi-messenger astronomy [2][3][4], in which GWs and electromagnetic signals are combined to unravel these highly energetic events.Indeed, the GW detection of GW170817 came along with electromagnetic signatures covering the whole frequency range from radio to gamma-rays [3,4].However, GW170817 was just the beginning.Two years later, the LIGO-Virgo-Kagra collaboration detected the BNS GW190425 [5] and with the increasing sensitivity of current GW detectors as well as planned new-generation detectors, it is expected that the number of BNS events observed will increase in the future [6][7][8][9].
One of the key scientific achievements of GW170817 was the improved knowledge of neutron star (NS) interior [2,10].The structure of the NS is primarily described by an equation of state (EoS) of neutron-rich supranuclear-dense matter [11,12].In fact, matter in the core of NSs is thought to be at least a few times denser than the nuclear saturation density ρ ∼ 10 14 g/cm 3 , and there are various (competing) theoretical and phenomenological models in nuclear physics (e.g., [11,[13][14][15][16][17][18][19]) describing this state of matter.
Given that the EoS describing the NS interior influences the tidal deformability, which is a measure of the star's deformation in response to an external tidal field (e.g., during the gravitational interaction with another compact object [20,21]), BNS merger observations will shed light on the behavior of matter at extreme densi-ties [2,5,20,22,23].
Physical parameters of the system, such as mass, spin, and tidal deformabilities, can be extracted from GW signals by comparing the observed data with theoretical predictions obtained by solving the Einstein field equations (EFEs), usually through Bayesian analysis [24,25].In principle, ab-initio, numerical-relativity (NR) simulations which solve the EFEs on supercomputers would be the method of choice for the description of a BNS waveform.However, such simulations come with high computational costs and typically only cover the late stage of the inspiral and merger [26][27][28][29][30][31][32][33].Furthermore, these simulations are also characterized by an intrinsic uncertainty due to the numerical discretization that is employed to solve the EFEs.
Over the years, the community has developed different waveform models for BNS systems.One class of models is based on post-Newtonian (PN) theory, which expands the EFEs in powers of v/c, where v is the binary characteristic velocity and c is the speed of light.These approximants include tidal interactions that start at the 5th PN order [34] (∝ (v/c) 10 ) up to 7.5PN order [35][36][37], which is the highest PN order that is currently known for tidal effects.However, even high-PN approximants are still not accurate enough in describing the full GW signal, particularly in the late inspiral shortly before the merger, when the velocities of the objects are larger, and the distance between the stars is smaller.One approach for improving the accuracy of PN predictions is made by using the effective-one-body (EOB) formalism [38][39][40][41][42][43].Over the years, two different 'families' of tidal EOB models have emerged: SEOBNRv4T [44,45] and TEOBResumS [46][47][48].Generally, tidal EOB models are more accurate than PN predictions and allow, in most cases, a reliable descrip-tion (within the uncertainty of NR simulations) of the GW signal up to the merger of the stars.To avoid the increase in computational costs, which comes naturally when using EOB (as these models solve ordinary differential equations) instead of PN approximations, people have introduced reduced-order models (ROM) [49][50][51][52], post-adiabatic approximations [53,54], machine learning techniques [55], or a combination of the stationary phase and post-adiabatic approximations [56].
Finally, phenomenological GW models [57,58] also exist, which incorporate EOB and/or NR data with the aim to model the tidal phase contribution during the coalescence as accurately as possible while still being fast and efficient (e.g. through the use of closed-form analytical expressions for the tidal contribution itself) [28,59].One such model is called NRTidal [59], which extracts information from NR, EOB, and PN and constructs an analytical tidal contribution that augments a given BBH waveform baseline, either from EOB models or phenomenological models (such as PhenomD and PhenomX).Two NRTidal versions (NRTidal [59] and NRTidalv2 [60]) have been implemented in the LIGO Algorithm Library (LAL) [61].The second iteration, NRTidalv2, builds upon NRTidal by using improved NR data, adding amplitude corrections to the GW signal, as well as incorporating spin effects [60].Though computationally efficient, both NRTidal models come with certain caveats.First, both versions only include equal-mass systems in the calibration of their fitting parameters.Second, the models consider the tidal bulge Q ij of the star to be directly related to the tidal field E ij generated by its companion via a parameter λ known as the tidal deformability, i.e., Q ij = λE ij .In this case, the former NRTidal models considered only adiabatic tides (λ = const.)throughout the duration of the inspiral [59,60].However, it has been shown that dynamical tides (implying a non-constant, frequency-dependent λ) can arise from these systems due to the quadrupolar fundamental oscillation mode of the stars [20,44,45,62], see also Refs.[63][64][65][66][67][68][69][70][71].
Given these limitations, the aim of this work is to extend the existing NRTidal model.In particular, this new model, called NRTidalv3, includes the following improvements: (i) a larger set of NR BNS waveforms; in total, we utilize 55 waveforms from BAM [31][32][33] and SACRA [28,72] including various total masses, EoSs, mass ratios; (ii) inclusion of mass-ratio dependence of the fitting parameters in the calibrated model; and (iii) dynamical tides [20,44,45,63].
The paper is organized as follows.In Sec.II, we discuss the employed NR data and the hybridization of the NR data with SEOBNRv4T to construct the hybrid waveforms that are used for the calibration.Sec.III introduces the time-domain NRTidalv3 phase, while Sec.IV discusses the frequency-domain phase.We discuss the implementation and tests for the validation of the model in LAL in Sec.V. We then conduct a parameter estimation analysis of the model with existing GW observations in Sec.VI.Finally, our main conclusion and recommendations are discussed in Sec.VII.Throughout the paper, we use geometric units, where G = c = 1, unless otherwise stated.For the individual components of the BNS, we define M A to be the primary mass (mass of the heavier companion), M B to be the mass of the secondary companion, each with tidal deformability Λ A and Λ B , respectively.The total mass is given by M = M A + M B and the mass ratio is defined as q = M A /M B ≥ 1.The aligned spin components are denoted by χ A and χ B .

II. NUMERICAL RELATIVITY DATA AND EXTRACTION OF TIDAL PHASE CONTRIBUTIONS
The NRTidal model is a closed-form GW model that describes the tidal effects in the inspiral part of BNS coalescences [59,60,73].We start from the GW strain where A(t) is the amplitude, and time-domain ϕ(t) is the phase.We assume that we can decompose the phase as where ω = M ω = M dϕ/dt = M (2πf ) is the rescaled GW frequency, ϕ 0 denotes the non-spinning point-particle contribution to the total phase, ϕ SO denotes the spinorbit (SO) coupling, ϕ SS denotes the spin-spin (SS) interactions (both self-spin and spin interactions), ϕ S 3 denotes contributions cubic in spin (S 3 ), and ϕ T is the tidal phase contribution [20].We are neglecting all other (higherorder) spin terms in Eq. (2).Similarly to the time domain, we can also write the frequency-domain strain as with The NRTidal model aims at modeling the tidal contributions ϕ T and ψ T , since, unlike the BBH case, tidal deformabilities are present in BNS and BHNS systems and provide valuable information about the internal composition of the individual stars.

A. Numerical relativity data set
Previous versions of the NRTidal model start with the tidal part of the GW constructed by combining EOB and NR waveforms [59,60].For the calibration of NRTidalv3, 46 NR waveforms simulated by the SACRA code [28,29,72], and 9 waveforms simulated with the  VII).The plot shows overlapping data points for systems with the same masses, but can have different tidal deformabilities due to the various EoSs that were employed.

B. Construction of hybrid waveforms
For each of the NR simulations that we employ for calibration, we compute corresponding configurations employing SEOBNRv4T [50] from the LALSuite [61] library.This means that we compute SEOBNRv4T waveforms including tidal contributions (which we refer to as EOB-BNS later) and without tidal effects, which we call EOB-BBH later.The EOB-BNS waveforms are used in the extraction of the tidal phase contributions, from the early inspiral up to the merger [60].The EOB-BBH and EOB-BNS waveforms are computed at a starting frequency of 15 Hz .These tidal phase contributions will then be used in the calibration of NRTidalv3.
For the construction of hybrid waveforms, we first determine the convergence order of the NR data set.If we find a clear convergence order, we construct a higherorder Richardson-extrapolated waveform assuming this convergence order [74,75].In our set of waveforms, this is the case for the BAM waveforms from Ref. [32,60].For all other data, we are simply employing the highest resolution for the hybrid construction.
In the next step, we align and hybridize the EOB-BNS with the NR waveforms.This is required since our NR data only cover the last few orbits before the merger, but we are planning to construct a closed-form approximant that is valid for the entire frequency range.
For the alignment, we minimize the following integral [26,27,76]: (5) over the chosen frequency interval (or hybridization window, typically near the beginning of the NR waveform) Once the EOB-BNS and the NR waveforms are aligned, we create the hybridized waveform through [26,27,76]: where is the Hann window function.This ensures a smooth transition between the EOB and NR waveforms.The result of this hybridization is shown in Fig. 2 for one example.As indicated above, the final waveforms are labeled EOB-NR.
From the expression of the hybrid waveform above in Eq. ( 6), we expect the tidal phase contributions of the hybrid to be identical to the one of the EOB-BNS waveform right up to the start of the hybridization window.The start of the hybridized waveform is chosen to be ω = 0.0015 (corresponding to f ∈ [17.7, 19.4]Hz).After the window, the tidal contribution would then be identical to that of the NR waveform.For the purposes of constructing NRTidalv3, we only consider the EOB-NR tidal phase contributions up to merger.The post-merger parts of the EOB-NR phase are not included.

C. Extracting the tidal phase
In the next step, we extract the tidal phase contribution ϕ T up to merger by subtracting the EOB-BBH phase from the EOB-NR phase, which (schematically) means: Fig. 3 shows the different tidal phase contributions.We find that the resulting ϕ T calculated from Eq. ( 8) has considerable noise, e.g., due to residual eccentricity or remaining density oscillations during the NR simulations.To avoid this, the noise is filtered using a Savitzsky-Golay filter [32,77] to smoothen the hybrid waveform.This leads to the final tidal phase that we will use as an input for the construction of our NRTidalv3 model.Our construction of the NRTidalv3 approximant begins with the analytical expression of the time-domain tidal phase contribution through 7.5PN order [36,37]: where x = (ω/2) 2/3 , and the c i 's are , where we also have A ↔ B and X A,B = M A,B /M .Note that the coefficients are different from Ref. [60] in NRTidalv2, which employed the PN expression derived in Ref. [35] 1 .The expression employed for NRTidalv3 uses the updated PN expression introduced in Ref [36,37].We also note that the tidal parameters are given by where C A = M A /R A is the compactness of star A in isolation, R A is the radius, and k A 2 is the tidal Love number of the star [21].
For NRTidalv3, we will incorporate dynamical tides, i.e., Λ A,B will be a function of the orbital frequency ω orb = ω/2 of the system [45,66].This stems from the quadrupolar oscillations of the star due to f-mode excitations, which can be represented by a dynamical quadrupole moment obeying the equation of motion of a tidally driven harmonic oscillator plus relativistic corrections such as redshift, frame-dragging, and spin [45,65,66].The dynamical Love number for nonspinning systems can be approximately expressed in terms of an enhancement factor k eff where 1 The PN coefficients in Ref. [35] were corrected in Ref. [36,37], though the actual differences in the computed ϕ PN T values are negligible.where ω 0ℓ is the fundamental mode frequency, Ω′ is the derivative (with respect to the gravitational radiation reaction) of the ratio of the mode and tidal forcing frequencies (which in the nonspinning case is Ω′ = −3/8), where µ = M A M B /M and the quadrupole coefficients are u 2 = 1/4 and v 2 = 3/4.The enhancement factor Eq. ( 11) results from a two-timescale approximation for the dynamical tidal quadrupole with a Newtonian estimate for the orbital evolution [45,66], which should be extended to higher PN orders in future work.
Since k eff 2 is a tidal enhancement factor, k eff 2 → 1.0 at ω orb → 0, which means that the tidal deformability of the NS is at its adiabatic value at large distances (i.e., r → ∞) from its partner.The fundamental frequency (in kHz) M 1.4 f 02 = (M/1.4M⊙ )ω 02 /(2π) rescaled to a 1.4M ⊙ NS can be obtained using the quasi-universal relation found in [78] where This relation is shown in Fig. 4 for various values of Λ, for stars A and B.Then, for ℓ = 2, our tidal parameters for NRTidalv3 contains the modification Note that, as constructed, κ recovers at large distances (very small frequency) the original (constant) adiabatic value (see Fig. 5).

B. Calibrating the time-domain, tidal phase
We then write the effective representation of the timedomain tidal phase in NRTidalv3 as: where we use the following functional form and the same with A ↔ B. The exact functional form employed in this work is based on numerous tests for various different fitting functions and provided overall the best performance.However, it is clearly an assumption, and also other forms could have been used.We also note that a polynomial form of the fitting function could describe the 55 EOB-NR tidal phase hybrids well, but it is not guaranteed to work for extreme cases, i.e., large masses and/or tidal deformabilities.For some of these configurations, the polynomial can become very large for some mass-and tidal-deformability-dependent combination of coefficients at low frequencies.
Taking the Taylor series expansion of Eq. ( 17) and comparing it with Eq. ( 9) allows us to enforce the following constraints to ensure consistency with the PN ex- and similar constraints for A → B. 2 This means that we are left with six unknown parameters, (n ), which can be fitted to the data.However, to make sure that the parameters remain symmetric with respect to stars A and B, we can impose additional constraints on these parameters.We find the following functional form to be sufficient for our model: This leaves us with the 13 parameters a i 's, d 1,i 's, α, and β that we now determine by fitting the function to the EOB-NR hybrid tidal phase data.We note that this parametrization, Eq. ( 20), should ensure that the coefficients of the Padé approximants themselves are functions of the mass-ratio and tidal deformability, thus making itself applicable to a wide range of EoSs.This is unlike the attempt that was made in Ref. [79] where each set of coefficients (n's and d's) of the Padé approximant in Eq. ( 18) was specifically determined for two equations of state describing two SpEC waveforms.The parameters for the time-domain NRtidal phase are shown in Table I.
We show in Fig. 6  The results of the fitting in the time-domain are shown in Fig. 7, together with the fractional differences from the hybridized data.We note the "curving up" of the fits near the merger part of the waveforms; this is primarily due to the dynamical nature of the tides that were incorporated here.For the fits, significant fractional differences ∆ϕ T /ϕ Data T = (ϕ NRT3 T − ϕ Data T )/ϕ Data T are found within the earlier frequencies (where very small phase magnitudes can correspond to large fractional differences) and in the vicinity of the hybridization window, where most of the noise from the NR data can persist even after filtering.

IV. FREQUENCY DOMAIN REPRESENTATION A. PN knowledge and SPA
We now construct the frequency-domain phase using Eq. ( 9) in the stationary phase approximation (SPA) [35,80] The 7.5PN order analytical expression [36,37] of the frequency-domain tidal phase with constant tidal de- formability is then with and similarly for A ↔ B. We use Eq. ( 22) in constraining the frequency-domain phase for NRTidalv3.

B. Dynamical tides in the frequency domain
One of the main differences between NRTidalv2 and NRTidalv3 is the inclusion of dynamical tides in the latter.Therefore, it is not straightforward to define a constant effective κ eff (as was done in NRTidal [59] and NRTidalv2 [60]) that can be used in both the time-and frequency-domain tidal phases, see also Ref. [63].We model the entire frequency-domain dynamical tidal deformability (or Love number) in the same manner as in Eq. ( 11), serving as a frequency-dependent correction/enhancement factor to the adiabatic Love number.This is unlike in Ref. [63], where the frequency-domain dynamical f -mode tidal effects were rewritten as an additive correction to the adiabatic component.Hence, we write the PN expression Eq. ( 22), with κ → κ(ω) as (for ℓ = 2) where keff 2 (ω) is the frequency-domain effective tidal enhancement factor.Note that in this case, the same keff 2 applies to both stars, which makes the calculation of the frequency-domain tidal phase more convenient.We can then solve for keff 2 by substituting Eq. ( 24) into the SPA in Eq. ( 21), yielding the following differential equation: We can then solve for keff 2 by imposing the following initial conditions: for it to behave similarly to the enhancement factor in the time-domain.Solving for keff 2 , however, for different values of the tidal parameter and mass can be computationally inefficient, especially when implemented in lalsuite.For this reason, we introduce a phenomenological representation of the enhancement factor of the frequency-domain Love number keff 2,rep as follows: where the parameters s i , (i = 1, 2, 3) are constrained using s i,j , (j = 0, 1, 2): and κ eff is the effective tidal parameter (29) The fitting parameters s i,j for s i are given in Table II.The fitting formula was chosen to ensure the monotonicity of the enhancement factor keff 2 up to some maximum at least near (or beyond) the merger frequency, and that it follows the initial conditions laid out in Eq. ( 26).The fitting formula satisfies the numerically calculated frequency-domain effective enhancement factor, with a maximum (deviation) error of 1.8% and mean error of 0.2% (see Appendix B) .
TABLE II.Fitting parameters si,j for si that comprise the representation of the frequency-domain Love number enhancement factor, given in Eq. ( 27).The previously derived representation keff 2,rep in Eq. ( 27), which we write from this point onward for convenience as just keff 2 , can finally be used for building the frequency-domain NRTidalv3 phase: where with The remaining parameters are given by pA,B The values of the fitting parameters are found in Table III.For comparison, recently, a phenomenological tidal approximant (with the Padé function up to the x 2 term in the denominator) was also developed in Ref. [42], where some coefficients of the Padé were constrained to the PN expressions for the tidal and self-spin contributions, while the free coefficients were fitted to the TEOBResumS waveform model; here, two sets of these free coefficents were obtained for q = 1 and q > 1.However,  for NRTidalv3, we consider the tides to be dynamical, and the free coefficients are themselves explicit functions of the individual masses and/or tidal deformabilities of the system (as indicated in Eq. ( 34)).
We show the results of the fitting in Fig. 8.In this case, the significant fractional differences with respect to the EOB-NR hybrid data ∆ψ T | = 0, denoted by a black dotted horizontal line.Above this line, NRTidalv2 has greater error with respect to the data than NRTidalv3 and below this line, NRTidalv3 has a greater error.We observe that NRTidalv3 performs overall better than NRTidalv2 for most of the configurations, both for q = 1.0 (yellow curves) and q > 1.0 (blue curves).
In addition, we plot the frequency-domain tidal phases of the different NRTidal models for four waveform hybrids in Fig. 10, representing q = [1.0,1.33, 1.75, 2.0].We note that NRTidalv3 describes the EOB-NR tidal phase contributions better than NRTidalv2 and 7.5PN approximant.We also indicate the merger frequency fit f est from the fitting function shown in Sec.IV G.
The difference between the absolute fractional differences of NRTidalv3 and NRTidalv2 with respect to the hybrid data.The curves for q = 1.0 and q > 1.0 are indicated.We observe NRTidalv3 to perform better (i.e., most curves are above the zero line) than NRTidalv2 for both equal and unequal mass ratios.

D. Spin effects
Spinning NS and BH have an infinite series of nonzero multipole moments [81,82] which will depend on the type of object and its internal structure.Hence, aside from the corrections in the tidal phase contributions that we discussed above, we now include spin effects in NRTidalv3.
Here, we follow the same approach as NRTidalv2 where EoS-dependent spin-squared terms up to 3.5PN were included, as well as leading order spin-cubed terms that enter at 3.5PN order [60,83].But note that the dynamicaltidal enhancement factor Eq. ( 11) is specialized to the nonspinning case in this work.
In terms of the spin-induced quadrupolar C A,B Q and octupolar deformabilities C A,B Oc for stars A and B, the self-spin terms in the phase that are added to the BBH baseline are given by (for aligned spins) where We subtract one to remove the BH multipole contribution that is already present in the baseline BBH phase.Meanwhile, the spincubed term is given by To reduce the number of free parameters, we link and

E. Tidal amplitude corrections
We also include amplitude corrections in the NRTidalv3 frequency-domain model following NRTidalv2 [60].
We incorporate these corrections as an ansatz whose form was taken from the frequencydomain representation of TEOBResumS-NR hybrids, as discussed in Ref. [60]: where D is the luminosity distance to the source.

F. Spin Precession Effects
As with NRTidalv2, we also consider BNS systems whose individual spins have an intrinsic rotation and where this rotation is not necessarily aligned with the orbital angular momentum, causing the orbital plane to precess.
We augment NRTidalv3 onto the BBH baseline IMRPhenomXP, as was done for IMRPhenomXP NRTidalv2 [85], already available in LALSuite.
The IMRPhenomXP baseline improves on IMRPhenomPv2 by incorporating double-spin effects in the twisting-up construction instead of a single-spin approximation.Moreover, the Euler angles are calculated via a precession-averaged treatment of the PN precession dynamics, and there is also an option offered to the user for a more accurate twisting-up prescription based on the solutions of the orbit-averaged SpinTaylorT4 equations implemented in LALSuite.Note that here, the spin-induced multipole corrections introduced in Section IV D are applied in the co-precessing frame and then twisted up.

G. Merger Frequency
As a final ingredient for the construction of the full model, we have to define the stopping criterion or the final frequency until which our model is applica-ble.For this purpose, we use the estimated merger frequency (peak in the GW amplitude) M f est /ν (with ν = X A X B = M A M B /M 2 ) from Ref. [33].The expression is then given by where the factor V M depends on the mass ratio, V S depends on the spin contributions (specifically the alignedspin components of the binary), and V T contains the tidal contributions [33].The factors are where with coefficients Since the Padé approximant is a rational function, asymptotes (corresponding to a zero denominator in the approximant) can appear depending on the specific source parameters.For this reason, we have performed careful checks to verify that no unphysical behavior occurs, i.e., that the estimated merger frequency of the system (as calculated in Eq. ( 41)) is less than the frequencies at which the asymptotes occur (f est < f asymp ).Practically speaking, we generated 30,000 random nonspinning configurations with masses M A , M B ∈ [0.5, 3.0], (corresponding to X A , X B ∈ [0.14, 0.86]) and Λ A , Λ B ∈ [0, 20000].We then performed another 30,000 random spinning configurations with aligned-spin components χ A,B = [−0.5,0.5].For all configurations, we verify that f est < f asymp if f asymp exists.Since NRTidalv3 only models the inspiral part up to the merger of the binary, the condition f est < f asymp suffices for the purposes of the model.
However, to ensure that the phase contribution remains smooth even after the merger, i.e., during the time interval when we taper the waveform, and to remove any asymptote, we make the tidal contribution constant if it reaches a local minimum, then connect the tidal phase from inspiral to merger given by ψ NRT3 T and the postmerger with the PN tidal contribution ψ PN T (which is polynomial in nature and is therefore smooth 3 ) to obtain the tidal function as a function of the frequency: where ψ NRT3a T includes the post-merger information from 3 In principle, any smooth function can do this, but we choose the usual PN for convenience.
PN, and σ(ω) is the Planck taper: for a frequency window after the estimated merger frequency [ω 1 , ω2 ] = [1.15ωest , 1.35ω est ] > ωest = M (2πf est ).The same Planck taper is used to taper the entire waveform abruptly up to 1.2ω est , to minimize the presence of any postmerger signal (since this is not part of the description of NRTidalv3).Then, for the rest of the discussion, we will just refer to the ψ NRT3a T as ψ NRT3 T , knowing that this implicitly includes the tapered phase beyond merger (in its LALSuite implementation).

V. IMPLEMENTATION AND VALIDATION
To make use of the newly developed NRTidalv3 model, we implement it into LALSuite [61] by adding the outlined corrections to several different BBH baseline models: the spinning, non-precessing models IMRPhenomD [86], IMRPhenomXAS [87], SEOBNRv5 ROM [88], and the spinning, precessing model IMRPhenomXP [89].We list all the new approximants in Table IV.The full approximant name is given by combining the BBH baseline name with the NRTidal version (e.g., IMRPhenomXAS NRTidalv3).Versions of NRTidal which are currently present in LALSuite are shown in Table IV, both previously existing models and the ones newly implemented in this work.We also show in Table IV their computation times employing different values of starting frequencies and using a sampling rate of 2 15 = 32768 Hz, for a system M A = M B = 1.35M ⊙ and Λ A = Λ B = 400.Overall, we find that NRTidalv3 has retained both speed and efficiency (particularly for starting frequencies f min ≥ 20 Hz) despite having a more accurate (and mathematically more complicated) description of the GW signal than NRTidalv2.

A. Time Domain Comparisons
In this section, we compare the waveform models utilizing NRTidalv3 with other waveforms, particularly their corresponding NRTidalv2 versions, as well as the SEOBNRv4T model using NR waveforms from BAM [31,33], SACRA [28,29,72], and SpEC [30].
The configurations of the NR waveforms used for comparison are found in Table V  The first column contains the BBH baseline model and the second column contains the versions of NRTidal that were added to the BBH model.Full approximant names are given by joining the BBH baseline name with the NRTidal version.We also include information about the corrections in the approximants (i.e., spin-spin, cubic-in-spin, tidal amplitude, and precession), as well as the computational time ∆T f min of the models at different lower frequencies fmin.The computational time was obtained using the configuration MA,B = 1.35M⊙,ΛA,B = 400, with no spins, and were simulated using an Apple M1 Pro processor.For the time-domain comparison, we are first aligning the waveforms using the same procedure as for the construction of the EOB-NR hybrids described in Sec.II.Fig. 11 shows our comparison for the BAM waveforms, where we also indicate the numerical uncertainty as shaded bands.In particular, we use a blue band for setups for which we just calculate the phase difference be-tween the two highest resolutions 4 .For NR simulations where we find a clear convergence order and employ the Richardson-extrapolated waveform, we are using a green band as an error measure given by the phase difference between the Richardson-extrapolated waveform and the FIG.11.Time-domain dephasing comparisons for the BAM, SACRA, and SpEC waveforms.For each NR waveform, the upper panel shows the real part of the gravitational wave strain as a function of the retarded time, while the bottom panel shows the phase difference between the waveform model and the NR waveform.We note that BAM:0001, BAM:0037, and both SACRA waveforms (indicated by blue frames in their plots) were also used in the calibration of NRTidalv3.highest resolution.

BBH baseline ψT
We observe from Fig. 11 that for the BAM NR data, both NRTidalv3 models perform well in terms of the dephasing with respect to the NR waveforms, and the NRTidalv3 models are consistent with the other models such that they generally fall within the estimated numerical uncertainties.The exception to this would be for BAM:0081 and BAM:0094 (the same can be said for the other waveform models), which are characterized by large mass ratios, and large tidal deformabilities which are far outside the space of calibration for NRTidalv3, though we still observe NRTidalv3 to perform a bit better than NRTidalv2.We also note that the employed EoS for these setups (MS1b) is already disfavored by the observation of GW170817 [2][3][4], i.e., this could be considered as a reasonable upper bound for a realistic BNS setup.We also present the comparison between the models and the SACRA and SpEC waveforms.Generally, due to the higher resolution, the computed uncertainties are generally smaller than for the BAM waveforms, however, we do not see a clear convergence order, which means that we cannot estimate the error based on Richardson extrapolation.For SACRA waveforms, we see that the NRTidalv3 models perform better than the other waveform models.However, we note that this might also simply be due to the fact that the same waveforms were also used in the calibration of NRTidalv3.Meanwhile, no model was able to capture the merger part of SXS:NSNS:0001 5 , while 5

An alternative time-domain dephasing comparison was done with
SEOBNRv4T, but not using a quasi-universal relation to compute the waveform parameters for the polytropic EoS (used for SXS:NSNS:0001), and the result stays the same as in Fig. 11.NRTidalv3 performed better than the others in terms of the dephasing for SXS:NSNS:002.

B. Mismatch against EOB-NR hybrid data
We further test the accuracy of our fits by comparing mismatches of the NRTidalv2 and NRTidalv3 models with the hybrid EOB-NR waveforms.The mismatch (or unfaithfulness) between two frequency-domain complex waveforms h1 and h2 is given by where the overlap is and the maximization of the overlap for some arbitrary phase ϕ c and time shift t c (in the time domain, corresponding to a phase and frequency shift in the frequency domain) ensures the alignment between the two waveforms.Here, we assume the spectral density of the detector as S n (f ) = 1.0 so that the computation is detectoragnostic.We set f max = f mrg .For the waveforms, we choose a sampling rate of 2 13 = 8192 Hz.For a full comparison, we construct full waveforms with corresponding BBH baselines (from the model we want to compare with) and add to its phase the EOB-NR hybrid tidal phase, i.e., we have IMRPhenomD Data, IMRPhenomXAS Data, and SEOBNRv5 ROM Data.The minimum frequency f min is chosen to be either 20 Hz or 350 Hz.The former takes into account almost the entire hybrid waveform, which is dominated in the very early inspiral by the EOB model used for constructing the EOB-NR hybrid (SEOBNRv4T), while the latter value considers mainly the NR contribution (typically where the hybridization starts).
We show results as mismatch histograms for IMRPhenomD NRTidalv3, IMRPhenomXAS NRTidalv3, and SEOBNRv5 ROM NRTidalv3 in Fig. 12. Taking a look, for example, at the histograms for IMRPhenomD NRTidalv3 with f min = 20 Hz (the uppermost left panel), we compute the mismatches between IMRPhenomD NRTidalv3 and IMRPhenomD Data, between IMRPhenomD NRTidalv3 and IMRPhenomD NRTidalv2, and between IMRPhenomD NRTidalv2 and IMRPhenomD Data.
For each histogram, the mean mismatch µ F is indicated by a vertical dashed line.The mismatches are larger with larger f min .We note the improvement in the accuracy of the NRTidalv3 waveforms with respect to NRTidalv2, as indicated by their lower mismatches with respect to the hybrid data.In general, NRTidalv3 has a mismatch of about half of an order of magnitude smaller than that of NRTidalv2.We also note that for f min = 350 Hz, NRTidalv2 has a tail in its distribution for F ≳ O(10 −2 ), while this is not present for NRTidalv3.While the reduction of the mismatch is to some extent also caused by the fact that all of the employed waveforms have been used in the calibration of NRTidalv3, it still shows the overall robustness and performance of the model.

C. Mismatches against other tidal waveform models
Finally, we want to compare our model with other tidal models in three cases: (a) nonspinning configurations, (b) aligned spins, and (c) precessing systems.
For the non-spinning configurations, we compute the mismatch between NRTidalv3 (specifically, IMRPhenomXAS NRTidalv3 and SEOBNRv5 ROM NRTidalv3) and TEOBResumS, SEOBNRv4T, and NRTidalv2.The mismatches computed in this section assume f min = 40 Hz, f max = 2048 Hz and a waveform sampling rate of 2 13 = 8192 Hz, except for when comparisons are done with SEOBNRv4T, where a sampling rate of 2 13 = 4096 Hz was used for the sake of efficiency.
a. Non-spinning setups: In the non-spinning case, the mismatches are computed for 4000 random configurations of masses M A,B = [0.5, 3.0]M ⊙ and tidal deformabilities Λ A,B = [0, 5000], for comparisons done with TEOBResumS and SEOBNRv4T.For Λ A,B ≥ 5000, and especially when combined with a high-mass-ratio or high-spin configuration, we find both TEOBResumS and SEOBNRv4T to fail to produce a waveform.For comparisons which do not involve these two waveform models, we use Λ A,B = [0, 20000].These mismatches are then plotted in a 2D scatter plot (in the x-y plane) with the masses and tidal deformabilites (Fig. 13).In addition, we also directly compare the mismatches between the NRTidalv3 models.We note the very small mismatches, with µ F = 3.608 × 10 −4 , indicating that the waveforms behave very similarly to each other.In general, we find very small mismatches for relatively low masses and low-mass ratios, with the largest mismatches occurring at higher masses (M A,B ≳ 2.0M ⊙ ) and high mass-ratios (q ≳ 1.5).Large mismatches also occur for large tidal deformabilities (Λ A,B ≳ O(10 3 )).For all comparisons in the non-spinning case (except IMRPhenomXAS NRTidalv3 vs SEOBNRv5 ROM NRTidalv3) µ F = O(10 −3 ) and max( F ) = O(10 −2 ) (see Table VI ).
b. Aligned-spin setups: For the aligned-spin configurations, we compute the mismatches between two NRTidalv3 variants, IMRPhenomXAS NRTidalv3 and SEOBNRv5 ROM NRTidalv3 and these NRTidalv3 variants against TEOBResumS and SEOBNRv4T.In addition to the 4000 random configurations of the masses and tidal deformabilities used in the non-spinning case, that we generated for the non-spinning case, we also generate 4000 random aligned spins χ A,B = [−0.5,0.5].In the mismatch scatter plots, we plot the masses, the tidal deformabilities, and the aligned-spin components in the x − y plane, and the mismatches in the color bar (see Fig. 14).When the NRTidalv3 variants are compared against TEOBResumS, large tidal deformabilities (Λ = O(10 3 )), mass-ratios (q ≳ 1.5), and masses (M A,B ≳ 2.0) result to large mismatches.In general, comparisons with SEOBNRv4T yield on average larger mismatches than that with TEOBResumS (see Table VI).The mismatches are largest (O(10 −1 )) for the configurations with very high spin magnitudes |χ A,B | ≳ 0.4.c.Precessing setups: For the precessing configurations, we compare IMRPhenomXP NRTidalv3 and IMRPhenomPv2 NRTidalv2 [60,95].Aside from the randomly generated masses and tidal deformabilities used in the non-spinning and aligned-spin cases, we also generate 4000 random values of the x, y, z components of the spin, χ x,y,z A,B = [−0.5,0.5].Here, for the spins, we plot in Fig. 15 on the x − y plane the effective spin precession χ p against the effective aligned spin χ eff , where where S p is the average magnitude of the spins: with K A = 2 + 3/(2q) (with M A > M B ), and S A,⊥ is the magnitude of the in-plane spin for body A [96].Mean-  16.The marginalized 1D and 2D posterior probability distributions for selected parameters of GW170817, obtained with IMRPhenomXP NRTidalv3 (black), IMRPhenomPv2 NRTidalv2 (blue), and IMRPhenomXP NRTidalv2 (orange).The parameters shown here are the individual star masses MA,B, binary tidal deformability Λ, luminosity distance D, and inclination angle ι0.The 68% and 90% confidence intervals are indicated by contours for the 2D posterior plots, while vertical lines in the 1D plots indicate 90% confidence interval.We note a narrow constraint on the tidal deformability for IMRPhenomXP NRTidalv3 compared to the other models, due to the updated tidal information that was used.
For all the configurations discussed above, we calculate the mean mismatch µ F and maximum mismatch max( F ) for all waveform comparisons done.The results are shown in Table VI.

VI. PARAMETER ESTIMATION
Finally, to test the applicability of the developed NRTidalv3 models, we will reanalyze the two BNS detections GW170817 [2] and GW190425 [5].For this purpose, we utilize parallel bilby [97,98], which performs GW parameter estimation using a parallelized nested sampler named dynesty [99].Parallel bilby uses Bayes' theorem: where P (θ|d, H) is the posterior probability distribution of the parameters θ given some data d (which consists the waveform) and hypothesis H, p(θ|H) is the prior probability distribution, and E(d|H) is the evidence, serving as normalization constant to P (θ|d, H).Meanwhile, L(d|θ, H) is the likelihood of obtaining the data d given that the parameters θ are under the hypothesis H.For further details, we refer to Refs.[97,100,101].
In our study, we use GW170817 and GW190425 as our GW events, and we use the following waveform models: IMRPhenomXP NRTidalv3, IMRPhenomPv2 NRTidalv2, and IMRPhenomXP NRTidalv2, and also compare their results with each other.We also use low-spin and highspin priors for both GW events following the standard LIGO Scientific-, Virgo, KAGRA Collaboration (LVK) analyses [2,5,10,23,100,101]. .000000.00025 0.00050 0.00075 0.00100 0.00125 Posterior probability; GW170817, low-spin prior FIG.17. Inferred spin parameters for GW170817 from a lowspin prior (χ ≤ 0.05) using IMRPhenomXP NRTidalv3.Plotted here are the probability densities for the dimensionless spin components χ1 (left hemisphere) and χ2 (right hemisphere) relative to the orbital angular momentum L and tilt angles (a tilt angle 0 • means that the spin is aligned with the L).The plot was done using a reference frequency of 20 Hz.

A. GW170817
We show the results of the inferred marginalized 1D and 2D posterior probability distribution of a selection of source parameters for GW170817 in Fig. 16, for a lowspin prior with χ ≤ 0.05.In this figure, we show the individual masses of the stars M A,B , the luminosity distance D in Mpc, and the inclination angle ι 0 .We also include here dimensionless tidal deformability defined in terms of the individual masses and tidal deformabilities of the stars [10]: 5  .
(53) From Fig. 16, we note that the performance of IMRPhenomXP NRTidalv3 in terms of the inferred parameters is consistent with the results of IMRPhenomPv2 NRTidalv2 [23].
The main difference is the slightly tighter constraint of IMRPhenomXP NRTidalv3 for Λ compared to the other models (which are based on NRTidalv2).This is due to the reduction of the secondary peak that NRTidalv2 models show at higher Λ.
Finally, we investigate the spin constraints of IMRPhenomXP NRTidalv3 on GW170817 by plotting the inferred spin component magnitudes and orientation (in terms of tilt angles).A tilt angle of 0 • means that the spin components are aligned with the orbital angular momentum L. The results are shown in Fig. 17, where the left hemisphere is for χ 1 and the right hemisphere corresponds to χ 2 .We note that large negative components and anti-aligned spins are ruled out, which is consistent with previous studies using IMRPhenomPv2 NRTidalv2 [23].
For the corresponding parameter estimation for the high-spin prior (χ ≤ 0.5), we refer to Appendix C.

B. GW190425
The results for the marginalized posterior distributions for GW190425 can be found in Fig. 18, where (as in the case for GW170817), we see a slightly tighter constraint on the tidal deformability Λ from IMRPhenomXP NRTidalv3, than the other waveform models.However, all results are consistent between the different GW models.
Finally, the spin magnitudes and orientation inferred for GW190425 using IMRPhenomXP NRTidalv3 in Fig. 19 shows that negative values of the spin component magnitudes are heavily disfavored, as well as orientation greater than 90 • .The corresponding parameter estimation for the high-spin prior (χ ≤ 0.5) is found in Appendix C.

C. Performance of NRTidalv3
To further investigate the performance of IMRPhenomXP NRTidalv3 we plot the posterior probability distribution of the natural logarithm of its likelihood (also known as the log-likelihood), ln L, and compare it with the other two waveform models, as shown in the top panel of Fig. 20.We note that for GW170817, the ln L of IMRPhenomXP NRTidalv3 is generally shifted towards slightly larger values when compared to that of IMRPhenomPv2 NRTidalv2 and IMRPhenomXP NRTidalv2.In addition, we compute the Bayes factors (which is the ratio of the evidences or posterior probabilities) with respect to the null hypothesis of a non-detection and find IMRPhenomXP NRTidalv3 has ln B = 526.726± 0.085 which is larger than that of IMRPhenomPv2 NRTidalv2 (ln B = 526.647± 0.084) and IMRPhenomXP NRTidalv2 (ln B = 526.668± 0.079).Although these numbers might indicate a slight preference for IMRPhenomXP NRTidalv3, they are not statistically significant.
Similar to GW170817, we also present the loglikelihood for GW190425 in the bottom panel of Fig. 20.
As before, we find slightly larger ln L for IMRPhenomXP NRTidalv3, and at the same time Bayes Factors of IMRPhenomXP NRTidalv3 ln B = 53.849± 0.079, which is larger than that of IMRPhenomPv2 NRTidalv2 (ln B = 53.756± 0.079) and IMRPhenomXP NRTidalv2 (ln B = 53.766± 0.079), though with overlapping uncertainties and not being statistically significant.The marginalized 1D and 2D posterior probability distributions for selected parameters of GW190425, obtained with IMRPhenomXP NRTidalv3 (black), IMRPhenomPv2 NRTidalv2 (blue), and IMRPhenomXP NRTidalv2 (orange).The parameters shown here are the individual star masses MA,B, binary tidal deformability Λ, luminosity distance D, and inclination angle ι0.The 68% and 90% confidence intervals are indicated by contours for the 2D posterior plots, while vertical lines in the 1D plots indicate 90% confidence interval.As with GW170817, note a narrow constraint on the tidal deformability for IMRPhenomXP NRTidalv3 compared to the other models, due to the updated tidal information that was used.

VII. CONCLUSIONS AND OUTLOOK
Summary: In this work, we introduced NRTidalv3 as a new description of the tidal phase contribution to the total GW phase of BNS systems.This model improves upon the previous version (NRTidalv2) by employing a larger set of NR waveforms across a wide range of tidal deformabilities with mass ratios up to q = 2.0.To construct the model, we employed a representation of the frequency-domain, dynamical enhancement factor for the Love number, as well as a Padé approximant, which imposes additional constraints pertaining to the masses and tidal deformabilities of the component stars.
The model was then augmented onto existing BBH baseline models in LALSuite (IMRPhenomD, IMRPhenomXAS, IMRPhenomXP, and SEOBNRv5 ROM).To test the performance of NRTidalv3 we calculated its dephasing relative to existing NR waveforms and found overall consistency between NRTidalv3 with respect to the uncertainties in the NR waveforms and with respect to other waveform models.The computed mismatches between NRTidalv3 and NR waveforms, as well as NR hybrids, were smaller than for the previously constructed model NRTidalv2 [60].With respect to other tidal waveform models, we observe the largest mismatches for high masses, mass ratios, tidal deformabilities, and spin magnitudes.However, we find that overall NRTidalv3 can be employed for masses M A,B ∈ [0.5, 3.0]M ⊙ , dimensionless spins with magnitude below |χ A,B | ≤ 0.5, and for dimensionless tidal deformalities Λ A,B ∈ [0., 20000].
To finalize our investigations, we performed parameter estimation analyses for the GW events GW170817 and GW190425 using IMRPhenomXP NRTidalv3, IMRPhenomPv2 NRTidalv2, and IMRPhenomXP NRTidalv2.We find a slightly tighter constraint on the tidal deformability Λ for NRTidalv3 than NRTidalv2.For both events, the NRTidalv3 model displays a slightly higher log Bayes factor, but the difference is too small to be of statistical significance.In general, the performance of NRTidalv3 is consistent with previous analyses done by LVK on these GW events [23,101].Outlook: A suitable extension to NRTidalv3 would be the inclusion of higher-order modes since NRTidalv3 only includes the (2,2)-mode in its tidal phase description.It has recently been shown that the higherorder modes for BNS can be rescaled with respect to the (2,2)-mode, in the same manner as found on the BBH waveform phase contributions of higher-modes [32].Moreover, the inclusion of higher-order spin contributions to the BBH baseline phase can also be investigated.Aligned with the inclusion of higher modes in NRTidalv3 models would be the extension of existing BHNS models such as IMRPhenomNSBH [102] or SEOBNRv4 ROM NRTidalv2 NSBH [103] with NRTidalv3 phase contributions.
Another possible extension of the model would be a generalization of the frequency dependence of the Love numbers to be able to set the fundamental f-mode frequency f 02 as a free parameter.In the current implementation, this cannot be straightforwardly done, as the time-domain effective enhancement factor k eff 2 (and subsequently the frequency-domain keff 2 ) depends on the quasi-universal relation Eq. ( 14) between f 02 and Λ 2 .
Finally, other NR waveforms may also be added in a future version of the model, such as BNS and/or NSBH waveforms with [104] and without [105] sub-solar-mass components, as well as waveforms whose NS have more exotic EoS.This will allow the model to accommodate an even wider range of neutron star properties and to constrain these EoSs with GW observations.We present in Table VII the configurations of all the 55 NR waveforms from SACRA [28,29,72] and CoRe (using the BAM code) [31][32][33] that were used in the calibration of NRTidalv3.We include in the table the individual masses of the stars M A,B (M ⊙ ), the total mass M (M ⊙ ), mass ratio q, dimensionless tidal deformabilities Λ A,B , and radii R A,B in km.The radii data for the SACRA waveforms were adapted from Ref. [28,29,72], while the radii for the CORE (BAM) waveforms were computed using the adiabatic Love number and tidal deformability, i.e. from Λ = (2/3)k 2 R 5 /M 5 .We show here a comparison between the frequencydomain effective enhancement factor keff 2 numerically calculated from Eq. ( 25) and the effective representation keff 2,rep in Eq. ( 27).The results are shown in Fig. 21.We also show the relative error between keff 2 and keff 2,rep .The fitting formula satisfies the numerically calculated keff 2 , with a maximum (deviation) error of 1.8% and mean error of 0.2%.

Appendix C: Parameter Estimation for high-spin priors
We present here the results for the parameter estimation runs for GW170817 and GW190425, as was done in Sec.VI, but this time, using high-spin priors.For both events, we use a spin prior of up to χ ≤ 0.50.We do not employ spins up to 0.89 as done in some LVK studies, given that this spin is well above the breakup spin of neutron stars and is unrealistically large for BNS systems.
a. GW170817.The inferred properties of GW170817 with a high-spin prior, together with the marginalized 1D and 2D posterior distributions are shown in Fig. 22.As with the low-spin prior case, the result is consistent with the other waveform models and with previous results [23], and we observe a slightly narrow constraint for IMRPhenomXP NRTidalv3 in the tidal deformability Λ.As for the inferred spin parameters, Fig. 23 shows that relatively large components of the spin which are aligned or anti-aligned with the orbital angular momentum L are heavily disfavored, and constraints are placed on the in-plane spin components.
b. GW190425.We show inferred properties of GW190425 with a high-spin prior, together with the marginalized 1D and 2D posterior distributions in Fig. 24.As with the low-spin prior case, the result is consistent with the other waveform models, and we observe a narrower constraint for IMRPhenomXP NRTidalv3 in the tidal deformability Λ. Fig. 25 shows similar behavior as in Fig. 23 such that relatively large components of the spin which are aligned or anti-aligned with the orbital angular momentum L are heavily disfavored, and constraints are placed on the in-plane spin components [5].
c. Performance.Finally, the posterior probability distribution of the log likelihood ln L of the three models used with GW170817 and GW190425 and shown in Fig. 26.We note that for both GW events, the ln L values for IMRPhenomXP NRTidalv3 are shifted towards larger values relative to the other two models.For GW170817, we find ln B = 525.[28,29,72], while the last nine are from the CoRe database simulated with the BAM code [31][32][33].For each waveform, we indicate its equation of state (EoS), the masses MA,B(M⊙) of the individual bodies, the total mass M (M⊙), the mass ratio q = MA/MB (MA > MB), the dimensionless tidal deformabilities ΛA,B, effective tidal deformability Λ, and the radii RA,B [km].

NR Waveform Name
EoS    22.The marginalized 1D and 2D posterior probability distributions for selected parameters of GW170817, with a high-spin prior (χ ≤ 0.5) obtained with IMRPhenomXP NRTidalv3 (black), IMRPhenomPv2 NRTidalv2 (blue), and IMRPhenomXP NRTidalv2 (orange).The parameters shown here are the individual star masses MA,B, binary tidal deformability Λ, luminosity distance D, and inclination angle ι0.The 68% and 90% confidence intervals are indicated by contours for the 2D posterior plots, while vertical lines in the 1D plots indicate 90% confidence interval.We note a narrow constraint on the tidal deformability for IMRPhenomXP NRTidalv3 compared to the other models, due to the updated tidal information that was used.  1 .5 3 .0 ι 0 GW190425, high-spin prior FIG.24.The marginalized 1D and 2D posterior probability distributions for selected parameters of GW190425, obtained with IMRPhenomXP NRTidalv3 (black), IMRPhenomPv2 NRTidalv2 (blue), and IMRPhenomXP NRTidalv2 (orange), for a high-spin prior (χ ≤ 0.50).The parameters shown here are the individual star masses MA,B, binary tidal deformability Λ, luminosity distance D, and inclination angle ι0.The 68% and 90% confidence intervals are indicated by contours for the 2D posterior plots, while vertical lines in the 1D plots indicate 90% confidence interval.We note a narrow constraint on the tidal deformability for IMRPhenomXP NRTidalv3 compared to the other models, due to the updated tidal information that was used.FIG. 26.The distribution of the log-likelihood ln L for GW170817 (top) and GW190425 (bottom) with high-spin priors (χ ≤ 0.5).The vertical lines indicate 90% confidence intervals for the models.We note the shifting of the distribution using IMRPhenomXP NRTidalv3 towards larger ln L.

FIG. 1 .
FIG. 1. Distribution of the masses (where we set MA > MB) and corresponding tidal deformabilities ΛA,B of the neutron stars used in the NR simulations (see TableVII).The plot shows overlapping data points for systems with the same masses, but can have different tidal deformabilities due to the various EoSs that were employed.

FIG. 2 .
FIG. 2.An example of a hybrid Waveform (SACRA:15H 135 135 00155 182 135).Here, we plot the real part of the (2,2) mode rh22 of the GW strain (where r is the extraction radius) as a function of the retarded time u.The blue curve is the NR waveform, the orange dashed curve is the EOB-BNS waveform from SEOBNRv4T, and the black, dotted curve is the hybridized waveform EOB-NR.The gray, dashed, vertical lines denote the boundaries of the hybridization window, [t1, t2] corresponding to ω ∈ [0.035, 0.04].The peak amplitude of the NR and EOB-NR hybrid, indicating merger, is set at u/M = 0.

FIG. 4 .
FIG. 4. The universal relation between the f -mode frequency (in kHz) and Λ used for stars A and B.

FIG. 7 .
FIG. 7. Top:The fits (dashed lines) plotted on top of the time-domain hybridized tidal phase data (solid lines).Bottom: Fractional differences between the fits and the data.The EOB-NR hybrid tidal phases are color-indicated by their EoS.

FIG. 8 .
FIG. 8. Top:The fits plotted on top of the frequency-domain hybridized tidal phase data.Bottom: Fractional differences between the fits and the data, where ∆ψT = ψ Fit T − ψ Data T .

T
are found again early in the frequencies and around the hybridization window.Furthering the investigation of the performance of NRTidalv3, we plot the difference between the absolute fractional differences of NRTidalv2 and NRTidalv3 with respect to the data, that is, |∆ψ NRT2 T | − |∆ψ NRT3 T |, against ω, and show the results in Fig. 9.We also indicate in the figure |∆ψ NRT2 T | − |∆ψ NRT3
FIG. 12. Mismatch histograms for different BBH baselines augmented with the NRTidalv3 model, versus the BBH baseline model added with the EOB-NR hybrid tidal phase.Vertical lines denote the mean of the histograms with the same color.The mismatches with NRTidalv2 are also shown for comparison.For each BBH baseline + NRTidalv3 model, mismatches are computed starting from fmin = 20 Hz in the left panel, indicating the mismatch for the full waveform where EOB contribution is present in the early inspiral.Mismatches are also computed for fmin = 350 Hz in the right panel, where typically the hybridization starts and where the NR contribution becomes dominant up to merger.

FIG. 13 .
FIG. 13.Mismatch Comparisons of NRTidalv3 with other tidal waveform models for non-spinning configurations.Each comparison with either TEOBResumS or SEOBNRv4T contains 4000 random configurations of mass MA,B = [0.5, 3.0] and tidal deformability ΛA,B = [0, 5000].The rest of the comparisons are done with ΛA,B = 20000.For each subfigure we include a color bar indicating the values of the mismatches.Note that the maximum value in the color bar for IMRPhenomXAS NRTidalv3 vs SEOBNRv5 ROM NRTidalv3 is different from the rest of the subfigures due to the NRTidalv3 models yielding very small mismatches with respect to each other.

10 F
FIG. 14. Mismatch Comparisons of NRTidalv3 with other tidal waveform models for aligned-spin configurations.Each subfigure has 4000 random configurations of mass MA,B = [1.0,3.0], tidal deformability ΛA,B = [0, 5000] and spin χA,B = [−0.5,0.5].We also put the mismatches with NRTidalv2 for comparison.For each subfigure we include a color bar indicating the values of the mismatches.
FIG.18.The marginalized 1D and 2D posterior probability distributions for selected parameters of GW190425, obtained with IMRPhenomXP NRTidalv3 (black), IMRPhenomPv2 NRTidalv2 (blue), and IMRPhenomXP NRTidalv2 (orange).The parameters shown here are the individual star masses MA,B, binary tidal deformability Λ, luminosity distance D, and inclination angle ι0.The 68% and 90% confidence intervals are indicated by contours for the 2D posterior plots, while vertical lines in the 1D plots indicate 90% confidence interval.As with GW170817, note a narrow constraint on the tidal deformability for IMRPhenomXP NRTidalv3 compared to the other models, due to the updated tidal information that was used.
FIG.19.Inferred spin parameters for GW190425 from a lowspin prior (χ ≤ 0.05) using IMRPhenomXP NRTidalv3.Plotted here are the probability densities for the dimensionless spin components χ1 (left hemisphere) and χ2 (right hemisphere) relative to the orbital angular momentum L and tilt angles (i.e. a tilt angle 0 • means that the spin is aligned with the L).The plot was done at a reference frequency of 20 Hz.

FIG. 20 .
FIG.20.The distribution of the log-likelihood ln L for GW170817 (top) and GW190425 (bottom).The vertical lines indicate 90% confidence intervals for the models.We note the shifting of the distribution using IMRPhenomXP NRTidalv3 towards larger ln L.
FIG.23.Inferred spin parameters for GW170817 from a highspin prior (χ ≤ 0.5) using IMRPhenomXP NRTidalv3.Plotted here are the probability densities for the dimensionless spin components χ1 (left hemisphere) and χ2 (right hemisphere) relative to the orbital angular momentum L and tilt angles (i.e. a tilt angle 0 • means that the spin is aligned with the L).The plot was done at a reference frequency of 20 Hz.
FIG.25.Inferred spin parameters for GW190425 from a highspin prior (χ ≤ 0.50) using IMRPhenomXP NRTidalv3.Plotted here are the probability densities for the dimensionless spin components χ1 (left hemisphere) and χ2 (right hemisphere) relative to the orbital angular momentum L and tilt angles (i.e. a tilt angle 0 • means that the spin is aligned with the L).The plot was done at a reference frequency of 20 Hz.

TABLE I .
(20)meters for the NRTidalv3 timed-domain approximant.The table lists the fitting parameters for the Padé approximant, given in Eqs.(18)and(20).

TABLE III .
(32)meters for the NRTidalv3 frequency-domain approximant.The table include the fitting parameters for the dynamical love number enhancement, as well as the respective fitting parameters for the Padé approximant given in Eq.(32).

TABLE IV .
Available NRTidal approximants in LAL.

TABLE V .
BNS configurations for the time-domain dephasing comparisons.We show the configurations of NR simulations using BAM, SACRA, and SpEC that are used in the comparison of waveforms in the time domain.The SACRA waveforms, as well as BAM:0001, BAM:0037 and BAM:0064, are also used in the calibration in NRTidalv3.For each NR simulation listed here, we indicated the EOS used, the individual masses MA,B, spins χA,B, tidal deformabilities ΛA,B, effective tidal deformability Λ of the system, resolution h fine , eccentricity e, and whether or not a Richardson extrapolated waveform was constructed from the simulations.

TABLE VI .
The mean mismatch µ F and maximum mismatch max( F ) for various waveform model comparisons for the nonspinning, aligned-spin, and precessing configurations.A '−' symbol indicates that the mismatches were not calculated between the two waveform models.

TABLE VII .
Properties of the 55 NR waveforms used in the calibration of NRTidalv3.The first 46 waveforms are from SACRA