Rising Tides: Analytic Modeling of Tidal Eﬀects in Binary Neutron Star Mergers

The gravitational waves produced by binary neutron star mergers oﬀer a unique window into matter behavior under extreme conditions. In this context, we model analytically the eﬀect of matter on the gravitational waves from binary neutron star mergers. We start with a binary black hole system, leveraging the post-Newtonian formalism for the inspiral and the Backwards-one-Body model for the merger. We combine the two methods to generate a baseline waveform and we validate our results against numerical relativity simulations. Next, we integrate tidal eﬀects in phase and amplitude to account for matter and spacetime interaction, by using the NRTidal model, and test its accuracy against numerical relativity predictions, for two equations of state, ﬁnding a mismatch around the merger. Subsequently, we lift the restriction on the coeﬃcients to be independent of the tidal deformability, and recalibrate them using the numerical relativity predictions. We obtain better ﬁts for phase and amplitude around the merger, and are able to extend the phase modeling beyond the merger. We implement our method in a new open-source Python code, steered by a Jupyter Notebook. Our research oﬀers new perspectives on analytically modeling the eﬀect of tides on the gravitational waves from binary neutron star mergers.


Introduction
Neutron stars are extremely dense remnants of massive stars, at the brink of collapsing into a black hole, with masses comparable to that of the sun, contained in only twenty kilometers diameter.They are captivating celestial objects, under the action of powerful gravitational and magnetic fields, which cannot be replicated on Earth.These exceptional properties present them as excellent astrophysical laboratories for investigating the behavior of matter in extreme conditions of density, pressure and temperature.When two neutron stars collide and merge, they emit gravitational waves, accompanied by electromagnetic radiation, matter, and neutrinos; carrying valuable information about their masses, sizes and interior structure.
The first direct detection of a gravitational wave (GW) signal from a binary neutron star (BNS) collision, named GW170817 [1], was accompanied by a gamma-ray burst [2] and ignited a new explosion, called a kilonova, powered by the radioactive decay of the heavy elements synthesized in the merger [3].This twin detection marked the onset of the golden era in neutron star research, alluding to the precious metals such as gold synthesized and expelled during the collision, and has inspired intense theoretical investigations.These studies provided us with a wealth of new insights into the yet unknown internal structure of the neutron stars [4][5][6], the characteristics of their magnetic fields [7], and the outflow of heavy matter during collision [8].
Unfortunately, simultaneous detections of gravitational waves and light are few and far between.Out of the almost 100 such events reported to date [9], the overwhelming majority came from binary black hole (BBH) collisions and only one other, named GW190425, was produced by an unusually heavy BNS collision [10,11].The two BNS mergers detected so far raised many questions, proving again that nature is more complicated than our models [12,13].It is thus imperative to revise our current understanding of the physics involved in those models, in order to gain insight on BNS mergers.This is a timely topic, because in the next decade hundreds of such GW events might be detected [14], with the advent of the new generation of gravitational wave observatories, including both ground-based detectors such as Cosmic Explorer [15] and the Einstein Telescope [16], and space-borne instruments such as LISA citearXiv:2001.09793,DECIGO [17] and TianQin [18].
In this study, we employ the NRTidal approximant [30,32,33], an elegant model that adjusts the analytically calculated binary black hole (BBH) waveforms by adding a closed-form expression to account for the tidal influences in the phase and amplitude of the GW.This model was used for estimating source properties [36] and constraining the equation of state for ultra-dense matter in the first two BNS detections [37], and it is implemented in the LIGO Scientific and Virgo Collaborations Algorithm Library Suite (LALSuite) [38].
Our goal is to analytically model the effect of tides on the GW signal during the inspiral and to accurately extend it through the merger of the BNS systems considered.We aim to create an open-source, easily replicable code that generates comprehensive analytical templates for gravitational radiation during BNS collisions.By doing so, we facilitate independent consistency checks and assist in defining the domains of validity for the approximation models employed.Ultimately, our work contributes to the development of a universal analytical model encapsulating the BNS dynamics and the effect of tides on GW emission.
This paper unfolds as follows: first, we present the point-particle system that characterizes BBH collisions, and calculate the GW by combining the post-Newtonian formalism for inspiral evolution with the Backwards-one-Body model for the merger.As baseline GW we use the fully analytical model for BBH collision we developed in [39] and expanded in [40], based on [41] for the inspiral and [42] for the merger.We prove our model's validity and efficiency through a comparison with numerical relativity (NR), utilizing the SXS:BBH:0180 template from the Simulating eXtreme Spacetimes (SXS) catalogue [43].We then incorporate the tidal deformability into the point-particle waveform, both as polynomial corrections and as rational functions approximations to the phase and amplitude, by using the coefficients proposed in [32,33].Next, we test our implementation against SXS:NSNS:0001/0002 for two equations of state [44].
Finally, we push this model past the merger, by performing a new fit to the numerical BNS data for the tidal phase and amplitude, from which we derive updated values for the polynomial coefficients, expanding thus its applicability.Upon determining the new coefficients, we reconstruct the tidal correction in phase and amplitude beyond the merger and analyze the tidal influence on the early stages of the collision, revealing the effect of the matter interaction on the system's orbits.
In this work we consider two neutron stars with total mass M = M A + M B , M A ≤ M B and mass ratio q = M B /M A ≥ 1.We express time, space and energy in geometric units (with G = c = 1), in terms of the binary mass M , written as a multiple of the sun mass M ⊙ .

The Baseline Model
We start with assembling the baseline analytical model for the calculation of the GWs from a BBH collision, by following the common procedure as we detailed in [39,40].First, we split the binary motion in two regions: the weak field, during the inspiral, and the strong field, during the merger, and apply different mathematical formalisms to obtain the waveform for each region.We then generate the complete GW template for the whole binary evolution by matching those two regions in frequency around the last stable orbit, and building the hybrid waveform.Lastly, we compare our model against a numerically generated equal mass BBH collision, to uphold its validity.

The Baseline Inspiral Model
Let us consider a tight binary system of separation r, in quasi-circular orbit.By defining the reduced mass µ = (M A M B )/M and the symmetric mass ratio η = µ/M , we reduce it further to a single particle of mass µ and position r, orbiting around the mass M located at the center of mass.We require the system to dissipate GWs and are led to the balance equation: This formula states that the GW flux F (t) is emitted at the expense of the orbital energy E(t), causing the orbit to shrink.Eq.( 1) can be rewritten in terms of a small factor x pN = (v/c) 2 called post-Newtonian parameter, where v is the orbital velocity of the binary and c is the speed of light (see [45][46][47][48]).
Using the post-Newtonian (pN) approximation, we expand the deviation from Newtonian gravity as a perturbation in power series of x pN .From the many different methods of solving eq.( 2), in this work we chose the TaylorT4 approximant, shown to agree best with numerical simulations [49].Within this method, eq.( 2) becomes [39,41]: where N/2 denotes the pN expansion order.In this work we go up to 3.5 in the leading pN order, adding self-force and hereditary correction terms up to 6pN order, to increase the accuracy in modeling the region near the merger [41].By integrating eq.( 3) with the coefficients ξ j we obtain the evolution of the pN parameter x pN .Next, we use Kepler's third law v 2 = (M Ω(t)) 2/3 where Ω is the angular orbital velocity, to obtain the equation for the orbital phase: We readily integrate eq.( 4) to find the time evolution of the orbital inspiral phase.We reuse Kepler's third law written as: v = M/r to extract the orbital separation as r = M/x pN , which we then expand in powers series of x pN for increased accuracy: Lastly, we calculate the GW amplitude for optimal orientation of the source, Now we have all the pieces necessary to construct the dimensionless strain, defined as: Mathematically, the strain is decomposed into two transverse quadrupolar polarization modes, with h + representing the real part of eq.( 7) and h × the imaginary part, This is because the GWs compress the space in one direction while simultaneously stretching it in the orthogonal direction, such that the signal goes twice through maxima and minima during one orbital cycle, making the frequency of the GW twice the orbital frequency.This technique, albeit powerful, is valid only if the gravitational field is sufficiently weak and the orbital velocity is smaller than the speed of light.

The Baseline Merger Model
Going beyond the pN approximation brings us up against the strong gravitational field around the merger.The transition between the weak and strong field is marked by the innermost stable circular orbit (ISCO), defined as the last stable orbit a particle would have when orbiting around a black hole.When the masses of the celestial objects are comparable, this location is not well defined [50], but can be approximated with the last stable photon orbit, called the light-ring (LR), which is close to the peak of the curvature potential [51].From there on, through the merger and ringdown, we use the Backwards-One-Body formalism [42].This technique starts from the perturbed final black hole resulting after collision and builds the GW signal back in time to the end of the inspiral, assuming we know the mass, spin and ringdown frequency of the remnant.In this model, the strain of the GW signal can be modeled analytically as exponentially decaying sinusoids that break free from the LR on null geodesics [52]: A lmn e iω lmn t e −t/τ lmn , where l is the principal, |m| ≤ l the azimuthal and n is the overtone index of a mode.This kind of perturbation is called quasinormal (QNM) ringing with frequency ω lmn and damping time τ lmn .The frequency of the lower mode is (l = 2, m = 1, n = 0) coincide with the orbital frequency and we will denote it as Ω QN M .The (l = 2, m = 2, n = 0) mode can be derived from this mode with the simple relation ω 22 = 2Ω QN M .This mode carries away most of the GW's energy (≈ 95%) [53], while the higher harmonics, being much quieter, can be usually neglected.We will drop the (l, m, n) indices and consider that the strain is well described by the dominant mode.The BoB model requires the QNM frequency and damping time, which are determined by the spin and mass of the final black hole.We will estimate the final spin of the resulting black hole with a polynomial of coefficients s ij , as given in [54,55]: We define an effective spin: with χ A,B = S A,B /M 2 A,B the dimensionless individual spins and S A,B the spin angular momentum of each black hole entering the merger.In our calculations we pick χ A,B = 0, as in the NR simulations we use for comparison, although the mean spin for lowspinning, astronomical neutron stars is χ A,B ≈ 2 × 10 −3 [56].
For the final mass we use the fit to NR given in [57] for comparable mass binaries, where ẼGW = E 0 + E 2 χ 2 f + E 4 χ 4 f is the dimensionless energy released in GWs.For the disk mass we take as upper limit of M disk ≈ 10 −2 M [58].The coefficients for E i are given in Table III of [57].
Next, we calculate the dominant resonant frequency with a polynomial fit: and use a similar formula for the quality factor: The pairs (f i , q i ) are taken from Table VIII of [53].Lastly, the damping time is: The orbital angular frequency Ω BoB (t) is given in this model by: Here, Ω i is the initial frequency, t 0 is the time at which the strain of the GW reaches its peak amplitude, and t i the initial time , marking the transition between the weak and strong regime.The parameter κ in eq.( 16) assures the continuity between the end of the inspiral and the beginning of the merger, and is given by: The essential variable in the BoB model is the initial time t i , that locks in the frequency Ω i at the beginning of the merger to the frequency at end of the inspiral: We obtain the phase Φ BoB (t) by integrating eq.( 16) between t i and a final time t f : The amplitude of the GW signal is modeled with the simple function: where A 0 is a scaling factor.This amplitude is taken to correspond to the Weyl scalar |Ψ 4 (t)|, which is related to the strain by the formula The model rests upon the assumption that the amplitude changes much slower compared to the phase during the merger, and uses the simple expression for the strain: where ω BoB (t) = 2Ω BoB (t) is the frequency of the dominant mode.Please note that in this model the influence of the proportionality coefficients arising from the double integration of eq.( 21) over time is overlooked.In reality, both the magnitude and the peak location of the amplitude are affected by this integration.As a result, we expect this model to yield a slight time deviation in predicting the merger position.
We assessed the strengths and weaknesses of the BoB model in a previous work [59].

The Complete Baseline Model
Let us now assemble the two forms of the GW strain by gluing the waveform around the LR, at the transition between the weak and strong field, estimated by [60]: The frequency and phase matching is robust against slight variations around this location, while the continuity of the amplitude is best ensured by choosing r LR ≈ 4M , representing the location of the retrograde photon orbiting a maximally rotating BH.
We stitch the amplitudes of the two models together by first normalizing the BoB amplitude with its peak value and then rescaling the pN amplitude with the normalized BoB amplitude at t i , by using the formula: where the quantities with an overhead bar represent their rescaled values.We build the complete hybrid model by gluing together the waveforms of the two domains at t i with the following simple piecewise (step) function: For example, the hybrid orbital frequency of the GW has the form: We form the hybrid amplitude, phase and strain of the GW with the same technique.Before moving forward, we must test the accuracy of our implementation with a numerical waveform.We chose SXS:BBH:0180, corresponding to an equal-mass, nonspinning BBH configuration of normalized mass and initial separation of 18M in geometric units [43].We use the same configuration, but start our analytic model at a separation of 15M , which for a regular equal-mass BNS with a total mass of 2.8M translates to a distance of abut 45km between the centers of mass.We start the BoB model with the initial orbital frequency (Ω i , Ωi ) of the pN model at the end of the inspiral marked by the LR, and build the hybrid baseline model.Next, we interpolate NR data onto a time array of constant time steps, as used by the analytical model, and align the strain at the beginning.This alignment starts at the initial frequency of the analytical model and extends over a duration equal to the initial complete period of the strain.
We plot in Fig. 1 a comparison between the evolution of the h + strain component obtained with our hybrid model, and the + polarization of the GW strain from the NR simulation, as well as the differences in phase and amplitude (lower insets).We observe a significant overlap between the strain from the analytical model and the NR simulation throughout the entire inspiral, extending up to the LR.The faster merger occurring in the analytical model is primarily attributed to the Post-Newtonian (PN) approximation's limited reach into the strong field domain.In subsequent studies we aim to refine this model by including higher-order post-Newtonian (pN) terms as they become available.The shift in the position of the merger is accentuated by the BoB model because, as we explained earlier, it assumes a simplified expression for the amplitude.We analyzed this deviation from the NR predicted waveform at the merger in [59], and in future work we intend to correct the merger amplitude by incorporating the neglected integration coefficients.The phase concordance between the analytical model and the NR data is consistent throughout the entire waveform during the inspiral phase, remaining within the sub-radian accuracy of the BBH numerical waveforms [44].We see an increase in the phase difference above 1 rad only at the LR and reaching 2.8 rad by the time of merger, which is within the numerical BNS dephasing due to tidal disruptions [44].Similarly, the amplitude correlates well up to the LR, with the discrepancy becoming more pronounced as the merger approaches.The addition of tidal effects to this analytic model will overwrite the errors around the merger, because a BNS system reaches merger earlier than its corresponding BBH.   1 Comparison between the numerical and the analytic plus component of the strain, matched over a period around the initial frequency.The upper inset shows that the strain peaks are slightly displaced, but the merger is well modeled if we overlap them (dotted red curve).The dashed vertical yellow line marks the LR.The two lower insets contain the difference in phase (lower left) and amplitude (lower right) between our model and SXS:BBH:0180.Only expected small differences before the merger are visible.

The BNS Model
Thus far, we have not taken into consideration the internal structure of the neutron stars, which is determined by the equation of state and accounts for the effect of the tidal interactions.It is expected for all neutron stars in the universe to be described by a unique equation of state, but its expression is not currently known, only estimated by theory [61] and experiments [62] in nuclear physics.The tidal effects become important in the late stage of the inspiral, inducing deformations in the stars and driving them to merge faster, thus affecting the emitted GW signal both in phase and in amplitude.It is this effect that will allow us to determine the equation of state corresponding to the dense matter inside neutron stars from the direct observations of the amplitude and phase of the gravitational waves emitted during BNS collisions.
We can encode the tidal effects into the tidal correction strain: where A T and φ T are the analytic tidal corrections to the GW amplitude and phase.We find the analytical GW strain of the BNS model by multiplying eq.( 27) to the strain of the baseline model [32]: This equation tells us that we must add the tidal phase correction to the baseline, while for the amplitude, we should multiply the tidal correction to the baseline: We implement the tidal effects as polynomial functions of the variable x = (M Ω hyb ) 2/3 and build the analytical expressions for the BNS phase and amplitude provided by the pN and NRTidal models, then we validate our implementation by comparing it two open-source BNS simulations from the SXS catalogue [44] .

The BNS Tidal Phase
As the stars get closer together, their tidal deformability increases, affecting the binding orbital energy.Looking back to eq.( 1), we should not be surprised to see that this will change the flux of gravitational waves emitted as the orbit shrinks.The tidal correction becomes significant only as the stars approach the merger, not entering in the pN decomposition until the 5 th pN correction term [63].The correction to the GW phase induced by the tidal deformation of one star due to the gravitational pull produced by its companion is known up to 7.5-pN [64], and can be written as: with x from the baseline BBH model and the pN tidal phase coefficients c [A,B],i depending on the mass ratio, as given in [33].The tidal effects are encoded in κ [A,B],eff , named the effective tidal coupling constant, introduced in [30]: where: ΛA = 16 13 ; ΛB = 16 13 are the symmetric mass-weighted tidal deformabilities [27], defined in terms of Λ [A,B] , the dimensionless tidal deformability of each star in the binary [65].Let us take a closer look at this quantity, because it is its dependence on the equation of state and mass that can reveal the internal structure of the neutron stars in the binary by measuring the phase of their GW emission and fitting it with analytical or numerical templates.In general, the nuclear models give numerical values for the dimensionless tidal deformability between 10 2 and 10 3 , varying inversely proportional with the mass and the compactness of the star [66].So called hard nuclear equations of state models describe less compact stars, predicting higher values for the tidal deformability, while nuclear models with soft equations of state favor more compact stars and tidal deformabilities towards lower values.The dimensionless tidal deformability calculated from the GW170817 event gave, for a neutron star of mass 1.4M ⊙ , a value somewhere in the middle range, namely Λ ≈ 400 [67].
We now return to the tidal correction to the GW phase, and, as mentioned before, choose the NRTidal model for supplying the expression of the tidal phase evolution.This analytical approximation, first introduced in [30] and subsequently improved in [32,33], is calibrated to NR, and because of the scarce availability of numerical simulations of BNS systems with unequal mass ratio, includes only non-spinning, equal-mass configurations.In this case, eq.( 31) simplifies to a polynomial of constant coefficients: In order to fit eq.(34) to NR data, NRTidal replaces the polynomial P (x) with a rational function R(x) given by a Padé approximant of constant coefficients: determined from the NR, and restricted to enforce consistency with the analytical coefficients entering the polynomial P (x).
Using one set of coefficients for c i and two sets of coefficients for (n i , d i ) from [30,32,33], we implement the tidal correction to our hybrid baseline model as: Note that, in contrast to the corresponding equations presented in [30,32,33], we found necessary to use a positive sign in eq.(36), in order to obtain the expected behavior for the BNS phase evolution when we use eq.( 29) to obtain the analytic BNS phase.
It is time to check our implementation and test how well does the NRTidal model predict the tidal phase against numerical simulations.We choose two equal-mass, nonspinning BNS systems with two different equations of state [44], publicly available in the SXS GW catalogue [43].The first numerical GW used is SXS:NSNS:0001, with mass M 1 = 2.8 and an ideal gas equation of state of polytropic index Γ = 2 corresponding to a dimensionless tidal deformability ΛΓ2 = 791.The second is SXS:NSNS:0002, with mass M 2 = 2.7 and a piecewise polytropic equation of state for cold dense matter called MS1b [68], with dimensionless tidal deformability ΛMS1b = 1540.While these tidal deformabilities surpass the one estimated from the GW170817 observation, they are not yet ruled out and are useful as an upper bound.We also note that in the pN and NRTidal models the tidal correction to the phase depends only linearly on the dimensionless tidal deformability, and the coefficients should remain independent of it.
We start with the assumption that at a considerable distance from the merger, the tidal effects are negligible and the BNS signal is virtually indistinguishable from that of the BBH.This allows us to rescale the phase and shift the time of the numerical BNS simulations to start from zero at the reference frequency of each simulation, marking when the initial burst of spurious transient radiation has dissipated in the NR simulation [43].Then, we find the time and phase corresponding to the same frequency in the numeric BBH model and align the numerical BNS waveforms with the corresponding values of the numeric BBH model.We take the difference in phase between the numeric BNS and BBH data, up to the merger frequency, taken where the numeric BNS amplitude reaches its peak.
We plot in Fig. 2 the phase comparison between our analytic baseline hybrid model, SXS:BBH:0180, SXS:NSNS:0001 (BNS1) and SXS:NSNS:0002 (BNS2).The insets in Fig. 2 contain the analytical models for the BNS tidal phase φ T,pN and φ T,NRfit from eq.( 36), compared to the numerical BNS tidal phase.We mark with a dashed yellow curve the phase of 2 rad, representing the upper limit of uncertainty in phase due to numeric errors [44].The maximum modeling error of our baseline BBH phase in the domain considered is below 0.5 rad.We observe that the NRTidal phase (dash-dotdot red and magenta) slightly overestimates both the pN (dotted maroon) and the NR phase, and is more accurate for the BNS1 case, corresponding to the smaller tidal deformabiliy.This suggests that the assumption of the coefficients being independent of the tidal deformability may be too restrictive.In Section 4 we will delve into the behavior of the phase around the merger and propose new coefficients for the tidal correction to push the model past the merger.

The BNS Tidal Amplitude
The amplitude of the GW is extremely small, and its increase is less dramatic than the phase accumulation observed up to the merger, defined as the point of maximum amplitude [69].This makes the task of accurately modeling analytic tidal contributions to the amplitude more difficult compared to the phase.Taking advantage of the relatively small change in the GW amplitude, we can effectively use a low-order Mϕ ϕNS2 F2 Fig. 2 Phase comparison between the baseline model (solid orange), SXS:BBH:0180 (solid black), SXS:NSNS:0001 (dashed blue) and SXS:NSNS:0002 (dash-dot green) overlapped at the reference frequencies of the NR simulations.The insets contain our analytical BNS phase calculated with the pN (dotted maroon) and NRTidal (dash-dot-dot red and magenta) corrections, plotted up to the merger against the numerical phase difference for the two cases.The horizontal dotted yellow curve marks the maximum uncertainty in phase of 2 rad.
polynomial within the context of the pN approximation, of the form [33,64]: where D L is the distance from the detector to the BNS system and ĉ[A,B],i are the pN tidal amplitude coefficients, given in [33,64].Let's assume again a similar-mass binary, where only the symmetric mass ratio η retains the information on the mass ratio.Now, the pN tidal correction in amplitude eq.( 37) simplifies to: As expected, this approximation becomes less reliable as the binary approaches the merger.The NRTidal model corrects it by adding a dependence of x: where the parameter d is fixed by the identity [33]: We compute the merger amplitude A mrg using an analytically predicted quasiuniversal relation valid at the moment of merger [70], that provides the peak amplitude only as function of the tidal coupling constant κ eff [33]: As a cautionary remark, eq.( 40) requires that we provide an analytical expression for x mrg as well, without relying on numerical data.A good approximation is to assume that the stars merge when they come in contact, and x mrg is given by: where R A,B are the radii of the stars, for which we use the value of 11.5km as specified in [44].Please note that this simplified expression does not account for the deformations of the stars.To address this, we use an alternative route, and calculate x mrg = (M Ω mrg ) 2/3 from the analytic expression for the merger frequency given in [71]: We chose to use the result calculated with eq.( 43) in our implementation.Working within the same assumption that far from the merger the BNS and BBH waveforms have indistinguishable amplitudes, we rescale the numeric BNS amplitude to coincide with the numeric BBH amplitude for a binary with the same total mass at the reference frequency.We align the BNS and BBH amplitudes and plot them in Fig. 3.This is not fulfilled exactly, as we see in the inset of Fig. 3, where we plot a comparison of their h + strain because at the BNS reference frequency, tidal effects are already present in the NR waveforms.
Let us now add the tidal amplitude correction to our baseline model, using eq.( 30).We assume that the amplitudes are indistinguishable at the reference frequency of the numerical BNS data, although the waveform is already in the regime where tidal effects are non-negligible.Nonetheless, we consider this error in amplitude as not essential for our analysis and consider that A T,ref ≈ 1.We capture this behavior for the evolution of our analytic BNS amplitude, by changing A T to: We plot in Fig. 4 the ratio of the amplitudes for SXS:NSNS:0001 (BNS1) and SXS:NSNS:0002 (BNS2) to the numeric BBH amplitude compared with the ratio of our hybrid analytic baseline amplitude o the numeric BBH amplitude.The insets in Fig. 4 contain the analytical approximations for the BNS tidal amplitude, as in eq.( 38) and eq.( 39), up to the corresponding BNS merger frequency.We observe that both analytical approximations underestimate the steep increase in amplitude around the merger.We will return to this comparison in Section 4, when we will make a new fit for the analytical amplitude to the numerical data and taper the amplitude past the merger using a Hanning window.Fig. 4 Amplitude comparison between the baseline model (solid black), SXS:NSNS:0001 (dashed blue) and SXS:NSNS:0002 (dash-dot green) overlapped at the reference frequencies of the NR simulations and divided by the amplitude for SXS:BBH:0180.The insets contain the amplitudes for the pN (dotted maroon) and NRTidal (dash-dot-dot red) approximations, plotted up to the merger against the numerical BNS amplitude for the two cases.

Modeling the BNS Merger
Up to this point we have limited ourselves with modeling the tidal effects up to the merger, and confirmed that both the pN and NRTidal approximations hold reasonably well in comparison with the two numerical relativity simulations considered, although they were not included in the calibrations [30,32,33].Let us now take a closer look at the behavior of the analytical models for the tidal phase and amplitude as we approach the merger.Indeed, we shall see a noticeable mismatch with the numerical simulations, and will attempt to improve the model by performing our own fit only up to the merger.We will see that because of the smooth evolution of the phase and amplitude even after the moment of the merger, we can use the new coefficients obtained from our fit to push the model past the location where the stars touch.We succeed in extending the model for the phase beyond the merger, and devise a method to determine how far we can reach, where we end it with a taper and continue with the baseline model.We carry the amplitude up to the merger and terminate it with a Hann taper, to ensure a smooth and continuous transition into the post merger.

New Fit for the Tidal Phase
Let us proceed by first taking the difference between the numerical BNS and BBH phase for the two system considered.This is the true tidal correction up to 2rad.numerical error [44]) that we subsequently use to compare with the analytical models for the tidal phase.We plot in the left side of Fig. 5 the difference between the numerical tidal phase and the three analytical approximations for φ T,pN , φ T,F 1 and φ T,F 2 with the coefficients from [30,32,33].We start at a time about 1000 M before the merger, where the differences between the analytic and numerical phase become noticeable.We then extend the model past the merger, until the approximation exhibits a sharp increase and breaks down.We observe an increased difference between the analytical and numerical tidal phase near the merger, above the uncertainty of 2rad.In the right side of Fig. 5 we plot the true tidal phase (upper plot), as well as the scaled tidal deformability (lower plot).The true tidal phase scales with the tidal deformability before the merger, as expected, but after the merger is larger for the first (BNS1) system compared to the second (BNS2), appearing to scale inversely with the tidal deformability.The analytical phase, which depends linearly on the tidal deformability, doesn't provide an accurate estimation after the merger.Specifically, the postmerger tidal phase is underestimated for the first system and overestimated for the second.We assume that the Padé approximant from eq.( 35) is complex enough to model the smooth increase of the tidal phase at the merger, and proceed with performing new curve fits for the analytical tidal phase to the numerical tidal phase, scaled by the tidal deformability.We use as initial guesses the coefficients from the pN [30], F 1 [32] and F 2 [33] NRTidal model and stop the fit at the merger.Irrespective of the initial fitting coefficients, we obtain comparable sets of new coefficients for the polynomial modeling the tidal interaction of a BNS system with a given deformability.However, in contrast to the original fits, we do not find any longer a one-size-fits-all set of coefficients, namely our coefficients do depend on the equation of state considered.To select the best fit among the three variations of the analytical model, we calculate and compare the sum of squared residuals (SSR) for each numerical dataset.For BNS1, we select new parameters obtained by fitting to the numerical data, using the F 2 coefficients as the initial guess.For BNS2, we choose the parameters derived from fitting the numerical data with the F 1 coefficients as initial guess.Subsequently, we compare the residuals for the chosen coefficients between the two datasets and select the optimal fitting parameters for BNS1 as the new coefficients.
We give in Table 1 the average values of the new coefficients obtained, in comparison with the F 2 NRTidal coefficients.We need to alert readers that the large differences between our new fitting coefficients and the original ones are primarily due to our efforts to model the highly nonlinear tidal interactions during the merger, which is a step beyond the bounds of the analytical approximation The assumptions made in the original model that the coefficients were independent of the equation of state, are no longer applicable.We use the new fitting coefficients to reconstruct the analytical tidal phase, then we incorporate it into the baseline framework, thus obtaining a new analytical model for the BNS phase that extends beyond the merger.Following this, we develop a procedure that determines the termination point of our phase, where we apply a Heaviside function, allowing only the phase of the baseline model to continue from then on.First, we calculate the difference between our new analytical tidal phase and the true tidal phase.After that, we determine the time derivative of this difference.Next, we identify potential cutoff points at locations where this derivative switches sign.Lastly, we select the final cutoff point as the last point for which the deviation from the numerical phase is less than 3%, to optimize the fit with the true tidal phase.Here we terminate our analytical tidal phase with a Heaviside function and continue smoothly with the baseline phase, obtaining a complete phase representation.
We plot in Fig. 6 our newly modeled phase.In the main plot we show our new analytic fit for the BNS phase tailored to the two BNS systems, against the true tidal phase.Although we applied the curve fit only up to the merger, we see that it is able to accurately follow the numerical tidal phase beyond the merger.We include in the two insets a comparison between our reassembled tidal phase for the two BNS systems, with the Heaviside taper applied at the cutoff point, and the numerical BNS phase.Comparison between our analytic new fit to the tidal phase (dashed blue and dash-dot green) to the numerical tidal phase (solid and dotted black), as well as the original F 2 fit (solid orange), for the two BNS systems, up to the merger, then continued beyond the merger.The insets contain the reconstructed BNS phase for the two BNS systems, extended beyond the merger, as well as the new fit with the chosen coefficients to BNS2 (dash-dot-dot light green).The vertical yellow dotted line marks the merger, and the vertical brown dotted line delineates the end point where tapered.

New Fit for the Tidal Amplitude
We use a similar procedure to model the analytical amplitude at the merger and start by calculating the numerical tidal amplitude as the ratio between the numerical BNS and BBH amplitudes for the two system considered.This is the true tidal amplitude that we use next to compare with the analytical approximations A T,pN and A T,d .We attempt to improve the amplitude modeling by applying a curve fit for eqs.(38,39) to the true tidal amplitude up to the merger, using as initial guess the coefficients given in [33].Again, we obtain similar sets of new coefficients regardless of the model we start with, but depending on the value of the tidal deformability.We apply again the same procedure of calculating and comparing the sum of squared residuals (SSR) and select the optimal fitting parameters for BNS1 amplitude as the new coefficients.We give in Table 2 the values of the new coefficients resulting from the new tidal amplitude fit, in comparison with the A T,d NRTidal coefficients.With these coefficients, we proceed to reconstruct the new analytical amplitude.This time we cannot track the amplitude beyond the merger, due to its unphysical steep increase.While the model is precise enough to capture the sharp rise in amplitude before the merger, it is too basic to accurately follow its evolution beyond the peak point.To complete the waveform beyond the merger, we taper the amplitude at the merger with a Hann window that terminates to zero at the phase cutoff time.
We plot in Fig. 7 the comparison between the numerical BNS amplitude and our new amplitude with tides.The two insets display our fit to the numerical tidal amplitude, contrasted with the pN and NRTidal tidal amplitudes for the two BNS systems.

Complete Analytic BNS Waveforms
With these new fits obtained for the phase and amplitude, we build the complete analytical GW strain for the BNS merger using eq.( 28).Although we overlooked the amplitude complexity after the merger, we modeled the phase accurately up to the cutoff time.This allows us to recompute the orbital frequency Ω BNS by taking the time derivative of the phase φ BNS .Then we recalculate x BNS and follow with the radial separation r BNS , using eq.( 5).We also reevaluate the orbital and radial velocities and compare their evolution through the merger.Looking at their ratio we observe that even at the merger, where the radial velocity reaches its maximum, it represents at most 6% of the orbital velocity, after which it falls abruptly.In contrast, the orbital velocity keeps increasing even after the merger, reaching its peak soon thereafter, after which it starts decreasing as well.The domain between merger and cutoff time Fig. 7 Comparison between the two numerical BNS amplitudes (solid and dotted black) and our new fit for the amplitudes (dashed blue and dash-dot green) up to the merger, then tapered by a Hann function between merger and the cutoff time for the phase.The insets include our fit to the numerical tidal amplitude, compared with A T,d .The vertical yellow dotted line marks the merger.
is very short, extending up to 240M for BNS1 and 180M for BNS2.For the systems considered, this is between 2 msec and 3 msec, not enough to inform us on the fate of the remnant.Most likely, we reach the early stages of a rapidly-rotating, tidally deformed remnant.The GW170817 event indicates that a remnant with a total mass of 2.8M ⊙ will settle as a neutron star, the collapse to a black hole being less likely [4].Note that our model does not account for the mass ejected due to the tidal interactions during coalescence, or the dynamically ejected mass during the collision-induced shock of the neutron star crust when the stars come in contact.The matter expulsion at the contact interface is called shock ejecta, and represents a significant source of ejecta for systems with similar masses, as considered in this work.Additionally, our model excludes the neutrino-driven wind ejecta emitted by the HMNS remnant as it cools down.The different types of ejecta produced during the BNS evolution form a rotating disk surrounding the remnant.The total mass ejected by the GW170817 event was estimated to be about 0.04M ⊙ , and this is the value we consider in our work.
We plot in Fig. 8 the h + component of the strain for the baseline model, in comparison with the BNS strain for the two equations of state considered.While our model for the amplitude past the merger is simplified, the phase modeling remains accurate.In the left plot of Fig. 8 we show the evolution of the separation between the neutron stars for the two equations of state considered, continued beyond the merger, depicted by the red dotted circle, in comparison to the separation of the baseline BBH model.We see how the value of the tidal deformability and thus of the equation of state influence the early stages before and beyond the collision, revealing the effect of the matter interaction on the orbits.A larger tidal deformability (dashdot cyan) speeds up the merger of the stars, leading to a larger-radius remnant, whereas a smaller value for it (dotted blue) prolongs the inspiral and yields a denser remnant with a smaller radius.Fig. 8 Our analytical strain (left plot) and separation (right plot) for the baseline BBH (solid black), the Λ = 791 BNS1 (dotted blue) and the Λ = 1520 BNS2 (dashdot cyan) for the last ≈ 10 orbits extended to the early stages of the collision.The dotted red line marks where the neutron stars touch.

Conclusion
The gravitational waves emitted during a BNS merger offer insights into the behavior of matter under extreme conditions.To understand the impact of matter on the evolution of such systems and on their gravitational waves signature, we embarked upon an analytical modeling journey.Initially, we focused on a BBH system, employing the post-Newtonian formalism for the inspiral and the Backwards-one-Body model for the merger.By combining them, we established a baseline waveform, which we validated against numerical relativity simulations, to ensure a robust foundation for our next step.To incorporate the effects of tides, we introduced corrections to the phase and amplitude of the point-particle waveform by using the polynomial expressions proposed by the NRTidal model.We then verified the model's accuracy and efficiency through careful comparison with numerical relativity data for two equations of state.However, we encountered a mismatch around the merger when solely relying on the NRTidal model.To address this limitation, we lifted the restriction on the coefficients' independence from tidal deformability and recalibrated them with the numerical relativity predictions.By performing new fits to the numerical BNS data, we obtained updated values for the polynomial coefficients.Armed with these new coefficients, we reconstructed the tidal corrections and achieved improved fits for the phase and amplitude, successfully extending the phase modeling beyond the merger.To achieve a comprehensive phase representation, we devised a method to determine its extent, and applied a taper at the end, seamlessly continuing with the baseline model.Regarding the amplitude modeling, we have successfully carried it up to the merger, where we employed a Hann taper to ensure a smooth and continuous transition into the post-merger.We ended by reconstructing the complete analytical BNS strain and by investigating the tidal influence on the system's orbits around the BNS collision.
We developed RisingTides, a Python code guided by a user-friendly Jupyter Notebook, for the analytical modeling the tidal effects on the gravitational waves emitted during BNS inspiral and merger.We make our implementation available to the scientific community, to foster collaboration and facilitate further future investigations.
In our future work, we aim to enhance our model by incorporating a broader range of numerical relativity BNS simulations encompassing diverse equations of state.We will seek to uncover a consistent pattern for the dependence of the phase on the tidal deformability around the merger, ultimately leading to a universal set of polynomial coefficients.Furthermore, we will explore universal relations for BNS systems, which may offer deeper insights into their behavior and characteristics.Additionally, we are committed to refining the accuracy of the amplitude modeling beyond the merger, thus increasing the overall predictive power of our approach.Our research will contribute to a better understanding of BNS mergers and their gravitational wave signatures.Supplementary information.We release RisingTides, an open-source Python code for the analytical modeling of tidal effects on gravitational waves from BNS mergers.Our implementation can be accessed on GitHub, or through Zenodo.
Data Availability.The datasets generated and analyzed during the current study are openly available at github.com/mbabiuc/RisingTides,linked to the open-access Zenodo repository, at DOI: 10.5281/zenodo.8206261.

Fig.
Fig.1Comparison between the numerical and the analytic plus component of the strain, matched over a period around the initial frequency.The upper inset shows that the strain peaks are slightly displaced, but the merger is well modeled if we overlap them (dotted red curve).The dashed vertical yellow line marks the LR.The two lower insets contain the difference in phase (lower left) and amplitude (lower right) between our model and SXS:BBH:0180.Only expected small differences before the merger are visible.

Fig. 3
Fig.3Amplitude comparison between SXS:NSNS:0001 (dashed blue) and SXS:NSNS:0002 (dash-dot green) rescaled at the reference frequencies of the NR simulations.The inset contain a comparison of the numerical h + between the BNS and BBH data.

Fig. 5
Fig.5Left plots: Difference between the numerical tidal phase correction and the analytic pN and the three NRTidal models (pN, F 1 and F 2) for the tidal phase for the two BNS systems, close to the merger, and extended beyond it.Right plots: true and rescaled tidal phase for the two BNS systems, overlapped at the merger time.The vertical yellow dotted line marks the merger.

Fig. 6
Fig.6Comparison between our analytic new fit to the tidal phase (dashed blue and dash-dot green) to the numerical tidal phase (solid and dotted black), as well as the original F 2 fit (solid orange), for the two BNS systems, up to the merger, then continued beyond the merger.The insets contain the reconstructed BNS phase for the two BNS systems, extended beyond the merger, as well as the new fit with the chosen coefficients to BNS2 (dash-dot-dot light green).The vertical yellow dotted line marks the merger, and the vertical brown dotted line delineates the end point where tapered.

Table 1
New merger fit for tidal phase, dependent on the tidal deformability.

Table 2
New merger fit for tidal amplitude, dependent on the tidal deformability