# Kilohertz Gravitational Waves From Binary Neutron Star Mergers:
Numerical-relativity Informed Postmerger Model

Matteo Breschi1 Sebastiano Bernuzzi1 Kabir Chakravarti1 Alessandro
Camilletti2,3 Aviral Prakash5,6 Albino Perego2,3 1Theoretisch-Physikalisches
Institut, Friedrich-Schiller-Universität Jena, 07743, Jena, Germany
2Dipartimento di Fisica, Università di Trento, Via Sommarive 14, 38123 Trento,
Italy  3INFN-TIFPA, Trento Institute for Fundamental Physics and Applications,
via Sommarive 14, I-38123 Trento, Italy  5Institute for Gravitation & the
Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
6Department of Physics, The Pennsylvania State University, University Park, PA
16802, USA

###### Abstract

We present NRPMw, an analytical model of gravitational-waves from neutron star
merger remnants informed using 618 numerical relativity (NR) simulations.
NRPMw is designed in the frequency domain using a combination of complex
Gaussian wavelets. The wavelet’s parameters are calibrated to equations of
state (EOS) insensitive relations from NR data. The NR simulations are
computed with 21 EOS (7 of which are finite-temperature microphysical models,
and 3 of which contain quark phase transitions or hyperonic degrees of
freedom) and span total binary masses $M\in[2.4,3.4]~{}{\rm M_{\odot}}$, mass
ratios up to $q=2$, and (nonprecessing) dimensionless spins magnitudes up to
${0.2}$. The theoretical uncertainties of the EOS-insensitive relations are
incorporated in NRPMw using recalibration parameters that enhance the
flexibility and accuracy of the model. NRPMw is NR-faithful with fitting
factors ${\gtrsim}0.9$ computed on an independent validation set of $102$
simulations.

###### pacs:

04.25.D-, 04.30.Db, 95.30.Sf, 95.30.Lz, 97.60.Jd

##  I Introduction

This work is the first of a series of papers that present a faithful and
complete (inspiral-merger-postmerger) model for gravitational-wave (GW)
signals from binary neutron star (BNS) mergers, and its application to GW
analyses with the third-generation Einstein Telescope (ET) detector Hild  _et
al._ (2008, 2011); Hild (2012); Punturo  _et al._ (2010a); Maggiore  _et al._
(2020). Our model builds on a state-of-art effective-one-body (EOB) approach
for the inspiral-merger regime Bernuzzi  _et al._ (2015a); Nagar  _et al._
(2018); Akcay  _et al._ (2021); Gamba  _et al._ (2021a) and on its numerical
relativity (NR) completion for the remnant’s emission Bernuzzi  _et al._
(2015b); Breschi  _et al._ (2019). Prospects applications to ET GW
observations include: the precision measurement of the neutron star (NS) tidal
polarizability parameters Damour and Nagar (2009); Damour  _et al._ (2012),
the determination of the remnant’s black hole (BH) collapse Agathos  _et al._
(2020); Kashyap  _et al._ (2021), constraints on the extreme density equation
of state (EOS) Prakash  _et al._ (2021); Breschi  _et al._ (2022), and multi-
messenger observations Breschi  _et al._ (2021a). These case studies will be
further discussed in companion papers in the context of a Bayesian analysis
framework Breschi  _et al._ (tion). Here, we start presenting NRPMw, a new
analytical model for the postmerger (PM) emission from merger’s remnant, that
improves over our previous NRPM Breschi  _et al._ (2019).

The PM GW emission from a merger’s remnant is predicted to have a peak
luminosity at frequencies of few kilohertz, e.g. Shibata and Taniguchi (2006);
Hotokezaka  _et al._ (2013); Takami  _et al._ (2015); Bauswein and Stergioulas
(2015); Zappa  _et al._ (2018); Breschi  _et al._ (2019). This high-frequency
GW transient can be robustly computed by means of NR simulations and it is key
to directly probe the nature of the remnant in a (possibly multi-messenger)
BNS merger observation. A GW observation from a merger remnant is also a
promising probe for the nuclear EOS at extreme densities, e.g. Radice  _et
al._ (2017); Bauswein  _et al._ (2019); Breschi  _et al._ (2022); Most and
Raithel (2021); Prakash  _et al._ (2021). Kilohertz PM transients are unlikely
to be captured by current ground-base detectors Abbott  _et al._ (2018a), and
no PM signal was detected for GW170817 Abbott  _et al._ (2017a, b, 2018b,
2019). However, they are a main target for third-generation observatories Hild
_et al._ (2011); Punturo  _et al._ (2010b); Torres-Rivas  _et al._ (2019);
Maggiore  _et al._ (2020); Martynov  _et al._ (2019) and for finely-tuned
instruments Ackley  _et al._ (2020). In view of these considerations, it is
essential to develop accurate PM models for Bayesian GW analyses.

Models of PM GWs were presented in Refs. Hotokezaka  _et al._ (2013); Takami
_et al._ (2015); Clark  _et al._ (2016); Chatziioannou  _et al._ (2017a);
Easter  _et al._ (2019); Bose  _et al._ (2018); Breschi  _et al._ (2019);
Tsang  _et al._ (2019); Easter  _et al._ (2020); Soultanis  _et al._ (2021);
Wijngaarden  _et al._ (2022); Whittaker  _et al._ (2022). These templates are
phenomenological models that capture the main PM spectral features but do not
attempt to model the underlining remnant’s dynamics. The complex spectral
frequencies are either inferred from the observations or (in part) fixed by
EOS-insensitive (quasiuniversal) relations that connect the main spectral
features to the binary parameters. Depending on whether the quasiuniversal
relations are employed or not during the GW data inference (and for which
quantities), the templates might be used in fully-informed, partially informed
or agnostic approach. Importantly, all approaches require the quasiuniversal
relations to extract astrophysical constraints, either a priori or a
posteriori.

Most of the PM templates are built from a simple ansatz made of few damped
sinusoids in the time domain, eventually represented in the frequency domain.
Notable exceptions to a sinusoids basis are the models proposed in Refs. Clark
_et al._ (2016); Easter  _et al._ (2019) where reduced basis were constructed
directly from NR data. Clark  _et al._ (2016) used a principal component
analysis and ${\sim}50$ non-spinning simulations (12 of which unequal masses)
to demonstrate faithfulnesses ${\gtrsim}0.9$ on a subsample of the data.
Easter  _et al._ (2019) used a hierarchical model trained on 35 non-spinning,
equal-mass NR simulations to demonstrate fitting factors up to $0.98$ on the
training set. However, similar fitting factors can be achieved with
significantly less modeling efforts in agnostic approaches based on wavelets
or sinusoids basis Chatziioannou  _et al._ (2017a); Tsang  _et al._ (2019);
Easter  _et al._ (2020); Wijngaarden  _et al._ (2022). Moreover, the finite
precision of NR simulations introduces uncertainties that impact the
faithfulnesses at ${\sim}0.9$ level Easter  _et al._ (2019); Breschi  _et al._
(2019). Hence, simpler analytical templates appear favored over more complex
statistical models. The agnostic approach utilized in Ref. Chatziioannou  _et
al._ (2017a); Wijngaarden  _et al._ (2022) delivers, on average, larger
fitting factors to numerical data when compared to fully or partially informed
approaches, e.g. Breschi  _et al._ (2019); Tsang  _et al._ (2019); Easter  _et
al._ (2020); Soultanis  _et al._ (2021). This suggests that agnostic
approaches are able to detect PM signals at lower signal-to-noise ratio (SNR)
because informed models are not sufficiently accurate. The two approaches,
however, appear comparable at SNR relevant for astrophysical parameter
estimation, and they deliver comparable constraints on the EOS. We stress that
faithfulnesses calculations are often presented on validation datasets of
different sizes and a detailed comparison is difficult. For example, Easter
_et al._ (2020) found faithfulnesses between 0.91-0.97 on a sample of 9
simulations; Tsang  _et al._ (2019) found faithfulnesses between 0.4-0.95 on a
sample of 60 simulations, and Breschi  _et al._ (2019) between 0.4-0.95 on a
sample of about $150$ simulations. A main motivation for (partially) informed
approaches is the possibility to design inspiral-merger-postmerger templates
by consistently extending inspiral-merger templates. In Ref. Breschi  _et al._
(2019), we developed the first model of this kind by completing the EOB
framework of Refs. Bernuzzi  _et al._ (2015a); Nagar  _et al._ (2018) with the
NRPM PM model.

The new NRPMw is a PM frequency-domain template that aims at striking a
balance between fully-informed and agnostic approaches. It is constructed by
superposing few Gaussian, frequency-modulated wavelets whose parameters are
informed by new EOS-insensitive relations. The latter build on the largest
public databases of NR simulations available to date. The theoretical
uncertainties of the EOS-insensitive relations are incorporated in the model
using recalibration parameters that are determined during the inference.
Hence, NRPMw performs best in a partially informed inference. The
recalibration enhances the flexibility of the template and improves the
fitting factors to a level similar to agnostic templates. Data analysis
applications of NRPMw are presented in a companion paper.

The rest of this paper is structured as follows. In Sec. II, we discuss the PM
waveforms’ phenomenology predicted by state-of-art NR simulations. The
modeling choices used in NRPMw are presented in Sec. III. The quasiuniversal
relations calibrated for NRPMw are discussed in Sec. IV. In Sec. V we validate
the model against NR data by calculating its faithfulness on an independent
validation set. We summarize our findings and conclude in Sec. VI. Moreover,
we include several Appendices in order to extend the discussions on the
waveform modeling and on the calibration of EOS-insensitive relations.

#### Conventions –

We use geometric units $c=G=1$ or explicitely state units. Masses are
expressed in solar masses ${\rm M_{\odot}}$. The GW polarizations $h_{+}$ and
$h_{\times}$, plus and cross respectively, are decomposed in $(\ell,m)$
multipoles as

$h_{+}-{\rm
i}h_{\times}=D_{L}^{-1}\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}h_{\ell
m}(t)\,{}_{-2}Y_{\ell m}(\iota,\varphi),$ (1)

where $D_{L}$ is the luminosity distance, ${}_{-2}Y_{\ell m}$ are the $s=-2$
spin-weighted spherical harmonics and $\iota,\varphi$ are respectively the
polar and azimuthal angles that define the orientation of the binary with
respect to the observer. Each mode $h_{\ell m}(t)$ is decomposed in amplitude
$A_{\ell m}(t)$ and phase $\phi_{\ell m}(t)$, as

$h_{\ell m}(t)=A_{\ell m}(t)\,{\rm e}^{-{\rm i}\phi_{\ell m}(t)}\,,$ (2)

with a related GW frequency,

$\omega_{\ell m}(t)=2\pi f_{\ell m}(t)=\frac{{\rm d}}{{\rm d}t}{\phi_{\ell
m}}(t)\,.$ (3)

The moment of merger is defined as the time of the peak of $A_{22}(t)$, and
referred simply as merger when it cannot be confused with the
coalescence/merger process. If the multipolar indeces $(\ell,m)$ are omitted
from a multipolar quantity, we implicitly refer to the dominant $(2,2)$ mode.
Note that the time $t$ refers to the retarded time in the case of NR data. We
define the Fourier transform $h_{\ell m}(f)$ of each multipolar mode as

$h_{\ell m}(f)=\int_{-\infty}^{+\infty}h_{\ell m}(t)\,{\rm e}^{-2\pi{\rm
i}ft}\,{\rm d}t\,.$ (4)

Analogously to the time-domain case, the frequency series $h_{\ell m}(f)$ is
decomposed in amplitude $A_{\ell m}(f)$ and phase $\phi_{\ell m}(f)$.

The binary mass is $M=m_{1}+m_{2}$, where $m_{1,2}$ are the masses of the two
stars, the mass ratio $q=m_{1}/m_{2}\geq 1$, and the symmetric mass ratio
$\nu=m_{1}m_{2}/M^{2}$. We define the parameter $X=1-4\nu$. The dimensionless
spin vectors are denoted with ${\boldsymbol{\chi}}_{i}$ for $i=1,2$ and the
spin component aligned with the orbital angular momentum L are labeled as
$\chi_{i}={\boldsymbol{\chi}}_{i}\cdot\textbf{L}/|\textbf{L}|$. The effective
spin parameter ${\chi_{\rm eff}}$ is the mass-weighted aligned spin, i.e.

${\chi_{\rm eff}}=\frac{m_{1}\chi_{1}+m_{2}\chi_{2}}{M}\,.$ (5)

Moreover, the quadrupolar tidal deformability parameters are defined as
$\Lambda_{i}=({2}/{3})\,k_{2,i}\,C_{i}^{-5}$ for $i=1,2$, where $k_{2,i}$ and
$C_{i}$ are respectively the $\ell=2$ gravitoelectric Love number and the
compactness of the $i$-th NS. The tidal coupling constant is Damour and Nagar
(2010)

$\kappa^{\rm
T}_{2}=3\nu\,\left[\left(\frac{m_{1}}{M}\right)^{3}\Lambda_{1}+(1\leftrightarrow
2)\right]\,,$ (6)

that, similarly to the reduced tidal deformability ${\tilde{\Lambda}}$ Favata
(2014), parametrizes the leading-order tidal contribution to the binary
interaction potential.

##  II Waveform morphology

Figure 1: Representative examples of BNS PM waveforms. The plot shows the plus
polarization $h_{+}(t)$ of the time-domain $\ell=m=2$ waveform (solid line)
and the instantaneous GW frequency $\omega_{22}(t)$ (dashed line). The NR
simulations are from the CoRe database and computed in Refs. Bernuzzi  _et
al._ (2015a); Dietrich  _et al._ (2017a).

The PM waveform morphology and its connection to the remnant’s dynamics
predicted by simulations was discussed in various papers, see e.g. Ref.
(Shibata and Uryu, 2002; Stergioulas  _et al._ , 2011; Bauswein and Janka,
2012; Bauswein  _et al._ , 2012; Hotokezaka  _et al._ , 2013; Takami  _et al._
, 2014; Bauswein and Stergioulas, 2015; Bernuzzi  _et al._ , 2015b; Radice
_et al._ , 2016a; Lehner  _et al._ , 2016; Dietrich  _et al._ , 2017a, b;
Zappa  _et al._ , 2018; Bernuzzi  _et al._ , 2020). We review here the main
aspects that are relevant for the GW model proposed in the rest of the paper.
Figure 1 shows the PM signal in exemplary cases; the time axis is shifted to
the moment of merger.

A merger remnant is a massive, hot and rotating NS whose mass is usually
larger than the maximum mass sustained by a cold, isolated Tolmann-
Oppenheimer-Volkoff (TOV) NS. It can either collapse to a BH or settle to a
stable rotating NS on secular timescales. Gravitational collapse to BH takes
place as the remnant reaches densities comparable to the TOV’s maximum density
Perego  _et al._ (2021) since the remnant’s core is very slowly rotating
Kastaun  _et al._ (2017). The remnant of a very massive BNS can promptly
collapse after the moment of merger and crucially before the first bounce of
the two cores Radice  _et al._ (2020); Bernuzzi (2020) 111 This definition of
prompt collapse implies negligible shocked dynamical ejecta because the bulk
of this mass ejection comes precisely from the first core bounce Radice  _et
al._ (2018a). Since it directly connects to the main dynamical feature of the
merger process (shock and bounce) and to related observables, it is preferable
to other empirical definitions based on collapse time from the moment of
merger..

In the case of a equal mass BNS, the prompt collapse is described by empirical
relations relating the binary mass to the TOV maximum mass and compactness
proposed in Ref. Hotokezaka  _et al._ (2011); Bauswein  _et al._ (2013) and
refined in various works, e.g. Agathos  _et al._ (2020); Kashyap  _et al._
(2021). For very asymmetric BNS, the tidal disruption of the secondary drives
the gravitational collapse Bernuzzi  _et al._ (2020) and it is mainly
controlled by the incompressibility parameter of nuclear matter around the TOV
maximum density Perego  _et al._ (2021). While a robust prompt collapse
criterion is not known in these conditions Bernuzzi  _et al._ (2020); Bauswein
_et al._ (2020); Perego  _et al._ (2021), tidal disruption effects are of the
order of current EOS effects in the equal-mass criterion, at least for mass-
ratio $q\lesssim 1.4$ Kashyap  _et al._ (2021); Perego  _et al._ (2021).
Prompt collapse mergers have the largest GW luminosities (at merger) Zappa
_et al._ (2018) but the PM signal is the rapidly damped ringdown of the BH and
it is practically negligible for the sensitivities of current and next-
generation detectors. A prompt collapse signal is showed in the top panel of
Figure 1.

The evolution of a NS remnant is driven by an intense emission of GWs lasting
${\sim}10{-}20$ milliseconds (GW-driven phase) Bernuzzi  _et al._ (2016);
Zappa  _et al._ (2018). During this phase, the remnant either collapses to BH
(short-lived remnant) or settles to an approximately axisymmetric rotating NS
(long-lived remnant) 222A commonly used terminology for short-lived remnant is
hypermassive NS. This name is not appropriate for remnants since it refers to
cold equilibrium. See Radice  _et al._ (2018b); Bernuzzi (2020) and references
therein.. The GW-driven phase is associated to a luminous GW transient at
frequencies ${\sim}2-4\,$kHz Shibata and Uryu (2002); Stergioulas  _et al._
(2011); Bauswein and Janka (2012); Hotokezaka  _et al._ (2013); Takami  _et
al._ (2014); Bernuzzi  _et al._ (2015b). The spectrum of this transient is
rather complex but has robust and well-studied features at a few
characteristic frequencies. Most of the power is emitted in the $\ell=m=2$ GW
mode at a nearly constant frequency $\omega_{22}(t)\approx 2\pi f_{2}$.
Examples of $\ell=m=2$ waveforms for short- and long-lived remnants are shown
in the three bottom panels of Figure 1. The $f_{2}$ frequency is easily
extracted from simulation data and it was shown to correlate with various
binary quantities in a EOS-insensitive way, e.g. Bauswein and Stergioulas
(2015); Takami  _et al._ (2014); Bernuzzi  _et al._ (2015b); Breschi  _et al._
(2019).

We stress that the PM spectrum is not composed of a discrete set of
frequencies: the presence of broad peaks with typical full width at half
maximum (FWHM) of $300{-}600$ Hz is simply a consequence of the efficiency of
the emission process. Indeed, inspection of the time-domain waveform’s
instantaneous frequency (see Fig. 1) shows that $\omega_{22}(t)$ increases as
the remnant becomes more compact and has a steep acceleration towards
gravitational collapse 333Considering gauge-invariant energetics it is
possible to associate to the remnant a dynamical frequency $\Omega$ such that
$f_{2}=\Omega/\pi$ and analogously for other modes., see e.g. Fig. 1 of
Bernuzzi  _et al._ (2015b). Moreover, the instantaneous GW frequency has
modulations with frequencies $f_{0}\sim O(1~{}{\rm kHz})$ that are stronger
for remnants closer to collapse. These modulations are associated to the
violent radial bounces of the remnant’s core prior to collapse. Other robust
features of the spectrum are two secondary peaks at frequencies $f_{2\pm 0}$,
respectively at larger and smaller frequencies than $f_{2}$. These features
are associated to hydrodynamical modes in the remnant, e.g. Shibata and Uryu
(2000); Stergioulas  _et al._ (2011); Bernuzzi  _et al._ (2014a) and have been
was tentatively interpreted as nonlinear coupling between $f_{2}$ and $f_{0}$
Stergioulas  _et al._ (2011), in analogy to perturbations of rotating NS
Dimmelmeier  _et al._ (2006); Passamonti  _et al._ (2007); Baiotti  _et al._
(2009). The remnant’s signal from asymmetric binaries with mass ratio
$q\gtrsim 1.5$ carries the imprint of the tidal disruption during merger. An
example is shown in the bottom panel of Figure 1. The PM amplitude can be
significantly smaller than in the equal-mass cases and the peaks at
frequencies $f_{2\pm 0}$ are typically suppressed.

The evolution of a NS remnant beyond the GW-driven phase is highly uncertain
at present. It requires detailed simulations of viscous and nuclear processes
on timescales beyond hundreds of milliseconds, for example to quantify
precisely the mass accreting or outflowing the central object. NS remnants
after the GW-driven phase have an excess of both gravitational mass and
angular momentum when compared to equilibrium configuration with the
corresponding baryon mass Radice  _et al._ (2018b); Nedora  _et al._ (2021).
Possible mechanisms to shed (part of) this energy are long-term GW
instabilities Chandrasekhar (1970); Friedman and Schutz (1978) including one-
arm instabilities East  _et al._ (2016); Radice  _et al._ (2016a), that would
lead to potentially detectable, long GW transients at ${\lesssim}1\,$kHz.

The PM model presented in the next sections describes the GW transient during
the GW-driven phase and it builds on our previous work in Refs. Bernuzzi  _et
al._ (2015b); Breschi  _et al._ (2019). In particular, we devise new EOS-
insensitive relations based on the tidal coupling constant $\kappa^{\rm
T}_{2}$ and incorporate them in a partially informed model. We do not use
empirical relations for modeling prompt collapse and instead design a model
capable of inferring a generic collapse time from the observational data (but
see e.g. Ref. Agathos  _et al._ (2020) for an application of prompt collapse
quasiuniversal realtions in data analysis context). Similarly, we account for
the theoretical uncertainties of the EOS-insensitive using recalibration
parameters inferred from the data.

##  III NRPMw design

In order to develop an analytical NR-informed PM model for BNS mergers in the
frequency-domain, we first introduce a truncated complex Gaussian wavelet
${W}(t)$,

${W}(t;\alpha,\beta,\gamma,\tau)=\begin{cases}\begin{split}{\rm e}^{\alpha
t^{2}+\beta t+\gamma}\quad&{\rm if}\,\,t\in[0,\tau]\\\ 0\quad&{\rm
otherwise}\end{split}\end{cases}$ (7)

where $\alpha,\beta,\gamma\in\mathbb{C}$ are time-independent parameters and
the real interval $[0,\tau]$ defines the non-vanishing support of ${W}$. The
coefficients $\\{\alpha,\beta,\gamma\\}$ can be interpreted as follows:
$\Re{(\alpha)}$ and $\Re{(\beta)}$ determine respectively the concavity and
the initial slope of the time-domain wavelet amplitude; $\Im{(\alpha)}$ and
$\Im{(\beta)}$ define respectively the slope and the initial value of the
time-domain frequency evolution; $\gamma$ is an overall factor determining
initial amplitude and phase.

The frequency-domain wavelet ${W}(f)$ can be analytically computed from Eq.
(7) using Gaussian integrals,

${W}(f)=\frac{{\rm e}^{\gamma}}{2}\sqrt{\frac{\pi}{\alpha}}\,{\rm
e}^{-z^{2}}\,\left[{\rm erfi}{\left(z+\sqrt{\alpha}\tau\right)}-{\rm
erfi}{\left(z\right)}\right]\,,$ (8)

where $z(f)$ encodes the frequency dependency,

$z(f)=\frac{\beta-2\pi{\rm i}f}{2\sqrt{\alpha}}\,,$ (9)

and ${\rm erfi}{(z)}$ is the imaginary error function. For $\alpha=0$, Eq. (8)
is not defined and it is directly replaced by a Lorentzian function. Moreover,
a direct implementation of Eq. (8) can lead to floating point overflow fin a
certain portion of the parameter space. In these cases, we employ the
analytical approximations discussed in App. A. Furthermore, we introduce a
global time-shift ${\tau_{0}}$ in order to allow the wavelet to move on the
time axis. The time-shift ${\tau_{0}}$ changes the wavelet support to
$[{\tau_{0}},\tau+{\tau_{0}}]$ and it is easily implemented by a unitary
factor, i.e. ${W}(f;{\tau_{0}})={W}(f)\,{\rm e}^{-2\pi{\rm i}f{\tau_{0}}}$.

The wavelet is the basic component of NRPMw. In the following paragraphs we
describe how different wavelets are combined based on the universal features
of the PM signal that are identified by characteristic times (“nodes”, Sec.
III.1). Then, we discuss the modeling of sub-dominant frequencies as
additional wavelet modulations in Sec. III.2. The basic construction of the
dominant $\ell=m=2$ mode is discussed Sec. III.3 and the modeling of higher-
order modes in Sec. III.4.

###  III.1 Nodal points

The time-domain strain has universal characteristic features at specific
times, as pointed out in Ref. Breschi  _et al._ (2019) (see also Fig. 2). We
call these times nodal points and indicate them as $\\{t_{i}\\}$ for
$i=0,1,2,3$. Nodal points are identified as stationary points of the strain’s
amplitude, that we indicate as $\\{A_{i}\\}$. Differently from Ref. Breschi
_et al._ (2019), we assume $t_{i+1}-t_{i}$ to be constant, for $i=0,1,2$.
Hence, the nodal points can be reduced to two independent parameters: the
moment $t_{0}$ of the first amplitude’s minimum after merger, and a
characteristic time-scale ${T_{0}}$ that is computed as the difference
$t_{3}-t_{1}$. The time-scale ${T_{0}}$ defines the subdominant frequency
$f_{0}\simeq{T_{0}}^{-1}$ that characterizes the modulations of the PM signal.
A further time-domain node is introduced for the time of the remnant collapse
$t_{\rm coll}$. Differently from Breschi  _et al._ (2019), here we do not
introduce $t_{4}$.

###  III.2 Amplitude and frequency modulations

Ampitude and frequency modulations (AMs, FMs) are prominent features of the PM
spectrum, as discussed in Sec. II. NR simulations show that the main GW
modulations are given in the $m=0$ channel, and are associated to the quasi-
radial density oscillations of the remnant Bernuzzi  _et al._ (2008). We
associate this mode to the fundamental frequency $f_{0}$ and, for the modeling
of the $(2,2)$ mode, we consider only the modulation couplings between $f_{2}$
and $f_{0}$ 444Our simulations indicate that the couplings between $(2,1)$ and
$(3,3)$ modes can also be relevant for unequal-mass BNS.. Moreover, we neglect
possible frequency evolution of the subdominant mode $f_{0}$, i.e. this
frequency component is assumed to be constant in time. Modulation effects
appear after the collision of the NS cores, for $t>t_{0}$, when the remnant is
strongly deformed and dynamically unstable.

AMs can be easily taken into account by employing a combination of wavelets.
Labeling the amplitude-modulated wavelet as ${\breve{W}}$, we can write

$\begin{split}{\breve{W}}(t)&={W}(t)\left[1+\Delta_{\rm
am}\sin\left(\Omega_{\rm am}t+\phi_{\rm am}\right)\right]\\\
&={W}(t)-\frac{{\rm i}\Delta_{\rm am}}{2}\sum_{k=\pm 1}k\,{W}(t)\,{\rm
e}^{{\rm i}k(\Omega_{\rm am}t+\phi_{\rm am})}\,,\end{split}$ (10)

where $\Delta_{\rm am}$ defines the magnitude, $\Omega_{\rm am}$ the
modulation frequency and $\phi_{\rm am}$ the initial phase of the AMs. Eq.
(10) can be transformed in the Fourier space obtaining

${\breve{W}}(f)={W}(f)-\frac{{\rm i}\Delta_{\rm am}}{2}\sum_{k=\pm
1}k\,{W}^{(k)}(f)\,,$ (11)

where

${W}^{(k)}(f)={W}(f;\alpha,\beta+{\rm i}k\Omega_{\rm am},\gamma+{\rm
i}k\phi_{\rm am},\tau)\,.$ (12)

Eq. (11) shows explicitly that an amplitude-modulated wavelet ${\breve{W}}$
can be easily written in terms of the Gaussian wavelets ${W}$ and it
introduces two subdominant contributions in the Fourier domain that are
displaced with respect to the primary peak of $\pm\Omega_{\rm am}$.

FMs affect the phase evolution of the time-domain wavelet. We implement a FM
wavelet ${\tilde{W}}$ defining the frequency evolution as

$\omega_{\tilde{W}}(t)=\omega_{W}(t)-\Delta_{\rm fm}{\rm e}^{-\Gamma_{\rm
fm}t}\sin(\Omega_{\rm fm}t+\phi_{\rm fm})\,,$ (13)

where $\omega_{\tilde{W}}$ is the instantaneous frequency of the frequency-
modulated wavelet ${\tilde{W}}$, $\omega_{W}$ is the instantaneous frequency
of the Gaussian wavelet ${W}$, and $\Delta_{\rm fm},~{}\Gamma_{\rm
fm},~{}\Omega_{\rm fm},~{}\phi_{\rm fm}\in\mathbb{R}$ are the parameters that
define the FM, i.e. $\Delta_{\rm fm}$ is the initial frequency displacement,
$\Gamma_{\rm fm}$ the inverse damping time, $\Omega_{\rm fm}$ the modulation
frequency and $\phi_{\rm fm}$ the initial phase. Using Taylor expansion, the
frequency-modulated wavelet ${\tilde{W}}$ can be rewritten in terms of the
frequency-domain Gaussian wavelet ${W}$. A detailed discussion on the analytic
form of ${\tilde{W}}(f)$ is provided in App. B. Note that, differently from
the AMs shown in Eq. (10), the FM contribution presented in Eq. (13) includes
damped behavior, i.e. $\Gamma_{\rm fm}\neq 0$ a priori. This term is needed to
properly characterize the different time-scales of the PM frequency components
$f_{2}$ and $f_{0}$.

Combining the definitions of ${\breve{W}}$, Eq. (11), and ${\tilde{W}}$ (see
Eq. (38)), it is possible to write a general modulated Gaussian wavelet,
labeled as ${\widetilde{W}}$. We consider AMs over the interval
$[t_{0},t_{3}]$ and FM for $t>t_{0}$. We fix the modulation frequencies to
$\Omega_{\rm am}=\Omega_{\rm fm}=2\pi f_{0}$. Then, the AM magnitude
$\Delta_{\rm am}$ and phase $\phi_{\rm am}$ are fixed by the values of the GW
amplitudes at the nodal points nodal points, i.e. $\\{t_{i},A_{i}\\}$ for
$i=1,2,3$. The FM inverse damping time $\Gamma_{\rm fm}$ is assumed to be
identically zero for $t<t_{1}$; then, it is fixed to a constant positive value
calibrated on NR data (see Sec. IV). Furthermore, NR simulations show that AMs
and FMs approximately fluctuate in opposite directions Breschi  _et al._
(2019); i.e. amplitude maxima occur at frequency minima and viceversa. The FM
phase $\phi_{\rm fm}$ is fixed in order to satisfy this requirement.

###  III.3 Wavelet combination

Figure 2: Exemplary case showing the morphology of NRPMw model. Different
wavelet components are reported with different colors: ${W}_{\rm fus}$ in
blue, ${\tilde{W}}_{\rm bnc}$ in orange, ${\widetilde{W}}_{\rm pul}$ in green,
and ${\tilde{W}}_{\rm peak}$ in purple. The top panel shows the time-domain
components and the overall GW amplitude $A(t)$ (black line) highlighting the
characteristic times with vertical lines, i.e. the time of the merger $t_{\rm
mrg}$, the nodal points $t_{i}$ for $i=0,1,2,3$ and the time of collapse
$t_{\rm coll}$. The bottom panel shows the Fourier spectra of each component,
the overall $h_{22}$ spectrum (black line) and the characteristic PM
frequencies (vertical lines), i.e. the merger frequency $f_{\rm mrg}$, the PM
peak $f_{2}$ and the subdominant couplings $f_{2\pm 0}=f_{2}\pm f_{0}$.

The NRPMw model is constructed by describing each part of the PM signal
between different nodal points with a modulated wavelet component. The overall
strain $h_{22}$ is computed summing all the contributions. The use of wavelets
allow us to assign a clear interpretation of each parameter employed in the
model. The combination of different wavelets can capture rather complex signal
morphologies.

In NRPMw, the physical quantities (times, amplitudes and frequencies) are
estimated using quasiuniversal relations calibrated on NR simulations (see
Sec. IV). This allows us to design a fully informed model that can connect the
signal’s morphology to the intrinsic parameters of the BNS system (masses,
spins and tidal parameters). However, some wavelet parameters could be let
unconstrained and directly inferred from observational data Easter  _et al._
(2020) or they could be reconstructed with regression methods directly from NR
simulations Easter  _et al._ (2019).

The time-domain $\ell=m=2$ mode is modeled employing a combination of four
different wavelet components,

$h(t)\approx{W}_{\rm fus}(t)+{\tilde{W}}_{\rm bnc}(t)+{\widetilde{W}}_{\rm
pul}(t)+{\tilde{W}}_{\rm peak}(t)\,,$ (14)

assuming continuity in amplitude and phase (except for a phase-shift
$\phi_{\rm PM}$, see later) for the time-domain counterpart. Detailed
expressions are given in App. C. The combination of wavelets includes the
following terms that are shown in color in Fig. 2:

  1. 1.

${W}_{\rm fus}$ describes the signal after merger and up to $t_{0}$,
corresponding to the fusion of the NS cores. The wavelet has an initial
frequency $f_{\rm mrg}$ and non-vanishing frequency drift that can be positive
or negative depending on the properties of the binary;

  2. 2.

${\tilde{W}}_{\rm bnc}$ describes the signal corresponding to the bounce after
the collision of the cores. The phase here has a discontinuity $\phi_{\rm PM}$
at $t_{0}$. Moreover, for $t>t_{0}$, all wavelets include FMs with the
subdominant frequency $f_{0}$;

  3. 3.

${\widetilde{W}}_{\rm pul}$ describes the emission up to $t_{3}$ during which
the remnant is typically highly dynamical. Since the largest amount of the GW
luminosity is emitted at times $\lesssim 5~{}{\rm ms}$ Zappa  _et al._ (2018),
this component also includes AMs with the subdominant frequency $f_{0}$;

  4. 4.

${\tilde{W}}_{\rm peak}$ describes the signal after the luminosity peak by a
damped sinusoidal with initial frequency $f_{2}$, a frequency evolution
parametrized by the drift $\alpha_{\rm peak}$ (also referred as
${\Im{(\alpha_{\rm peak})}}$ in App. C). This component characterizes the
dominant Fourier peak and it lasts until the time of collapse $t_{\rm coll}$.

Additionally, the GWs emitted by the collapse and BH ringdown can be modeled
as a fifth term in Eq. (14), ${W}_{\rm coll}$ (see App. C for a detailed
discussion). Knowing the properties of the final BH, this component could be
modeled with the quasi-normal modes of the remnant Berti  _et al._ (2009);
Berti and Klein (2014). For simplicity, however, we set here ${W}_{\rm
coll}=0$.

Figure 2 shows an example of the discussed contributions in time- and
frequency-domain, with the different terms appearing in Eq. 14 shown in
different colors. The overall spectrum shows the typical PM patterns: a
dominant quasi-Lorentzian peak, a weaker peak at lower frequencies
corresponding to the merger dynamics and sub-dominant peaks due to AMs and
FMs. The superposition of the wavelet components generates several local
minima and maxima in the overall $h_{22}$ spectrum. Moreover, the destructive
interference of the wavelets originates a local minimum typically located
between $f_{\rm mrg}$ and $f_{2}$. This feature is also generally observed in
BNS PM spectra extracted from NR simulations. Moreover, the sharp cut at
$t_{\rm coll}$ in time-domain waveform originates the ringing effects observed
in the $h_{22}$ spectrum 555This can be easily seen performing the convolution
product of a sinusoidal wavelet with a Heaviside function.. The further
inclusion of $W_{\rm coll}$ will mitigate this effect, yielding to a smoother
waveform representation.

Overall, the model is characterized by 17 parameters, that are the
characteristic frequencies, amplitudes, times and phases that define
instantaneous GW amplitude and frequency (see App. C). Most of these
quantities can be related to the binary properties using NR information.

###  III.4 Higher-order modes

NR simulations show prominent coupling effects in higher mode (HM) terms of
BNS PM transients, similarly to what we discussed for the dominant $(2,2)$
mode. Also for this reason, the power of HM contributions in BNS PM radiation
is considerably larger compared with the pre-merger dynamics Calderón Bustillo
_et al._ (2021); Bernuzzi  _et al._ (2020). For typical BNS systems, these
contributions cover a relatively broad spectrum, roughly from ${\sim}500~{}$Hz
to 5–7 kHz.

In general, HM contributions can be modeled as a combination of wavelets with
different frequencies imposing continuity in amplitude and phase. For $m\neq
0$, the characteristic peak frequencies of HMs can be approximated using the
quadrupolar term employing the multipolar scaling, i.e. $f_{\ell
m}\simeq({m}/{2})\,f_{2}$. However, the hierarchy of frequency couplings is
not fully resolved. A detailed analysis of these subdominant features might
require better resolved simulations to robustly identify the trend in the
spectra. We remand the inclusion of HM PM characteristic properties to a
future study.

##  IV NR calibration

The NRPMw model has 17 parameters, i.e.

$\begin{split}\boldsymbol{\theta}_{\rm PM}=\\{&\phi_{\rm PM},\phi_{\rm
fm},t_{0},t_{\rm coll},\\\ &A_{\rm mrg},A_{0},A_{1},A_{2},A_{3},\\\ &f_{\rm
mrg},f_{2},f_{0},\Delta_{\rm fm},\Gamma_{\rm fm},\\\ &{\Re{(\beta_{\rm
peak})}},{\Im{(\alpha_{\rm fus})}},\alpha_{\rm peak}\\}\,,\end{split}$ (15)

that can be mapped to the binary parameters,

$\boldsymbol{\theta}_{\rm
bin}=\left\\{m_{1},m_{2},\Lambda_{1},\Lambda_{2},\chi_{1},\chi_{2}\right\\}\,,$
(16)

using NR simulations. We chose to map only a subset of
$\boldsymbol{\theta}_{\rm PM}$ and let some other parameters to be determined
by the inference or any other minimization procedure with given data. In
particular, we map the following 13 parameters

$\begin{split}{\boldsymbol{\theta}_{\rm fit}}=\\{&A_{\rm
mrg},A_{0},A_{1},A_{2},A_{3},f_{\rm mrg},f_{2},f_{0},t_{0},\\\
&{\Re{(\beta_{\rm peak})}},{\Im{(\alpha_{\rm fus})}},\Delta_{\rm
fm},\Gamma_{\rm fm}\\}\,,\end{split}$ (17)

we fix $\phi_{\rm fm}$ by the AMs and the FMs as discussed in Sec. III.2, and
we leave three additional degrees of freedom,

${\boldsymbol{\theta}_{\rm free}}=\left\\{\phi_{\rm PM},t_{\rm
coll},\alpha_{\rm peak}\right\\}\,.$ (18)

This choice is motivated by the fact that these three parameters cannot be
robustly mapped using NR data. The PM phase $\phi_{\rm PM}$ shows a strong
dependence on the simulation’s grid resolutions and on the physical models,
e.g. Bernuzzi  _et al._ (2012); Dietrich  _et al._ (2018a); Bernuzzi and
Dietrich (2016). The time of collapse $t_{\rm coll}$ is difficult to robustly
determine from simulations due its dependence on grid resolution Breschi  _et
al._ (2019); moreover, it strongly depends on the properties of the nuclear
EOS and might be biased by the relatively small EOS set available Radice  _et
al._ (2017); Prakash  _et al._ (2021); Fujimoto  _et al._ (2022). The
frequency drift $\alpha_{\rm peak}$ is also connected to the collapse dynamics
and, as such, it can be affected by various processes, especially in long-
lived remnants. For example, we discuss in App. D the dependency of
$\alpha_{\rm peak}$ on the turbulent viscosity in a subset of simulations.

The calibration set of binaries includes the public available non-precessing
NR simulations of the CoRe Dietrich  _et al._ (2018b); cor  and the SACRA
Kawaguchi  _et al._ (2018); Kiuchi  _et al._ (2017, 2019) databases, plus
additional data from simulations of Ref. Bernuzzi  _et al._ (2020); Prakash
_et al._ (2021); Camilletti  _et al._ (2022) with the BLh and BLQ EOS. The
CoRe database includes data computed with two different NR codes, BAM Brügmann
_et al._ (2008); Thierfelder  _et al._ (2011) and THC Radice and Rezzolla
(2012) and simulate microphysics, neutrino transport (with various schemes)
and turbulent viscosity. The final dataset is composed by 618 simulations and
it includes 190 different binary configurations computed with three
independent NR codes and 21 different EOSs. The finite temperature,
composition-dependent EOSs are BHB$\Lambda\phi$ Banik  _et al._ (2014), BLh
Bombaci and Logoteta (2018); Logoteta  _et al._ (2020), BLQ Bombaci and
Logoteta (2018); Logoteta  _et al._ (2020); Prakash  _et al._ (2021), HS(DD2)
(Typel  _et al._ , 2010; Hempel  _et al._ , 2012, DD2 hereafter), LS220
Lattimer and Swesty (1991), SFHo Steiner  _et al._ (2013), SRO(SLy) (Schneider
_et al._ , 2017, SLy hereafter); the EOSs in piecewise polytropic forms are:
ALF2 Alford  _et al._ (2005), ENG Engvik  _et al._ (1996), MPA Müther  _et
al._ (1987), MS1 Müller and Serot (1996), MS1b Müller and Serot (1996), SLy
Douchin and Haensel (2001), 2B, 2H, 15H, 125H, B, H, H4, HB from Ref. Lackey
_et al._ (2006); Read  _et al._ (2009) and the $\Gamma{=}2$ ideal gas EOS. We
remark that ALF2 and BLQ include a phase transition to deconfined quark
matter, and BHB$\Lambda\phi$ takes into account the appearance of hyperons at
high densities. The intrinsic parameters of the data cover the ranges
$M\in[2.4,3.4]~{}{\rm M_{\odot}}$, $q\in[1,2]$, $\kappa^{\rm
T}_{2}\in[22,458]$ and ${\chi_{\rm eff}}\in[-0.14,+0.22]$. Among the
considered dataset, 80 simulations (${\sim}13\%$ of the sample) resulted in
prompt collapse and ${\sim}40\%$ of the total data is composed by equal-mass
non-spinning binaries. We include all available resolutions for every binary
configuration and we treat each point as an independent measure in order to
improve the characterization of NR uncertainties. The quasiuniversal relations
presented in this work extend those in Ref. Breschi  _et al._ (2019) including
effects of large mass ratios, i.e. $q>1.5$, and aligned spins with
$|{\chi_{\rm eff}}|\lesssim 0.2$. In App. F, we present a recalibration of the
quasiuniversal relations between the PM peak frequency $f_{2}$ and the NS
radius that is not used in NRPMw but often employed in GW inference.

The mapping between binary and NRPMw parameters is performed on the mass-
rescaled PM parameters using a factorized fitting function (for any quantity
$Q$),

$Q^{\rm fit}=a_{0}\,\,Q^{\rm M}(X)\,\,Q^{\rm S}(\hat{S},X)\,\,Q^{\rm
T}(\kappa^{\rm T}_{2},X)\,,$ (19)

where $Q^{\rm M}=1+a_{1}^{\rm M}X$ includes the mass ratio contributions in
terms of the $X=1-4\nu$ parameter; $Q^{\rm S}=1+p_{1}^{\rm S}\hat{S}$ takes
into account spin corrections in terms of the spin parameter Nagar  _et al._
(2019)

$\hat{S}=\left(\frac{m_{1}}{M}\right)^{2}\chi_{1}+\left(\frac{m_{2}}{M}\right)^{2}\chi_{2}\,,$
(20)

and $p_{1}^{\rm S}=a_{1}^{\rm S}(1+b_{1}^{\rm S}X)$. The term

$Q^{\rm T}=\frac{1+p_{1}^{\rm T}\,\kappa^{\rm T}_{2}+p_{2}^{\rm
T}\,{\kappa^{\rm T}_{2}}^{2}}{1+p_{3}^{\rm T}\,\kappa^{\rm T}_{2}+p_{4}^{\rm
T}\,{\kappa^{\rm T}_{2}}^{2}}\,,$ (21)

takes into account tidal effects in terms of $\kappa^{\rm T}_{2}$ and with
$p_{i}^{\rm T}=a_{i}^{\rm T}(1+b_{i}^{\rm T}X)$. The coefficients
$\\{a_{i}^{\cdot},b_{i}^{\cdot}\\}$ are determined fitting the NR data. We
note that the choice of the fitting function in Eq. (19) might be not unique
nor optimal; we have experimented with few functions and found Eq. (19)
sufficiently simple, general and accurate for our purposes. The choice of a
rational function for $Q^{\rm T}(\kappa^{\rm T}_{2})$ is instead motivated by
previous works Bernuzzi  _et al._ (2014b, 2015b); Zappa  _et al._ (2019);
Breschi  _et al._ (2019). Finally, we stress the importance of using mass-
rescaled quantities in quasiuniversal relations Bernuzzi  _et al._ (2015b);
Breschi  _et al._ (2019); App. F demonstrates that factorizing the (trivial)
binary mass scale is key to obtain EOS-insensitive relations.

The fitting is performed using a least squared method. Denoting by $Q^{\rm
NR}_{i}$ any NR quantity of interest extracted from the $i$-th NR simulation
and $Q^{\rm fit}_{i}$ its fit, we define the relative residual of the $i$-th
NR simulation,

$r_{i}=\frac{Q^{\rm fit}_{i}-Q^{\rm NR}_{i}}{Q^{\rm fit}_{i}}\,,$ (22)

and minimize $\chi^{2}=\sum_{i}r_{i}^{2}$. For each calibrated PM parameter,
Table 1 reports the calibrated coefficients and the associated relative error,
defined as the standard deviation of the relative residuals, i.e. $\sqrt{{\rm
Var}(r_{i})}$. For later purposes (see Sec. V.1), we report in Table 1 the
Kullback–Leibler divergence $D_{\text{KL}}$ between the distribution of the
residuals $r_{i}$ and a normal distribution with zero mean and variance ${\rm
Var}(r_{i})$. This quantity allows us to verify the Gaussian character of the
residuals.

In the following, we discuss the fit results, i.e. empirical relations for the
merger properties (Sec. IV.1), for the characteristic PM frequencies and
amplitudes (respectively Sec. IV.2 and Sec. IV.3) and for the late-time
properties (Sec. IV.4).

Figure 3: Quasi-universal relation for the PM peak frequency $f_{2}$ as
function of the tidal polarizability $\kappa^{\rm T}_{2}$. Top panel:
calibrated relations (black lines) compared to NR data (colored dots)
extracted from the CoRe and the SACRA databases. Each color corresponds to a
different EOS. NR medians and error-bars are reported averaging over different
numerical resolutions (when available) for the same binary configuration.
Bottom panel: Relative residuals between the calibrated relation and the NR
data validation set. The gray areas show the 50% (dark) and 90% (light)
credible regions of the residuals.  Table 1:  Summary of the quasiuniversal
relations for the PM parameters ${\boldsymbol{\theta}_{\rm fit}}$ as functions
of the inspiral parameters $\boldsymbol{\theta}_{\rm bin}$. The first column
report the quantity of interest and the second column shows the range spanned
by the available NR data. From the third to the fourteenth column, we report
the calibrated coefficients of the quasiuniversal relations. The last three
columns show respectively the $\chi^{2}$, the relative standard deviation of
the fit and the Kullback–Leibler divergence $D_{\text{KL}}$ between the
recovered residuals and a normal distribution.

$Q^{\rm fit}$ | Range | $a_{0}$ | $a_{1}^{\rm M}$ | $a_{1}^{\rm S}$ | $b_{1}^{\rm S}$ | $a_{1}^{\rm T}$ | $a_{2}^{\rm T}$ | $a_{3}^{\rm T}$ | $a_{4}^{\rm T}$ | $b_{1}^{\rm T}$ | $b_{2}^{\rm T}$ | $b_{3}^{\rm T}$ | $b_{4}^{\rm T}$ | $\chi^{2}$ | Error | $D_{\text{KL}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$A_{\rm mrg}/M$ | [ 0.159 , 0.313 ] | $0.3948$ | $-1.133$ | $-0.02992$ | $-2.593$ | $0.03902$ | $5.1846{\times 10^{{-5}}}$ | $0.06033$ | $1.380{\times 10^{{-4}}}$ | $10.41$ | $54.51$ | $10.83$ | $54.54$ | $0.189$ | $1.8\%$ | $0.26$
$A_{0}/M$ | [ $2.04{\times 10^{{-4}}}$ , $0.0699$ ] | $0.02356$ | $0$ | $1.077$ | $260.4$ | $-1.318{\times 10^{{-3}}}$ | $0$ | $0$ | $0$ | $-4.314$ | $0$ | $0$ | $0$ | $267$ | $66\%$ | $0.14$
$A_{1}/M$ | [ $0.0262$ , $0.238$] | $-0.05641$ | $-5$ | $-1.135$ | $146.8$ | $-0.8343$ | $3.882{\times 10^{{-4}}}$ | $0.2464$ | $0$ | $-5$ | $0$ | $0$ | $0$ | $13.2$ | $15\%$ | $0.18$
$A_{2}/M$ | [ $4.76{\times 10^{{-4}}}$ ,$0.175$ ] | $0.1667$ | $-5.135$ | $-3.796$ | $-28.47$ | $0$ | $0$ | $5.774{\times 10^{{-3}}}$ | $0$ | $0$ | $0$ | $4.027{\times 10^{{-8}}}$ | $0$ | $78.2$ | $38\%$ | $0.077$
$A_{3}/M$ | [ $5{\times 10^{{-3}}}$ , $2.04{\times 10^{{-1}}}$ ] | $0.1662$ | $0.1072$ | $-2.046$ | $-45.06$ | $-7.06{\times 10^{{-5}}}$ | $0$ | $1.354{\times 10^{{-3}}}$ | $0$ | $-1423$ | $0$ | $284.7$ | $0$ | $38.4$ | $26\%$ | $0.044$
$Mf_{\rm mrg}/\nu$ | [ $0.0554$, $0.141$] | $0.2276$ | $0.9233$ | $0.5938$ | $-1.994$ | $0.03445$ | $5.58{\times 10^{{-6}}}$ | $0.08405$ | $1.133{\times 10^{{-4}}}$ | $13.83$ | $517.4$ | $12.75$ | $139.8$ | $0.431$ | $2.6\%$ | $0.45$
$Mf_{2}$ | [ $0.0216$ , $0.0512$] | $0.0881$ | $22.81$ | $0.2925$ | $25.0$ | $0.007023$ | $-1.782{\times 10^{{-6}}}$ | $0.02587$ | $6.58{\times 10^{{-6}}}$ | $5.428$ | $0$ | $39.29$ | $0$ | $0.814$ | $3.9\%$ | $0.029$
$Mf_{0}$ | [ $1.86{\times 10^{{-3}}}$, $0.0441$] | $0.02734$ | $19.32$ | $-1.857$ | $-75.77$ | $-2.967{\times 10^{{-3}}}$ | $8.484{\times 10^{{-6}}}$ | $8.584{\times 10^{{-3}}}$ | $0$ | $20.49$ | $21.5$ | $10.47$ | $0$ | $107$ | $45\%$ | $0.10$
$M/t_{0}$ | [ $8.09{\times 10^{{-3}}}$, $0.0288$] | $0.03265$ | $0.2994$ | $-0.2329$ | $4.768$ | $3.584{\times 10^{{-3}}}$ | $0$ | $0.01053$ | $0$ | $-11.96$ | $0$ | $-3.22$ | $0.0$ | $5.11$ | $9.2\%$ | $0.72$
$M\,{\Re{(\beta_{\rm peak})}}$ | [ $6{\times 10^{{-4}}}$ ,$6.58{\times 10^{{-3}}}$ ] | $0.1912$ | $4.074$ | $-1.573$ | $100$ | $0.05884$ | $0$ | $3.896$ | $0$ | $-5.293$ | $0$ | $0$ | $0$ | $37.7$ | $27\%$ | $0.042$
$M^{2}\,\Im(\alpha_{\rm fus})/\nu$ | [ $-3.3{\times 10^{{-4}}}$, $5{\times 10^{{-3}}}$ ] | $0.003721$ | $-1.799$ | $0.3555$ | $-7.167$ | $0.0139$ | $-2.425{\times 10^{{-5}}}$ | $0.05883$ | $1.882{\times 10^{{-4}}}$ | $-28.64$ | $-36.18$ | $19.53$ | $7.089$ | $346$ | $75\%$ | $3.7$
$M\,\Delta_{\rm fm}$ | [ $1.5{\times 10^{{-4}}}$, $0.0423$ ] | $0.05139$ | $0.4944$ | $-3.734$ | $-145$ | $-6.25{\times 10^{{-3}}}$ | $1.728{\times 10^{{-5}}}$ | $0.01944$ | $0$ | $-7.936$ | $1.882$ | $100$ | $0$ | $278$ | $74\%$ | $0.041$
$M\,\Gamma_{\rm fm}$ | [ $0$, $0.05$] | $0.1637$ | $209.3$ | $-0.2997$ | $24.5$ | $0.02195$ | $0$ | $0.3528$ | $0$ | $-0.5111$ | $0.$ | $74.72$ | $0$ | $326$ | $98\%$ | $1.05$

###  IV.1 Merger properties

Among all the quantities of interest, the amplitude and the frequency at
merger, respectively $A_{\rm mrg}$ and $f_{\rm mrg}$, are properties that can
be extracted with good accuracy from NR data Bernuzzi  _et al._ (2014b);
Breschi  _et al._ (2019). Our new relations have 1-$\sigma$ errors smaller
than $3\%$, as shown in Tab. 1. These relations are constructed to match the
binary black hole (BBH) values for $\kappa^{\rm T}_{2}\to 0$; the limiting
values are taken from the EOB model of Ref. Nagar  _et al._ (2019).

The slope parameter ${\Im{(\alpha_{\rm fus})}}$ characterizes the derivative
of the GW frequency immediately after merger, i.e. ${\Im{(\alpha_{\rm
fus})}}\propto\left({{\rm d}f}/{{\rm d}t}\right)_{\rm mrg}$. For every NR
simulation, we estimate this property from the $(2,2)$ time-domain waveform as
the mean value of ${{\rm d}f}/{{\rm d}t}$ taken in the range $[t_{\rm
mrg},t_{0}]$. The calibrated relation for ${\Im{(\alpha_{\rm fus})}}$ shows
larger uncertainties compared to $f_{\rm mrg}$, as reported in Tab. 1.
However, the presented relation shows clear trends in the tidal parameter and
in the mass ratio. In particular, large-mass-ratio binaries (i.e. $q\gtrsim
1.5$) show ${\Im{(\alpha_{\rm fus})}}\lesssim 0$ due to tidal disruption.

Another early-PM quantity is the time of the first amplitude minimum $t_{0}$.
This quantity is extracted from the time-domain waveform and can be well
captured by the our relations within ${\sim}10\%$. NR simulations of binaries
with $q\gtrsim 1.5$ generally show increasing $t_{0}$ due to tidal disruption.
Also the calibrated relations for ${\Im{(\alpha_{\rm fus})}}$ and $t_{0}$
include a robust BBH limit for $\kappa^{\rm T}_{2}\to 0$ within the nominal
error bars.

###  IV.2 PM frequencies

We extract the main PM frequency $f_{2}$ from NR PM spectra of the $(2,2)$
mode. Generally, the $f_{2}$ frequency is estimated as the global maximum of
the PM spectrum; however, when modulations are prominent and the PM portions
are short (i.e. $\lesssim 8~{}{\rm ms}$), the $f_{2}$ contribution is no
longer the dominant peak and it needs to be identified in the local maxima. As
shown in Tab. 1, the quasiuniversal relation for $f_{2}$ is accurate to
${\sim}4\%$ at $1{-}\sigma$ level ($6{-}7\%$ at $90\%$ credibility level),
that corresponds to an error of about $100~{}{\rm Hz}$ ($200~{}{\rm Hz}$). The
latter is typically smaller than the FWHM of the spectrum peaks. Figure 3
shows this quasuniversal relation: the frequency $Mf_{2}$ primarily correlates
with the tidal polarizability $\kappa^{\rm T}_{2}$, while mass ratio and spin
contributions mildly affect the overall value of this quantity.

The bottom panel of Fig. 3 shows data points with deviations larger than
$2{-}\sigma$. Around $\kappa^{\rm T}_{2}\simeq 207$, it is possible to
identify a cluster of NR data corresponding to spinning unequal-mass H4
binaries $1.65{+}1.10~{}{\rm M_{\odot}}$ with different combinations of spins
Dietrich  _et al._ (2017a, b). For these large mass-ratio cases, the spin
correction employed in Eq. (19) will be improved in a future work when more
data will be available. The largest residual (${\sim}15\%$) is given by the
non-spinning equal-mass binaries BHB$\Lambda\phi$ $1.50{+}1.50~{}{\rm
M_{\odot}}$ Radice  _et al._ (2017) and BLQ $1.40{+}1.40~{}{\rm M_{\odot}}$
Prakash  _et al._ (2021). In both cases, the remnant collapses into BH shortly
after merger, i.e. $t_{\rm coll}\simeq 3~{}{\rm ms}$, and the determination of
the peak and secondary frequencies from this signal is rather delicate due to
the short duration of the transient. From the Fourier spectra, it is possible
to identify two dominant broad peaks at frequencies $Mf_{2-0}\simeq 0.036$ and
$Mf_{2}\simeq 0.048$ for the BHB$\Lambda\phi$ binary and $Mf_{2-0}\simeq
0.036$ and $Mf_{2}\simeq 0.047$ for the BLQ binary. These values agree with
the estimate of $Mf_{0}$ coming from the instantaneous GW frequency; however,
the peak widths vary depending on the window used to smooth the NR data and it
is not possible to clearly identify a carrier frequency and a modulation
magnitude from the time-domain waveform. Consistently with Ref. Breschi  _et
al._ (2019) we chose to identify the second peak with $f_{2}$ and
conservatively include it in the determination of the quasiuniversal relation.
In contrast, the choice of the first peak as $f_{2}$ would be consistent with
Ref. Radice  _et al._ (2017); Prakash  _et al._ (2021), and the datapoints
would not be outliers in the residual plot.

The value of the frequency $f_{0}$ is estimated as
$f_{0}={T_{0}}^{-1}\simeq(t_{3}-t_{1})^{-1}$ (Sec. III.1). The frequency
$f_{0}$ shows a non-monotonic dependency on the tidal coupling $\kappa^{\rm
T}_{2}$ for $q\simeq 1$. The relative error associated with $f_{0}$ is
${\sim}60\%$, that is considerably larger than the error on the peak frequency
$f_{2}$. This uncertainty can be related to the method used to estimate
$f_{0}$ and to the numerical error that affects amplitude fluctuations. In
principle, the frequency $f_{0}$ can be also extracted from the $(\ell=2,m=0)$
mode of the GW waveform. However, numerical errors appear to be larger for HM
components, due to the lower magnitude of the strains, and the corresponding
spectra do not show neat and unambiguous Fourier peaks, yielding to less
accurate calibrated relations.

###  IV.3 PM amplitudes

The PM amplitudes $A_{1}$, $A_{2}$ and $A_{3}$ are extracted from the time-
domain NR data and they show a decreasing trend for increasing $\kappa^{\rm
T}_{2}$ and for increasing mass ratio, similarly to Ref. Breschi  _et al._
(2019). This can be understood as the effects of stiffer EOSs and larger mass
ratios that produce less violent dynamics in the remnant (for a fixed $M$). As
a consequence of tidal disruption, the first amplitude $A_{0}$ increases with
increasing mass ratio. Overall, these quantities show errors between $15\%$
and $40\%$, except for $A_{0}$, which shows an error $\gtrsim 60\%$ since this
quantity is comparable in magnitude to NR errors.

###  IV.4 Late-time features

The damping time ${\Re{(\beta_{\rm peak})}}$ of the decaying tail in NRPMw is
estimated from NR data using the approximation for exponential sinusoidal
functions, i.e. ${\Re{(\beta_{\rm peak})}}\simeq{\rm max}(A(t))/[2\,{\rm
max}(A(f))]$, where ${\rm max}(A(t))$ is the maximum amplitude of the time-
domain waveform and ${\rm max}(A(f))$ is the maximum amplitude of the
frequency-domain spectrum. Despite errors of ${\sim}30\%$, the calibrated
relation has a physically reasonable trend. For example, ${\Re{(\beta_{\rm
peak})}}$ decreases for increasing mass ratios, in agreement with the tidally
disruptive dynamics of high-mass ratio mergers.

The FM displacement $\Delta_{\rm fm}$ is estimated from the time-domain NR
waveforms as the largest displacement in the instantaneous GW frequency
$f_{22}(t)$ from the PM peak $f_{2}$. The $\Delta_{\rm fm}$ predictions show
similar trends and comparable values to $f_{0}$ for equal-mass binaries. More
significant differences emerge instead as the mass-ratio get larger. The FM
damping time $\Gamma_{\rm fm}$ is also estimated from the time-domain NR data
fitting a damped sinusoidal to the instantaneous GW frequency. This quantity
has the less accurate relation among the presented cases (${\sim}90\%$) due to
the large errors introduced by the extraction method.

##  V Validation

We validate the NRPMw model by computing its faithfulness ${{\mathcal{F}}}$
against 102 NR waveforms of Refs. Radice (2017); Radice  _et al._ (2018a);
Breschi  _et al._ (2019); Nedora  _et al._ (2021); Bernuzzi  _et al._ (2020);
Camilletti  _et al._ (2022) that were not used for the calibration. Among the
considered simulations, 12 binaries show prompt collapse into BH. The
validation set is composed by NR simulations of non-spinning BNS performed
with THC Radice and Rezzolla (2012) that include different neutrino
treatments, turbolent viscosity schemes and five EOSs, i.e. BHB$\Lambda\phi$
Banik  _et al._ (2014), DD2 Typel  _et al._ (2010), LS220 Lattimer and Swesty
(1991), SFHo Steiner  _et al._ (2013) and SLy Douchin and Haensel (2001). The
intrinsic binary properties cover the ranges $M\in[2.6,3.4]~{}{\rm
M_{\odot}}$, $q\in[1,1.8]$ and $\kappa^{\rm T}_{2}\in[47,199]$. The
unfaithfulness ${\bar{\mathcal{F}}}=1-{{\mathcal{F}}}$ between two waveform
templates, say $h_{1}$ and $h_{2}$, is defined as

${\bar{\mathcal{F}}}(h_{1},h_{2})=1-\max_{t_{\rm mrg},\phi_{\rm
mrg}}\frac{(h_{1}|h_{2})}{\sqrt{(h_{1}|h_{1})(h_{2}|h_{2})}}\,,$ (23)

where the maximization is performed over the coalescence time and phase,
respectively $t_{\rm mrg}$ and $\phi_{\rm mrg}$. The inner product
$(h_{1}|h_{2})$ is

$(h_{1}|h_{2})=4\Re\int\frac{h_{1}^{*}(f)\,h_{2}(f)}{S_{n}(f)}\,{\rm d}f\,,$
(24)

where $S_{n}(f)$ is the power spectral density (PSD) of the detector. We
employ the PSD curve of the next-generation detector ET Hild  _et al._ (2011);
Hild (2012) (configuration D). The unfaithfulness is computed between the PM
part of the NR waveform and the NRPMw for the same intrinsic parameters, i.e.
${\bar{\mathcal{F}}}(h_{\rm NR},h_{\text{{\tt NRPMw}}})$, over the frequency
range $[1,8]~{}{\rm kHz}$.

Moreover, we compare the NR faithfulness of NRPMw to that of the time-domain
NRPM model introduced in Ref. Breschi  _et al._ (2019). Note that NRPM can be
also enhanced with the parameters $\\{\alpha,\beta,\phi_{\rm PM}\\}$ Breschi
_et al._ (2022), that are analogous to $\\{\alpha_{\rm peak},t_{\rm
coll},\phi_{\rm PM}\\}$ for ${{\tt NRPMw}}$ and further discussed in App. E.
The main differences between the two models are the following. The frequency
evolution of NRPMw around merger is fully calibrated on NR data, while NRPM
uses a post-Newtonian approximation. The quasiuniversal relations used in NRPM
are not calibrated on SACRA data, although they are compatible with the new
ones computed here for ${{\tt NRPMw}}$. Moreover, NRPMw includes a full
description of damped FM effects and it permits the calibration of the
collapse time $t_{\rm coll}$, that improves the characterization of the
$f_{2}$ peak.

In the following sections, we discuss the introduction of the recalibration
parameters for both time- and frequency-domain models (Sec. V.1) and we
present the unfaithfulness results (Sec. V.2) computed on the independent
validation set of NR data.

###  V.1 Recalibrations

Figure 4: Effect of recalibration terms on NRPMw waveform. The figures show
exemplary templates of the GW plus polarization $h_{+}$ in the time-domain
(left) and in the frequency-domain (right). The template has been computed for
the parameters $M=2.5~{}{\rm M_{\odot}}$, $q=1.08$, $\kappa^{\rm T}_{2}=102$,
$t_{\rm coll}=8~{}{\rm ms}$, $\alpha_{\rm peak}=0.013~{}{\rm kHz}^{2}$,
$\phi_{\rm PM}=\pi/2$ and locating the source at a luminosity distance of
$40~{}{\rm Mpc}$. Black lines show the exact NRPMw predictions, i.e. the
recalibration parameters are identically zero, ${\boldsymbol{\delta}_{\rm
fit}}=0$. The colored lines show three exemplary cases where the values of the
recalibrations ${\boldsymbol{\delta}_{\rm fit}}$ have been randomly extracted
from a zero-mean normal distribution with variance prescribed by the errors of
the residuals.

The EOS-insensitive relations developed in Sec. IV carry intrinsic
uncertainties due to small violations of universality (EOS dependence) and/or
fitting inaccuracies. Calibration errors of the empirical relations should be
taken into account every time such mappings are employed, in particular during
the calculation of fitting factors and during parameter estimation, in order
to perform robust predictions. This can be done by introducing appropriate
parameters associated with the fluctuation of the residuals. A by-product of
this process is that the model can improve its performance in describing the
data.

Labeling $Q$ a generic quantity estimated from a quasiuniversal relation
calibrated on NR data, we introduce an associated recalibration $\delta_{Q}$
that affects the prediction $Q^{\rm fit}$ of the EOS-insensitive relation as

$Q=Q^{\rm fit}\,(1+\delta_{Q})\,.$ (25)

The recalibration $\delta_{Q}$ corresponds to a fractional displacement from
the prediction $Q^{\rm fit}$ of the quasiuniversal relation. The recalibration
procedure employed here is similar to the spectral calibration envelopes used
in GW analyses Vitale  _et al._ (2012). However, here we aim to integrate the
model’s uncertainties in the inference rather than the instrumental errors. A
similar approach has been used in Breschi  _et al._ (2021a) (Sec. 5).

In GW inference applications, the recalibrations of each calibrated PM
property are treated as standard parameters. In this context, it is key that
the prior distribution used in the inference is a good representation of the
residuals of the EOS-insentive relation. This allows us to perform a rigorous
marginalization on the theoretical uncertainties of the model, delivering more
robust and conservative estimates. Interestingly, under the assumption that
the NR error is subdominant compared to the physical breaking of
quasiuniversality, the measurement of the recalibration parameters from the
data could also be used to distinguish between different EOSs and
observatively probe the breaking of quasiuniversality.

A robust characterization of the NR errors is needed in order to employ a
coherent prior distribution for the recalibration parameters in the GW
inference routines. In principle, the uncertainties associated to an EOS-
insensitive relation can be estimated as functions of the employed parameters
using regressive methods or parameter estimation techniques. Following the
methods of Ref. Vitale  _et al._ (2012); Breschi  _et al._ (2021a), an
alternative and simpler approach is to consider a normally distributed prior
distribution with variance prescribed by the errors of the residuals (see Tab.
1). Thus, the relative errors of the EOS-insensitive relations are key
quantities, since they define the theoretical uncertainties of the model.

Figure 4 illustrates the use of recalibrations in NRPMw for an examplary case.
The recalibration parameters ${\boldsymbol{\delta}_{\rm
fit}}=\\{\delta_{i}\\}$ are considered for each element of
${\boldsymbol{\theta}_{\rm fit}}$. These additional degrees of freedom mildly
affect the merger portion, i.e. $t<t_{0}$ due to the accuracy of the empirical
relations close to merger. However, the recalibration coefficients have a
larger effect on the late-time PM features whose EOS-insensitive relations
introduce larger uncertainties. Analogously, the recalibrations can be
introduced for NRPM. This additional flexibility is expected to significantly
improve the data fitting by adjusting the PM morphology of the template to
match the targeted signal, similarly to agnostic approaches (e.g.
Chatziioannou  _et al._ , 2017a; Easter  _et al._ , 2020).

###  V.2 Unfaithfulness

Figure 5: Recovered unfaithfulness ${\bar{\mathcal{F}}}$ between PM models and
NR data of the validation set Radice (2017); Radice  _et al._ (2018a); Breschi
_et al._ (2019); Nedora  _et al._ (2021); Bernuzzi  _et al._ (2020);
Camilletti  _et al._ (2022) employing ET-D sensitivity Hild  _et al._ (2011);
Hild (2012). For NRPM Breschi  _et al._ (2019) (thin lines), we compute
${\bar{\mathcal{F}}}$ with the standard model (a), including PM parameters (b)
and also the recalibrations (c). Analogously, the ${\bar{\mathcal{F}}}$
recovered for NRPMw (thick lines) include the PM parameters (d) and also the
recalibrations (e). The dashed histogram shows the ${\bar{\mathcal{F}}}$ for
case (e) computed over the calibration set.  Figure 6: Comparison between PM
models and exemplary NR data of the validation set. Colored lines show the
spectra for the different models, analogously to Figure 5. Solid lines are
used for NRPMw spectra and dashed lines are employed for NRPM. NR spectra are
reported with black solid lines. The plot includes also the corresponding
unfaithfulnesses estimated with NRPMw model, i.e. case (d) in blue and case
(e) in green.

We compare here the NR faithfulness results for NRPMw and NRPM. In Figure 5,
we report histograms of the unfaithfulness computed on the validation NR
sample of

  * (a)

NRPM without resorting to minimization methods;

  * (b)

NRPM minimizing over the additional PM parameters $\\{\alpha,\beta,\phi_{\rm
PM}\\}$ and setting ${\boldsymbol{\delta}_{\rm fit}}=0$;

  * (c)

NRPM with recalibration parameters ${\boldsymbol{\delta}_{\rm fit}}$ and
minimizing over ${\boldsymbol{\delta}_{\rm fit}}$ and
$\\{\alpha,\beta,\phi_{\rm PM}\\}$;

  * (d)

NRPMw minimizing over the additional PM parameters ${\boldsymbol{\theta}_{\rm
free}}$ and setting ${\boldsymbol{\delta}_{\rm fit}}=0$;;

  * (e)

NRPMw with recalibration parameters ${\boldsymbol{\delta}_{\rm fit}}$ and
minimizing over ${\boldsymbol{\delta}_{\rm fit}}$ and
${\boldsymbol{\theta}_{\rm free}}$.

In particular, the minimization procedure is performed as follows. For each NR
waveform, we compute the corresponding NRPMw (or NRPM) template fixing the
intrinsic parameters $\boldsymbol{\theta}_{\rm bin}$ to the values of the NR
simulation and estimating the additional parameters
(${\boldsymbol{\theta}_{\rm free}}$ and ${\boldsymbol{\delta}_{\rm fit}}$)
minimizing the unfaithfulness ${\bar{\mathcal{F}}}$, i.e. Eq. (23), using a
differential evolution method Storn and Price (1997). For each case and for
each NR data, the additional degrees of freedom are independently varied over
a physically-motivated range 666For NRPMw, we set the time of collapse $t_{\rm
coll}\geq t_{0}$, the frequency drift $M^{2}\alpha_{\rm
peak}\in[-10^{-5},10^{-5}]$, the PM phase $\phi_{\rm PM}\in[0,2\pi]$ and the
recalibrations $\delta_{i}\in[-4\sigma_{i},+4\sigma_{i}]$, where $i$ runs over
the calibrated PM quantities and $\sigma_{i}$ is the corresponding standard
deviation of the NR residuals (see Tab. 1).  in order to estimate the minimum
${\bar{\mathcal{F}}}$.

Case (a) gives results comparable to Breschi  _et al._ (2019), with median
value ${\bar{\mathcal{F}}}$ equal to $0.45$ and few cases with
${\bar{\mathcal{F}}}\leq 0.1$ ($2\%$). Indeed, the only differences between
this work and Ref. Breschi  _et al._ (2019) are the PSD and the different
validation set. In case (b), the inclusion of the free parameters
$\\{\alpha,\beta,\phi_{\rm PM}\\}$ improves the faithfulness of NRPM by
shifting the median value to ${\bar{\mathcal{F}}}\simeq 0.27$, but the
majority of the recovered values ($97\%$) lies above
${\bar{\mathcal{F}}}=0.1$. The additional inclusion of the recalibration
parameters, shown in case (c), considerably enhances the quality of the
recovered waveforms, since ${\bar{\mathcal{F}}}$ decreases with median
${\bar{\mathcal{F}}}={0.06}$ and down to ${\bar{\mathcal{F}}}\sim O(10^{-3})$,
corresponding to short-lived remnants and prompt BH collapses. The fraction of
cases with ${\bar{\mathcal{F}}}<0.1$ corresponds to $83\%$ and we recovered
${\bar{\mathcal{F}}}<0.2$ for all binaries in the validation set.

Moving to the novel NRPMw model, case (d) show an overall improvement in the
faithfulness compared to the equivalent case (b), with median
${\bar{\mathcal{F}}}{\sim}0.13$ and a fraction of $38\%$ with
${\bar{\mathcal{F}}}<0.1$. We attributed this enhancement to the modeling
choices employed in NRPMw, since the number of parameters minimized
(${\boldsymbol{\theta}_{\rm free}}$) is the same as case (b). Moreover, case
(d) shows a small cluster with ${\bar{\mathcal{F}}}\lesssim 3{\times}10^{-2}$
(${\sim}20\%$), mainly populated by short-lived remnant and prompt BH
collapses. In case (e), the additional inclusion of recalibration terms
considerably improves the agreement of NRPMw to the NR data. We obtain a
median ${\bar{\mathcal{F}}}$ of $2.5{\times}10^{-2}$ and report $94\%$ of the
validation set with ${\bar{\mathcal{F}}}<0.1$. We recover similar statistics
applying case (e) over the six-hundred NR simulations of the calibration set,
shown with dashed line in Figure 5. Moreover, the histogram (e) shows that the
cluster constituted by short-duration signals moves toward
${\bar{\mathcal{F}}}=10^{-2}$ and we recovered values comparable to or smaller
than ${\bar{\mathcal{F}}}=3{\times}10^{-2}$ for several long-duration
transients, such as SLy $1.30{+}1.30~{}{\rm M_{\odot}}$, and unequal-mass
binaries, such as DD2 $1.50{+}1.25~{}{\rm M_{\odot}}$. The overall improvement
with respect to the comparable case (c) is roughly half order of magnitude.

The recovered results validate the modeling choices, suggesting that the
primary contributions of the theoretical errors are the inaccurate predictions
of the EOS-insensitive relations. Considering the faithfulness condition
proposed in Ref. Damour  _et al._ (2011); Chatziioannou  _et al._ (2017b);
Gamba  _et al._ (2021b) and fixing $N=9$ as number of intrinsic parameters
$\\{\boldsymbol{\theta}_{\rm bin},{\boldsymbol{\theta}_{\rm free}}\\}$, the
recovered upper-bound accuracy ${\bar{\mathcal{F}}}\simeq 10^{-1}$ of NRPMw in
case (e) can be translated into a model robustness threshold of SNR ${\sim}7$.
Above this threshold, systematic waveform errors can become relevant. The
threshold moves to SNR ${\sim}11$ if we include the recalibrations
${\boldsymbol{\delta}_{\rm fit}}$ as intrinsic parameters, i.e. $N=22$. On the
other hand, employing the recovered median value ${\bar{\mathcal{F}}}\simeq
2.5{\times}10^{-2}$, we estimate a faithfulness threshold SNR equal to $13$
for $N=9$ and $21$ for $N=22$. Considering an averaged threshold of SNR
${\sim}10$, this limit matches the requirements imposed by ET detector for
(optimally-oriented) sources located at luminosity distances $\gtrsim
40~{}{\rm Mpc}$ Tsang  _et al._ (2019); Breschi  _et al._ (2019); Easter  _et
al._ (2020); Breschi  _et al._ (2021b, 2022).

Notably, the ${\bar{\mathcal{F}}}$ values computed on simulations with
different grid resolution or physical schemes suffer from considerable
fluctuations for some binaries. Some examples are: LS220 $1.47{+}1.27~{}{\rm
M_{\odot}}$ that gives $\log_{10}{\bar{\mathcal{F}}}=-0.84$ at standard
resolution without turbolent viscosity and
$\log_{10}{\bar{\mathcal{F}}}=-1.42$ at high resolution with turbolent
viscosity; and LS220 $1.35{+}1.35~{}{\rm M_{\odot}}$ (with turbolent
viscosity) that gives $\log_{10}{\bar{\mathcal{F}}}=-1.05$ at standard
resolution and $\log_{10}{\bar{\mathcal{F}}}=-1.79$ at low resolution. These
results suggest that the largest ${\bar{\mathcal{F}}}$ might be related to an
inaccurate modeling of the late-time features or to an excess of numerical
error in the data. On the other hand, the accuracy of NR templates computed
from different grid resolutions spans the range ${\bar{\mathcal{F}}}\simeq
0.6$ to ${\bar{\mathcal{F}}}\simeq 10^{-2}$ Breschi  _et al._ (2019),
comparably to case (c) and (e). These non-negligible errors originate from
finite resolution of numerical data.

Figure 6 shows the comparison between the PM model spectra and NR data for
four exemplary cases extracted from the validation set. The first case is DD2
$1.509{+}1.235~{}{\rm M_{\odot}}$ which generates a long-lived remnant,
$t_{\rm coll}\sim O(100~{}{\rm ms})$. The NRPM model (a) predicts an erroneous
$f_{2}$ peak, which biases the estimation of the damping time in case (b). The
result improves to $\log_{10}{\bar{\mathcal{F}}}=-1.3$ in case (c). The novel
NRPMw matches well the NR data, delivering $\log_{10}{\bar{\mathcal{F}}}=-1.3$
in case (d) and $\log_{10}{\bar{\mathcal{F}}}=-1.6$ in case (e). The second
case is LS220 $1.635{+}1.146~{}{\rm M_{\odot}}$ with tidal-disruptive behavior
that collapses into BH ${\sim}12~{}{\rm ms}$ after merger. For this
simulation, NRPMw (d) is not capable to match the dominant PM peak returning
${\bar{\mathcal{F}}}\simeq 0.15$. Then, the recalibrations (e) strongly
improve the agreement to NR data, yielding
$\log_{10}{\bar{\mathcal{F}}}=-1.7$. The third case is SFHo
$1.364{+}1.364~{}{\rm M_{\odot}}$ which generates a short-lived remnant with
$t_{\rm coll}\simeq 4~{}{\rm ms}$. This spectrum highlights the relevance of
modulation effects in PM signals. The comparison shows the flexibility of the
recalibrated NRPMw (e) in capturing the several Fourier peaks, delivering
$\log_{10}{\bar{\mathcal{F}}}=-1.9$. The last case is SLy
$1.364{+}1.364~{}{\rm M_{\odot}}$ with $t_{\rm coll}\simeq 12~{}{\rm ms}$ that
shows prominent modulations in the spectrum. NRPM does not match well the
prominent subdominant peaks returning $\log_{10}{\bar{\mathcal{F}}}=-0.7$ in
case (c). This result is similar to NRPMw (d) but considerably improved with
the inclusion of recalibrations (e) to $\log_{10}{\bar{\mathcal{F}}}=-1.8$.

##  VI Conclusions

This paper presents NRPMw, a frequency-domain model for PM GW from BNS
remnants calibrated with EOS-insensitive relations from the largest publicly-
available set of NR simulations. NRPMw is designed to be employed in fully or
partially informed Bayesian inference from GW data. NRPMw includes the
dependency on the intrinsic binary parameters through the EOS-insensitive
relations, thus allowing (i) the direct astrophysical inference of all the BNS
parameters without assuming a pre-merger signal/detection, and at the same
time (ii) a phase-coherent attachment with pre-merger templates Breschi  _et
al._ (2019). The current uncertainties of EOS-insensitive relations can be
taken into account in a partially informed approach using recalibrations
parameters. This enhances the flexibility of the model in capturing the
complex morphology of PM signals and improves fitting factors. We stress that
a recalibration procedure similar to the one introduce here should be employed
every time EOS-insensitive are applied to any type of data.

NRPMw was validated with an independent set of 102 NR simulations. The fitting
factors favorably compares against the results obtained with similar
frequency-domain models Tsang  _et al._ (2019) and significantly improve those
we obtained with NRPM Breschi  _et al._ (2019). The improvement is mainly
related to a more accurate modeling of the merger features and to an improved
description of the FMs when compared to NRPM. The faithfulness of the
recalibrated NRPMw is comparable to that obtained using unmodeled (non-
informed) templates and agnostic approaches Easter  _et al._ (2020); Soultanis
_et al._ (2021). This comes at the cost of 13 recalibration parameters and
three free parameters, compared to the typical $O(10)$ parameters of unmodeled
templates. However, differently from the latter, NRPMw delivers complete
posteriors for the BNS parameters, including mass, mass ratio, etc. The NRPMw
faithfulnesses are comparable to the accuracy of current NR templates for BNS
remnants. The further development of high-precision NR simulations is key for
the design of robust PM models.

NRPMw builds on a new set of EOS-insensitive relations for the PM spectra. We
have focused the development of quasiuniversal relation that employ the tidal
coupling constant $\kappa^{\rm T}_{2}$ in view of utilizing the model as an
EOB completion Bernuzzi  _et al._ (2014b, 2015b); Zappa  _et al._ (2019);
Breschi  _et al._ (2019). The most robust relations we obtained are, not
surprisingly, the merger amplitude, the merger frequency and the dominant
$f_{2}$ peak. The $1{-}\sigma$ uncertainties of these relations are of the
order of 4% due to either uncertainties of NR data or EOS-dependent features.

We found that the presence of softening effects due to quark deconfinement or
hyperonic degrees of freedom at high densities does not introduce significant
deviations above the $2{-}\sigma$ credibility level in the EOS-insensitive
relation for $f_{2}$ developed here. Hence, the observational imprint of EOS
softening might be better revealed from an earlier BH collapse phenomenology,
e.g. Sekiguchi  _et al._ (2011); Radice  _et al._ (2017); Prakash  _et al._
(2021); Huang  _et al._ (2022); Fujimoto  _et al._ (2022), rather than from
the measurement of PM frequencies (under the assumption that our sample of
models adequately represent the “true” EOS). However, the cases
BHB$\Lambda\phi$ $1.50{+}1.50~{}{\rm M_{\odot}}$, BLQ $1.40{+}1.40~{}{\rm
M_{\odot}}$ 777As discussed in Sec. IV.2, these massive binaries of Ref.
Radice  _et al._ (2017); Breschi  _et al._ (2019); Prakash  _et al._ (2021)
are marginal cases very close to prompt collapse and its interpretation is not
fully clear with the present data. and other literature results Bauswein  _et
al._ (2019); Chatziioannou and Han (2020); Most  _et al._ (2020) suggest that
the EOS-insentive relations for $Mf_{2}(\kappa_{2}^{\rm T})$ (and in principle
for other quantities) might break for particular binary masses and in presence
of “strong” phase transitions 888 The term “strong” is often used in the
literature but it does not have any precise meaning; in this context it is
used as a tautology to indicate that the EOS model breaks the quasiuniversal
relation. . This opens the possibility of using NRPMw to identify this new
extreme matter physics via Bayesian analyses following the method of Breschi
_et al._ (2019); Wijngaarden  _et al._ (2022). We stress that these types of
analyses strictly probe only the violation of the particular quasiuniversal
relation that is assumed in a model. Hence, the robust construction of EOS-
insensitive relations and the use of recalibration parameters are key for the
interpretation of the inference results. Future analysis must incorporate the
here proposed recalibration parameters in the inference, because assessing the
breaking of the quasiuniversality requires the knowledge of the theoretical
uncertainty of the EOS-insentive relation.

The use of dimensionless and mass-rescaled quantities is a key aspect in
building EOS-insensitive relations. An example illustrating this fact is the
breaking of the quasiuniversal relations claimed in Ref. Raithel and Most
(2022). The latter refers to relations of type $f_{2}(R)$ that are different
from those employed here. In Appendix F, we verified that those $f_{2}(R)$
relations are broken also by some data of the CoRe database. The additional
term proposed in Raithel and Most (2022) does not fix the breaking of some
CoRe data with softening effects at high densities. However, we verified that
the use of mass-rescaled quantities, i.e. $Mf_{2}(R/M)$, leads to more robust
EOS-insensitive relations. Hence, considering $Mf_{2}(R/M)$ or $f_{2}(R)$ in a
Bayesian analysis of the same data would incorrectly lead to two different
conclusions about the EOS. In a companion paper we report a study on the
application of NRPMw to detection and Bayesian parameter estimation of PM
signal with ET. As anticipated by the faithfulness calculations presented
here, NRPMw can improve the performances of NRPM Breschi  _et al._ (2019),
yielding to threshold SNRs comparable to those of unmodeled analyses
Chatziioannou  _et al._ (2017a); Easter  _et al._ (2020). For example, the
model can be used to infer the dynamical frequency evolution of the remnant
and the time of BH collapse already at the minumum SNR treshold. Under the
important caveat on the robusteness of the assumed EOS-insensitive relations
discussed above, these observables can provide insight into the properties of
matter under extreme conditions. Moreover, NRPMw can be employed together with
inspiral-merger templates to characterize the full GW spectrum of BNS mergers
following Breschi  _et al._ (2019, 2022). Further studies in this directions
will be reported in a third paper of this series.

## Acknowledgments

MB, SB and KC acknowledge support by the EU H2020 under ERC Starting Grant,
no. BinGraSp-714626. MB acknowledges support from the Deutsche
Forschungsgemeinschaft (DFG) under Grant No. 406116891 within the Research
Training Group RTG 2522/1. SB acknowledges the hospitality of KITP at UCSB and
partial support by the National Science Foundation under Grant No. NSF
PHY-1748958 during the conclusion of this work. AC and AP acknowledge PRACE
for awarding them access to Joliot-Curie at GENCI@CEA. They also acknowledge
the usage of computer resources under a CINECA-INFN agreement (allocation
INF20_teongrav and INF21_teongrav). The computations were performed on ARA, a
resource of Friedrich-Schiller-Universtät Jena supported in part by DFG grants
INST 275/334-1 FUGG, INST 275/363-1 FUGG and EU H2020 BinGraSp-714626. and on
the Tullio sever at INFN Turin.

The waveform model developed in this work, NRPMw, is implemented in bajes and
the software is publicly available at:

https://github.com/matteobreschi/bajes.

##  Appendix A Wavelet approximations

We discuss the approximations employed to compute the frequency-domain wavelet
${W}(f)$ Eq. (8) for different values of the $\alpha$ parameter.

When $\alpha$ is identically zero, Eq. (7) reduces to a damped sinusoidal
function, but Eq. (8) leads to an indeterminate form. Then, the latter can be
replaced by

${W}(f)={\rm e}^{\gamma}\,\left(\frac{{\rm e}^{\tau\zeta}-1}{\zeta}\right)\,,$
(26)

where $\zeta(f)=\beta-2\pi{\rm i}f$. Eq. (26) represents a good approximation
also when the wavelet is strongly damped, i.e. $\Re{(\beta)}$ dominates over
the quadratic contributions. This is crucial to simplify the computations for
higher-order FM terms, whose damping time decrease linearly with the
approximation order (see App. (B)).

For $|\alpha|\ll 1$, arithmetic overflows arise in numerical computations.
Then, for these cases, we expand Eq. (7) around small values of $\alpha
t^{2}$, i.e.

${W}(t)={\rm e}^{\beta t+\gamma}\sum_{n=0}^{\infty}\frac{\left(\alpha
t^{2}\right)^{n}}{n!}\,.$ (27)

Each term in Eq. (27) can be analytically integrated, leading to a well-
defined solution of the Fourier counterpart. In particular,

${W}(f)={\rm e}^{\gamma}\,\sum_{n=0}^{\infty}\int_{0}^{\tau}\frac{\left(\alpha
t^{2}\right)^{n}}{n!}\,{\rm e}^{(\beta-2\pi{\rm i}f)t}\,{\rm d}t\,,$ (28)

from which it follows

${W}(f)=\frac{{\rm
e}^{\gamma}}{\sqrt{\pi}}\sum_{n=0}^{\infty}\left(4\alpha\right)^{n}\,\Gamma{\left(n+\frac{1}{2}\right)}\,\frac{G_{2n}(-\zeta\tau)-1}{\zeta^{2n+1}}\,,$
(29)

where $\Gamma(n)$ is the gamma function and $G_{n}(x)$ corresponds to

$G_{n}(x)={\rm e}^{-x}\,\sum_{k=0}^{n}\frac{x^{k}}{k!}\,.$ (30)

We observe that $G_{n}\to 1$ for $n\to\infty$. Limiting the series to $n=0$,
Eq. (29) leads to the damped sinusoidal case, i.e. Eq. (26). In our
implementation, we use Eq. (29) as approximation of ${W}(f)$ for
$|\alpha|\tau^{2}\lesssim 0.1$ accounting up to $n=4$.

##  Appendix B FM approximation

In this appendix, we discuss the approximation performed in order to reach an
analytical form for the FM effects in terms of ${W}(f)$, i.e. Eq. (8).

Let us start considering a generic non-modulated wavelet ${W}(t)$, as the one
in Eq. (7). This term can be decomposed in amplitude and phase, analogously to
Eq. (2), from which we can compute the frequency, that reads

$\omega_{W}(t)=-2\Im{(\alpha)}t-\Im{(\beta)}\,.$ (31)

In order to include damped FMs, we generalize the notion of ${W}$ introducing
${\tilde{W}}$, such that

$\omega_{\tilde{W}}(t)=\omega_{W}(t)-\Delta_{\rm fm}{\rm e}^{-\Gamma_{\rm
fm}t}\sin(\Omega_{\rm fm}t+\phi_{\rm fm})\,,$ (32)

where $\Delta_{\rm fm},\Gamma_{\rm fm},\Omega_{\rm fm},\phi_{\rm
fm}\in\mathbb{R}$ define the modulation, i.e. the frequency displacement
$\Delta_{\rm fm}\geq 0$, the inverse damping time $\Gamma_{\rm fm}$, the
modulation frequency $\Omega_{\rm fm}$ and the initial phase $\phi_{\rm fm}$.
Integrating Eq. (32), the frequency-modulated wavelet ${\tilde{W}}(t)$ can be
rewritten in the time-domain as

${\tilde{W}}(t)={W}(t;\alpha,\beta,\gamma,\tau)\,{\rm e}^{-{\rm
i}{F}(t;\Delta_{\rm fm},\Gamma_{\rm fm},\Omega_{\rm fm},\phi_{\rm fm})}\,,$
(33)

where ${F}(t)$ corresponds to

${F}(t)=\frac{\Delta_{\rm fm}\,{\rm e}^{-\Gamma_{\rm fm}t}}{\Gamma_{\rm
fm}^{2}+\Omega_{\rm fm}^{2}}\,\big{[}\Gamma_{\rm fm}\sin(\Omega_{\rm
fm}t+\phi_{\rm fm})+\Omega_{\rm fm}\cos(\Omega_{\rm fm}t+\phi_{\rm
fm})\big{]}-{F_{0}}\,,$ (34)

with

${F_{0}}=\frac{\Delta_{\rm fm}}{\Gamma_{\rm fm}^{2}+\Omega_{\rm
fm}^{2}}(\Gamma_{\rm fm}\sin\phi_{\rm fm}+\Omega_{\rm fm}\cos\phi_{\rm
fm})\,.$ (35)

Notice that ${F}(t)\in\mathbb{R}$ and ${\rm e}^{-{\rm i}{F}(t)}$ is a unitary
complex factor for every given $t$.

Due to the oscillatory nature of ${F}(t)$, the frequency-domain wavelet
${\tilde{W}}(f)$ cannot be analytically computed using Gaussian integration
rules. Then, we rewrite ${F}(t)$ in terms of exponential functions,

${F}(t)=\frac{{\rm i}\Delta_{\rm fm}}{2|\beta_{\rm fm}|^{2}}\left(\beta_{\rm
fm}^{*}{\rm e}^{-\beta_{\rm fm}t-{\rm i}\phi_{\rm fm}}-\beta_{\rm fm}{\rm
e}^{-\beta_{\rm fm}^{*}t+{\rm i}\phi_{\rm fm}}\right)-{F_{0}}\,,$ (36)

where $\beta_{\rm fm}=\Gamma_{\rm fm}+{\rm i}\Omega_{\rm fm}$. Subsequently,
we expand the exponential ${\rm e}^{-{\rm i}{F}}$, i.e.

${\rm e}^{-{\rm i}{F}(t)}=\sum_{n=0}^{\infty}\frac{\left[-{\rm
i}{F}(t)\right]^{n}}{n!}\,.$ (37)

Combining Eq. (33), (36) and (37), we can write ${\tilde{W}}(t)$ in terms of
${W}(t)$ and perform an analytical Fourier transform. In particular,

${\tilde{W}}(f)={\rm e}^{{\rm
i}{F_{0}}}\,\sum_{n=0}^{\infty}\,\left(\frac{\Delta_{\rm fm}}{2|\beta_{\rm
fm}|^{2}}\right)^{n}\,\frac{w_{n}(f)}{n!}\,,$ (38)

where

$w_{n}(f)=\sum_{k=0}^{n}{n\choose k}\,(\beta_{\rm fm}^{*})^{k}\,(-\beta_{\rm
fm})^{n-k}\,{W}(f;\alpha,\beta_{n,k},\gamma_{n,k},\tau)\,,$ (39)

with

$\begin{split}\beta_{n,k}&=\beta-k\beta_{\rm fm}-(n-k)\beta_{\rm fm}^{*}\\\
&=\beta-n\Gamma_{\rm fm}+{\rm i}(n-2k)\Omega_{\rm fm}\,,\\\
\gamma_{n,k}&=\gamma_{\rm fm}+{\rm i}(n-2k)\phi_{\rm fm}\,.\end{split}$ (40)

and $\\{\alpha,\beta,\gamma,\tau\\}$ are the parameters of the corresponding
non-modulated wavelet.

Eq. (38) generates several Fourier contributions centered around the
frequencies $\Im{(\beta)}\pm n\Omega$, as expected from FM effects. A second
order approximation gives good agreement for small modulation indices, i.e.
$\Delta/\Omega\ll 1$; however, when $\Delta$ is comparable to $\Omega$,
additional terms need to be taken into account for an accurate description.
The maximum order of approximation $n_{\rm max}$ is estimated using an
empirical rule of thumb, $n_{\rm max}\approx
2\,\left(1+{\Delta}/{\Omega}\right)$.

##  Appendix C Choices for wavelet composition

The first contribution, ${W}_{\rm fus}$, corresponds to the fusion of the NS
cores. This term is modeled with a Gaussian wavelet (i.e. ${\rm
i}\beta\in\mathbb{R}$), with initial amplitude, frequency and phase defined by
the values at merger, respectively $A_{\rm mrg}$, $f_{\rm mrg}$ and $\phi_{\rm
mrg}$. The width of the amplitude is fixed as follows,

${\Re{(\alpha_{\rm fus})}}=\frac{\log(A_{0}/A_{\rm mrg})}{t^{2}_{0}}\,,$ (41)

while the frequency slope ${\Im{(\alpha_{\rm fus})}}$ is directly estimated
from NR data. The fusion wavelet ${W}_{\rm fus}$ is truncated at $t_{0}$. It
follows that

$\begin{split}{W}_{\rm fus}(f)={W}(f;\,\,&\alpha={\Re{(\alpha_{\rm
fus})}}-{\rm i}{\Im{(\alpha_{\rm fus})}},\\\ &\beta=-2\pi{\rm i}f_{\rm
mrg},\\\ &\gamma=\log(A_{\rm mrg})-{\rm i}\phi_{\rm mrg},\\\ &\tau=t_{0},\\\
&{\tau_{0}}=0)\,.\end{split}$ (42)

Subsequently, we include an intermediate wavelet ${\tilde{W}}_{\rm bnc}$ that
characterizes the bounce of the remnant after the collision of the NS cores,
corresponding to the time interval $[t_{0},t_{1}]$. The initial amplitude is
determined in order to match the $A_{0}$ and the phase is computed from the
wavelet ${W}_{\rm fus}$ including an additional phase-shift $\phi_{\rm PM}$,
shown by NR simulations Kastaun  _et al._ (2016); Dietrich  _et al._ (2018b),
i.e.

$\phi_{\rm bnc}=\phi_{\rm mrg}+\phi_{\rm PM}+2\pi f_{\rm
mrg}t_{0}+{\Im{(\alpha_{\rm fus})}}t_{0}^{2}\,.$ (43)

The amplitude coefficients, ${\Re{(\alpha_{\rm bnc})}}$ and ${\Re{(\beta_{\rm
bnc})}}$, are chosen such that the amplitude peaks in the first local
amplitude maximum, i.e. $t_{1}$, with value $A_{1}$;

$\displaystyle{\Re{(\alpha_{\rm bnc})}}$
$\displaystyle=\frac{\log(A_{0}/A_{1})}{(t_{1}-t_{0})^{2}}\,,$ (44)
$\displaystyle{\Re{(\beta_{\rm bnc})}}$
$\displaystyle=\frac{2\log(A_{1}/A_{0})}{t_{1}-t_{0}}\,.$ (45)

The frequency is kept constant with value ${\Im{(\beta_{\rm bnc})}}=-2\pi{\rm
i}f_{2}$. Then, including FM effects as discussed in Sec. III.2, we get

$\begin{split}{\tilde{W}}_{\rm bnc}(f)={W}(f;\,\,&\alpha={\Re{(\alpha_{\rm
bnc})}},\\\ &\beta={\Re{(\beta_{\rm bnc})}}-2\pi{\rm i}f_{2},\\\
&\gamma=\log(A_{0})-{\rm i}\phi_{\rm bnc},\\\ &\tau=t_{1}-t_{0},\\\
&{\tau_{0}}=t_{0},\\\ &\Delta_{\rm fm}=\Delta_{\rm fm},\\\ &\Gamma_{\rm
fm}=0,\\\ &\Omega_{\rm fm}=2\pi f_{0},\\\ &\phi_{\rm fm}=\phi_{\rm
fm})\,.\end{split}$ (46)

After $t_{1}$, the remnant is strongly deformed and the quadrupolar radiation
is affected by couplings with subdominant modes, that introduce AMs.
Physically, this phenomenon can be naively interpreted with the presence of
radial pulsation in the mass distribution of the remnant object Bauswein and
Stergioulas (2015). We limit ourselves to the modeling of AMs in the region
$[t_{1},t_{3}]$ taking into account the coupling with the $(2,0)$ mode,
analogously to Ref. Breschi  _et al._ (2019). This pulsating portion of signal
can be approximated using a wavelet ${\widetilde{W}}_{\rm pul}$ of the form,

${\widetilde{W}}_{\rm pul}(t)=A_{1}\left[1-\Delta_{\rm am}\sin^{2}\left(\pi
f_{0}t\right)\right]{\rm e}^{\left[{\Re{(\beta_{\rm pul})}}-2\pi{\rm
i}f_{2}\right]t-{\rm i}\phi_{\rm pul}}\,,$ (47)

where the initial amplitude and phase are chosen to match values of
${\tilde{W}}_{\rm bnc}$ at $t_{1}$, in particular the phase $\phi_{\rm pul}$
corresponds to

$\phi_{\rm pul}=\phi_{\rm bnc}+2\pi f_{2}(t_{1}-t_{0})\,,$ (48)

the coefficient ${\Re{(\beta_{\rm pul})}}$ is defined by the amplitudes
$A_{1,3}$ as

${\Re{(\beta_{\rm
pul})}}=\frac{\log\left({A_{3}}/{A_{1}}\right)}{t_{3}-t_{1}}\,,$ (49)

and the coefficient $\Delta_{\rm am}$ defines the magnitude of AMs,

$\Delta_{\rm
am}=1-\frac{A_{2}}{A_{1}}\left(\frac{A_{1}}{A_{3}}\right)^{\frac{t_{2}-t_{1}}{t_{3}-t_{1}}}=1-\frac{A_{2}}{\sqrt{A_{1}\,A_{3}}}\,,$
(50)

where we made use of the definition of $t_{i}$ (Sec. III.1) in the second
equality. Then, Eq. (47) can be rewritten in terms of frequency-domain
wavelets, Eq. (8), as

$\begin{split}{\widetilde{W}}_{\rm pul}(f)=\left(1-\frac{\Delta_{\rm
am}}{2}\right){W}(f;\,\,&\alpha=0,\\\ &\beta={\Re{(\beta_{\rm pul})}}-2\pi{\rm
i}f_{2},\\\ &\gamma=\log(A_{1})-{\rm i}\phi_{\rm pul},\\\
&\tau=t_{3}-t_{1},\\\ &{\tau_{0}}=t_{1},\\\ &\Delta_{\rm fm}=\Delta_{\rm
fm},\\\ &\Gamma_{\rm fm}=\Gamma_{\rm fm},\\\ &\Omega_{\rm fm}=2\pi f_{0},\\\
&\phi_{\rm fm}=\phi_{\rm fm})\\\ +\,\,\frac{\Delta_{\rm
am}}{4}\,\,{W}(f;\,\,&\alpha=0,\\\ &\beta={\Re{(\beta_{\rm pul})}}-2\pi{\rm
i}(f_{2}-f_{0}),\\\ &\gamma=\log(A_{1})-{\rm i}\phi_{\rm pul},\\\
&\tau=t_{3}-t_{1},\\\ &{\tau_{0}}=t_{1},\\\ &\Delta_{\rm fm}=\Delta_{\rm
fm},\\\ &\Gamma_{\rm fm}=\Gamma_{\rm fm},\\\ &\Omega_{\rm fm}=2\pi f_{0},\\\
&\phi_{\rm fm}=\phi_{\rm fm})\\\ +\,\,\frac{\Delta_{\rm
am}}{4}\,\,{W}(f;\,\,&\alpha=0,\\\ &\beta={\Re{(\beta_{\rm pul})}}-2\pi{\rm
i}(f_{2}+f_{0}),\\\ &\gamma=\log(A_{1})-{\rm i}\phi_{\rm pul},\\\
&\tau=t_{3}-t_{1},\\\ &{\tau_{0}}=t_{1},\\\ &\Delta_{\rm fm}=\Delta_{\rm
fm},\\\ &\Gamma_{\rm fm}=\Gamma_{\rm fm},\\\ &\Omega_{\rm fm}=2\pi f_{0},\\\
&\phi_{\rm fm}=\phi_{\rm fm})\,.\end{split}$ (51)

Subsequently, we model the signal with a damped tail related to the
quadrupolar deformations of the rotating remnant. The corresponding wavelet
${\tilde{W}}_{\rm peak}$ is modeled in the range $[t_{3},t_{\rm coll}]$. If
the remnant is a stable NS configuration, $t_{\rm coll}\to\infty$. The initial
amplitude and phase are chosen to match the values of ${\widetilde{W}}_{\rm
pul}$ at $t_{3}$,

$\displaystyle\phi_{\rm peak}$ $\displaystyle=\phi_{\rm pul}+2\pi
f_{2}(t_{3}-t_{1})\,,$ (52) $\displaystyle A_{\rm peak}$
$\displaystyle=A_{3}\,.$ (53)

The frequency evolution is characterized by the typical $f_{2}$ peak, i.e.
${\Im{(\beta_{\rm peak})}}=-2\pi f_{2}$, with a non-vanishing slope
${\Im{(\alpha_{\rm peak})}}$ (also referred as $\alpha_{\rm peak}$ in the
manuscript to lighten the notation). Then,

$\begin{split}{\tilde{W}}_{\rm peak}(f)={W}(f;\,\,&\alpha=-{\rm
i}\Im{(\alpha_{\rm peak})},\\\ &\beta={\Re{(\beta_{\rm peak})}}-2\pi{\rm
i}f_{2},\\\ &\gamma=\log(A_{3})-{\rm i}\phi_{\rm peak},\\\ &\tau=t_{\rm
coll}-t_{3},\\\ &{\tau_{0}}=t_{3},\\\ &\Delta_{\rm fm}=\Delta^{\prime}_{\rm
fm},\\\ &\Gamma_{\rm fm}=\Gamma_{\rm fm},\\\ &\Omega_{\rm fm}=2\pi f_{0},\\\
&\phi_{\rm fm}=\phi_{\rm fm})\,,\end{split}$ (54)

where $\Delta^{\prime}_{\rm fm}=\Delta_{\rm fm}\,\exp[\Gamma_{\rm
fm}(t_{3}-t_{1})]$.

When $t_{\rm coll}$ is finite and the remnant collapse into BH, NR simulations
show an increasing frequency and a damping amplitude similarly to a BH
ringdown. This evolution can be captured with the inclusion of an additional
wavelet component, i.e. ${W}_{\rm coll}$. However, this contribution is
expected to be relatively weak in terms of GW luminosity with respect to the
previous dynamics Zappa  _et al._ (2018). Moreover, the characteristic BH
frequencies for this kind of systems lie in a very high frequency range,
roughly $\gtrsim 6$ kHz Berti and Klein (2014), where the sensitivities of the
detectors are generally poor. It follows that the collapse portion of the
signal is expected to have negligible effect on the overall GW power and, for
these reasons, we approximate ${W}_{\rm coll}=0$.

The overall model includes 17 parameters: the merger amplitude $A_{\rm mrg}$
and frequency $f_{\rm mrg}$; the frequency drift at merger ${\Im{(\alpha_{\rm
fus})}}$; the characteristic PM frequencies $f_{2}$ and $f_{0}$; the frequency
drift ${\Im{(\alpha_{\rm peak})}}$ (or $\alpha_{\rm peak}$); the time of the
first nodal point $t_{0}$ and the corresponding phase-shift $\phi_{\rm PM}$;
the amplitude values at the different nodal points $\\{A_{i}\\}$, for
$i=0,1,2,3$; the inverse damping time of the Lorentzian tail ${\Im{(\beta_{\rm
peak})}}$; the time of collapse $t_{\rm coll}$; and the FM properties, i.e.
$\Delta_{\rm fm}$ and $\Gamma_{\rm fm}$ and$\phi_{\rm fm}$. Finally, we
observe that ${\tilde{W}}_{\rm bnc}(f)$, ${\widetilde{W}}_{\rm pul}(f)$ and
${\tilde{W}}_{\rm peak}(f)$ are chosen to be identically zero for
$\Lambda_{1}=\Lambda_{2}=0$.

##  Appendix D Viscosity impact on frequency drift

Figure 7: GW data extracted from NR simulations with neutrino reabsorption of
a BNS merger with $M=2.7~{}{\rm M_{\odot}}$, $q=1$ and LS220 EOS Radice
(2017); Radice  _et al._ (2018a). Blue curve refers to data with no turbolent
viscosity, orange curve refer to fixed $\ell_{\rm mix}=5~{}{\rm m}$, green
curve refers to fixed $\ell_{\rm mix}=25~{}{\rm m}$ and purple curve refers to
$\ell_{\rm mix}$ calibrated on Kiuchi  _et al._ (2018). Solid lines show the
GW frequency $Mf$ and shaded lines show the GW waveform $h/M$. The instant
$t=0$ corresponds to the merger.

In this Appendix, we show the effects of different viscosity schemes on the
dynamical evolution of the GW frequency. This discussion aims to motivate the
introduction of the frequency drift $\alpha_{\rm peak}$ as free parameter.

Figure 7 shows three NR simulations extracted from Radice (2017); Radice  _et
al._ (2018a) and computed with identical grid resolution. The data correspond
to a BNS system with $M=2.7~{}{\rm M_{\odot}}$ and $q=1$ with matter
properties described by the same EOS, i.e. LS220 Lattimer and Swesty (1991).
Moreover, all cases include neutrino reabsorption scheme Radice  _et al._
(2016b). The blue curves refer to binaries with no turbolent viscosity. The
orange and green curves include turbolent viscosity with a fixed mixing length
respectively equal to $\ell_{\rm mix}=5~{}{\rm m}$ and $\ell_{\rm
mix}=25~{}{\rm m}$. The mixing length $\ell_{\rm mix}$ represents the
characteristic scale over which turbulence acts Radice (2017). Finally, the
purple data are simulated with a turbolent viscosity scheme calibrated on
high-resolution magneto-hydro-dynamical simulations of BNS mergers Kiuchi  _et
al._ (2018).

Over the domain $t/M\lesssim 300$, the different cases show a similar
behavior. However, for later times, the frequency drift significantly differs,
with a more pronounced slope for the $\ell_{\rm mix}=25~{}{\rm m}$ case. The
same binary is the one that shows the earliest BH collapse. Notably, as shown
in Radice (2017), the frequency slope is softer for $\ell_{\rm mix}=50~{}{\rm
m}$ with respect to the $\ell_{\rm mix}=25~{}{\rm m}$ case; however, the
assumption $\ell_{\rm mix}=50~{}{\rm m}$ appears to be physically disfavored
from studies of magnetorotational instability turbolence in BNS simulations
Kiuchi  _et al._ (2018). On the other hand, the simulation with calibrated
$\ell_{\rm mix}$ shows an initial trend similar to the $\ell_{\rm
mix}=5~{}{\rm m}$ case; later, for $t/m\gtrsim 1000$, the two GW frequencies
differ and the calibrated-viscosity case shows a BH collapse. Interestingly,
the binaries with the steepest frequency drifts tend to generate shorter GW
bursts due to earlier BH collapse. These physical effects cannot be described
by the binary properties only (i.e. masses, spins and tides). Thus, it is
necessary to rely on additional coefficients that aim to characterize the
physical information on the matter dynamics encoded in the PM transients.

##  Appendix E Additional parameters for NRPM

In Breschi  _et al._ (2022) and in Sec. V.2, we introduced the additional
parameters $\\{\alpha,\beta,\phi_{\rm PM}\\}$ for the NRPM model Breschi  _et
al._ (2019). This Appendix aims to expand this discussion, specifying the role
of each term with reference to Breschi  _et al._ (2019).

The additional phase $\phi_{\rm PM}$ affects the phase evolution of NRPM
introducing a phase discontinuity in $t_{0}$, as previously discussed for
NRPMw. The damping time $\alpha$, defined in Eq. (11) of Breschi  _et al._
(2019), is promoted to additional parameter, since it improves the fitting of
the characteristic peak for NRPM. This is also motivated by the relation
between the time of collapse and the high-density EOS properties Radice  _et
al._ (2017); Prakash  _et al._ (2021); Fujimoto  _et al._ (2022). Finally, the
$\beta$ parameter is inspired by the template model introduced in Easter  _et
al._ (2020) and it takes into account linear deviations from the peak
frequency $f_{2}$. In particular, it modifies Eq. (7d) of Breschi  _et al._
(2019) as

$\hat{\omega}(\hat{t}>\hat{t}_{3})=\hat{\omega}_{2}\cdot\left[1+\beta(\hat{t}-\hat{t}_{3})\right]\,.$
(55)

As discussed for NRPMw and in App. D, general PM GW transients show non-
vanishing frequency slope and this feature appears to correlate with the
viscosity scheme employed in the NR simulation.

##  Appendix F Quasiuniversal relations of type $f_{2}(R)$

Figure 8: Predicted values from the calibrated relations Eq. (56) compared to
the respective NR observed quantities. Top panels show $X=1.4$ and bottom
panels show $X=1.8$. Left panels show non-mass-scaled $f_{2}$ and right panels
show mass-scaled dimensionless $Mf_{2}$. The diagonal (black line) represents
the case in which predictions and observations match and the gray area is the
90% credibility level. The CoRe data are reported with circles colored
according to $R_{1.4}/R_{1.8}$ and magenta crosses are the data extracted from
Most and Raithel (2021).

Figure 9: Predicted values from the calibrated relations Eq. (57) compared to
the respective NR observed quantities. Top panels show $X=1.4$ and bottom
panels show $X=1.8$. Left panels show non-mass-scaled $f_{2}$ and right panels
show mass-scaled dimensionless $Mf_{2}$. The diagonal (black line) represents
the case in which predictions and observations match and the gray area is the
90% credibility level. The CoRe data are reported with circles colored
according to $R_{1.4}/R_{1.8}$ and magenta crosses are the data extracted from
Most and Raithel (2021). Table 2: Summary of the calibrated relations for the
PM peak frequency $f_{2}$ as function of the NS radii $R_{1.4}$ and $R_{1.8}$.
The first column shows the calibrated quantity of interest; the calibrated
values of the empirical coefficients are reported from the second to the fifth
column.

$Q^{\rm fit}$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$
---|---|---|---|---
$f_{2}(R_{1.4})$ | $5.42$ | $0.0449$ | $-0.0198$ | –
$f_{2}(R_{1.8})$ | $11.5$ | $-0.990$ | $0.0233$ | –
$Mf_{2}\left({R_{1.4}}/{M}\right)$ | $0.2$ | $-0.0762$ | $0.0078$ | –
$Mf_{2}\left({R_{1.8}}/{M}\right)$ | $0.236$ | $-0.103$ | $0.0125$ | –
$f_{2}\left(R_{1.4},{R_{1.4}}/{R_{1.8}}\right)$ | $6.4$ | $-1.33$ | $0.0381$ | $6.97$
$f_{2}\left(R_{1.8},{R_{1.4}}/{R_{1.8}}\right)$ | $9.99$ | $-1.24$ | $0.0349$ | $2.76$
$Mf_{2}\left({R_{1.4}}/{M},{R_{1.4}}/{R_{1.8}}\right)$ | $0.162$ | $-0.115$ | $0.0145$ | $0.0919$
$Mf_{2}\left({R_{1.8}}/{M},{R_{1.4}}/{R_{1.8}}\right)$ | $0.213$ | $-0.105$ | $0.013$ | $0.0241$

In this appendix, we discuss the quasiuniversal relations between the PM peak
frequency $f_{2}$ and the NS radius at fiducial values in light of the results
of Most and Raithel (2021); Raithel and Most (2022). In particular, we
calibrate the relations $f_{2}(R_{1.4})$ and $f_{2}(R_{1.8})$ including the
CoRe data Dietrich  _et al._ (2018b); Radice  _et al._ (2018a); Bernuzzi  _et
al._ (2020); Nedora  _et al._ (2021); Prakash  _et al._ (2021), where
$R_{1.4}$ ($R_{1.8}$) is the radius of a $1.4~{}{\rm M_{\odot}}$ ($1.8~{}{\rm
M_{\odot}}$) NS computed from the TOV equations.

Following Raithel and Most (2022), we employ a quadratic relation for the
calibration of the PM peak as function of the NS radius, including linear
corrections in the ratio ${R_{1.4}}/{R_{1.8}}$, i.e.

$\displaystyle f_{2}(R_{X})$
$\displaystyle=a_{0}+a_{1}R_{X}+a_{2}R_{X}^{2}\,,$ (56) $\displaystyle
f_{2}\left(R_{X},\frac{R_{1.4}}{R_{1.8}}\right)$
$\displaystyle=a_{0}+a_{1}R_{X}+a_{2}R_{X}^{2}+a_{3}\frac{R_{1.4}}{R_{1.8}}\,,$
(57)

for $X=1.4,1.8$, where $f_{2}$ is measured in kHz and $R_{X}$ in km.
Subsequently, we fit the NR data scaling the calibrated quantities by the
total mass $M$ of the system, i.e. $f_{2}\mapsto Mf_{2}$ and $R_{X}\mapsto
R_{X}/M$. Our final calibration set is composed by $65\%$ by binaries with
$R_{1.4}/R_{1.8}>1$. Table 2 shows the values of the calibrated coefficients
$\\{a_{i}\\}$ for the different quantities. Figure 8 and Figure 9 show the NR
data $f_{2}^{\rm NR}$ plotted against the predictions $f_{2}^{\rm fit}$ of the
calibrated relation and the statistical quantities of interest; i.e. the
$\chi^{2}$ (defined in Sec. IV), the adjusted coefficient of determination
$\mathcal{R}^{2}$, the Bayesian information criterion (BIC) and the Akaike
information criterion (AIC). For the computation of the BIC and the AIC, we
define a log-likelihood from Eq. (22) equal to $-\frac{1}{2}\chi^{2}$.

From our analysis, the calibrations performed with the mass-scaled quantities
show improved trends with respect to the analogous non-mass-scaled case. This
is due to the factorization of the total binary mass $M$, as expected by basic
arguments in general relativity. Moreover, the additional contribution
$R_{1.4}/R_{1.8}$ appears to be more relevant for the calibration of low-
density properties, i.e. $f_{2}(R_{1.4})$, in agreement with Raithel and Most
(2022). However, the BIC and the AIC do not favor the introduction of these
additional term, even if the $\chi^{2}$ of the calibrated relation
$Mf_{2}(R_{1.4}/M)$ slightly decreases including $R_{1.4}/R_{1.8}$ in the fit.
We find that the most robust and reliable quasiuniversal relation is the
mapping $Mf_{2}(R_{1.8}/M)$, which reinforces the hypothesis that PM
quantities correlate with high-density EOS properties Breschi  _et al._
(2022). The differences between our results and the findings of Raithel and
Most (2022) might be related to the different size and composition of the NR
set, such as a different set of EOSs, or to different definitions in the
statistical quantities 999We verified that our results are stable when
employing a Gaussian likelihood or the standard Pearson’s $\chi^{2}$
statistic.. Notably, simulations of the CoRe database show deviations
comparable to the cases presented in Most and Raithel (2021); Raithel and Most
(2022) that cannot be fully cured with the introduction of the additional term
$R_{1.4}/R_{1.8}$.

## References

  * Hild  _et al._ (2008) S. Hild, S. Chelkowski, and A. Freise, (2008), arXiv:0810.0604 [gr-qc] .
  * Hild  _et al._ (2011) S. Hild _et al._ , Class. Quant. Grav. 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Hild (2012) S. Hild, Class.Quant.Grav. 29, 124006 (2012), arXiv:1111.6277 [gr-qc] .
  * Punturo  _et al._ (2010a) M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, _et al._ , Class.Quant.Grav. 27, 194002 (2010a).
  * Maggiore  _et al._ (2020) M. Maggiore _et al._ , JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Bernuzzi  _et al._ (2015a) S. Bernuzzi, A. Nagar, T. Dietrich, and T. Damour, Phys.Rev.Lett. 114, 161103 (2015a), arXiv:1412.4553 [gr-qc] .
  * Nagar  _et al._ (2018) A. Nagar _et al._ , Phys. Rev. D98, 104052 (2018), arXiv:1806.01772 [gr-qc] .
  * Akcay  _et al._ (2021) S. Akcay, R. Gamba, and S. Bernuzzi, Phys. Rev. D 103, 024014 (2021), arXiv:2005.05338 [gr-qc] .
  * Gamba  _et al._ (2021a) R. Gamba, S. Bernuzzi, and A. Nagar, Phys. Rev. D 104, 084058 (2021a), arXiv:2012.00027 [gr-qc] .
  * Bernuzzi  _et al._ (2015b) S. Bernuzzi, T. Dietrich, and A. Nagar, Phys. Rev. Lett. 115, 091101 (2015b), arXiv:1504.01764 [gr-qc] .
  * Breschi  _et al._ (2019) M. Breschi, S. Bernuzzi, F. Zappa, M. Agathos, A. Perego, D. Radice, and A. Nagar, Phys. Rev. D100, 104029 (2019), arXiv:1908.11418 [gr-qc] .
  * Damour and Nagar (2009) T. Damour and A. Nagar, Phys. Rev. D80, 084035 (2009), arXiv:0906.0096 [gr-qc] .
  * Damour  _et al._ (2012) T. Damour, A. Nagar, and L. Villain, Phys.Rev. D85, 123007 (2012), arXiv:1203.4352 [gr-qc] .
  * Agathos  _et al._ (2020) M. Agathos, F. Zappa, S. Bernuzzi, A. Perego, M. Breschi, and D. Radice, Phys. Rev. D101, 044006 (2020), arXiv:1908.05442 [gr-qc] .
  * Kashyap  _et al._ (2021) R. Kashyap _et al._ , (2021), arXiv:2111.05183 [astro-ph.HE] .
  * Prakash  _et al._ (2021) A. Prakash, D. Radice, D. Logoteta, A. Perego, V. Nedora, I. Bombaci, R. Kashyap, S. Bernuzzi, and A. Endrizzi, Phys. Rev. D 104, 083029 (2021), arXiv:2106.07885 [astro-ph.HE] .
  * Breschi  _et al._ (2022) M. Breschi, S. Bernuzzi, D. Godzieba, A. Perego, and D. Radice, Phys. Rev. Lett. 128, 161102 (2022), arXiv:2110.06957 [gr-qc] .
  * Breschi  _et al._ (2021a) M. Breschi, A. Perego, S. Bernuzzi, W. Del Pozzo, V. Nedora, D. Radice, and D. Vescovi, Mon. Not. Roy. Astron. Soc. 505, 1661 (2021a), arXiv:2101.01201 [astro-ph.HE] .
  * Breschi  _et al._ (tion) M. Breschi _et al._ , (in preparation).
  * Shibata and Taniguchi (2006) M. Shibata and K. Taniguchi, Phys.Rev. D73, 064027 (2006), arXiv:astro-ph/0603145 [astro-ph] .
  * Hotokezaka  _et al._ (2013) K. Hotokezaka, K. Kiuchi, K. Kyutoku, T. Muranushi, Y.-i. Sekiguchi, _et al._ , Phys.Rev. D88, 044026 (2013), arXiv:1307.5888 [astro-ph.HE] .
  * Takami  _et al._ (2015) K. Takami, L. Rezzolla, and L. Baiotti, Phys.Rev. D91, 064001 (2015), arXiv:1412.3240 [gr-qc] .
  * Bauswein and Stergioulas (2015) A. Bauswein and N. Stergioulas, Phys. Rev. D91, 124056 (2015), arXiv:1502.03176 [astro-ph.SR] .
  * Zappa  _et al._ (2018) F. Zappa, S. Bernuzzi, D. Radice, A. Perego, and T. Dietrich, Phys. Rev. Lett. 120, 111101 (2018), arXiv:1712.04267 [gr-qc] .
  * Radice  _et al._ (2017) D. Radice, S. Bernuzzi, W. Del Pozzo, L. F. Roberts, and C. D. Ott, Astrophys. J. 842, L10 (2017), arXiv:1612.06429 [astro-ph.HE] .
  * Bauswein  _et al._ (2019) A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Phys. Rev. Lett. 122, 061102 (2019), arXiv:1809.01116 [astro-ph.HE] .
  * Most and Raithel (2021) E. R. Most and C. A. Raithel, (2021), arXiv:2107.06804 [astro-ph.HE] .
  * Abbott  _et al._ (2018a) B. P. Abbott _et al._ (VIRGO, KAGRA, LIGO Scientific), Living Rev. Rel. 21, 3 (2018a), [Living Rev. Rel.19,1(2016)], arXiv:1304.0670 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Astrophys. J. 851, L16 (2017b), arXiv:1710.09320 [astro-ph.HE] .
  * Abbott  _et al._ (2018b) B. P. Abbott _et al._ (LIGO Scientific, Virgo),  (2018b), 10.3847/1538-4357/ab0f3d, arXiv:1810.02581 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Punturo  _et al._ (2010b) M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, _et al._ , Class.Quant.Grav. 27, 084007 (2010b).
  * Torres-Rivas  _et al._ (2019) A. Torres-Rivas, K. Chatziioannou, A. Bauswein, and J. A. Clark, Phys. Rev. D99, 044014 (2019), arXiv:1811.08931 [gr-qc] .
  * Martynov  _et al._ (2019) D. Martynov _et al._ , Phys. Rev. D99, 102004 (2019), arXiv:1901.03885 [astro-ph.IM] .
  * Ackley  _et al._ (2020) K. Ackley _et al._ , Publ. Astron. Soc. Austral. 37, e047 (2020), arXiv:2007.03128 [astro-ph.HE] .
  * Clark  _et al._ (2016) J. A. Clark, A. Bauswein, N. Stergioulas, and D. Shoemaker, Class. Quant. Grav. 33, 085003 (2016), arXiv:1509.08522 [astro-ph.HE] .
  * Chatziioannou  _et al._ (2017a) K. Chatziioannou, J. A. Clark, A. Bauswein, M. Millhouse, T. B. Littenberg, and N. Cornish, Phys. Rev. D96, 124035 (2017a), arXiv:1711.00040 [gr-qc] .
  * Easter  _et al._ (2019) P. J. Easter, P. D. Lasky, A. R. Casey, L. Rezzolla, and K. Takami, Phys. Rev. D 100, 043005 (2019), arXiv:1811.11183 [gr-qc] .
  * Bose  _et al._ (2018) S. Bose, K. Chakravarti, L. Rezzolla, B. S. Sathyaprakash, and K. Takami, Phys. Rev. Lett. 120, 031102 (2018), arXiv:1705.10850 [gr-qc] .
  * Tsang  _et al._ (2019) K. W. Tsang, T. Dietrich, and C. Van Den Broeck, Phys. Rev. D100, 044047 (2019), arXiv:1907.02424 [gr-qc] .
  * Easter  _et al._ (2020) P. J. Easter, S. Ghonge, P. D. Lasky, A. R. Casey, J. A. Clark, F. H. Vivanco, and K. Chatziioannou, Phys. Rev. D 102, 043011 (2020), arXiv:2006.04396 [astro-ph.HE] .
  * Soultanis  _et al._ (2021) T. Soultanis, A. Bauswein, and N. Stergioulas, (2021), arXiv:2111.08353 [astro-ph.HE] .
  * Wijngaarden  _et al._ (2022) M. Wijngaarden, K. Chatziioannou, A. Bauswein, J. A. Clark, and N. J. Cornish, (2022), arXiv:2202.09382 [gr-qc] .
  * Whittaker  _et al._ (2022) T. Whittaker, W. E. East, S. R. Green, L. Lehner, and H. Yang, (2022), arXiv:2201.06461 [gr-qc] .
  * Damour and Nagar (2010) T. Damour and A. Nagar, Phys. Rev. D81, 084016 (2010), arXiv:0911.5041 [gr-qc] .
  * Favata (2014) M. Favata, Phys.Rev.Lett. 112, 101101 (2014), arXiv:1310.8288 [gr-qc] .
  * Dietrich  _et al._ (2017a) T. Dietrich, M. Ujevic, W. Tichy, S. Bernuzzi, and B. Brügmann, Phys. Rev. D95, 024029 (2017a), arXiv:1607.06636 [gr-qc] .
  * Shibata and Uryu (2002) M. Shibata and K. Uryu, Prog. Theor. Phys. 107, 265 (2002), arXiv:gr-qc/0203037 .
  * Stergioulas  _et al._ (2011) N. Stergioulas, A. Bauswein, K. Zagkouris, and H.-T. Janka, Mon.Not.Roy.Astron.Soc. 418, 427 (2011), arXiv:1105.0368 [gr-qc] .
  * Bauswein and Janka (2012) A. Bauswein and H.-T. Janka, Phys.Rev.Lett. 108, 011101 (2012), arXiv:1106.1616 [astro-ph.SR] .
  * Bauswein  _et al._ (2012) A. Bauswein, H. Janka, K. Hebeler, and A. Schwenk, Phys.Rev. D86, 063001 (2012), arXiv:1204.1888 [astro-ph.SR] .
  * Takami  _et al._ (2014) K. Takami, L. Rezzolla, and L. Baiotti, Phys.Rev.Lett. 113, 091104 (2014), arXiv:1403.5672 [gr-qc] .
  * Radice  _et al._ (2016a) D. Radice, S. Bernuzzi, and C. D. Ott, Phys. Rev. D94, 064011 (2016a), arXiv:1603.05726 [gr-qc] .
  * Lehner  _et al._ (2016) L. Lehner, S. L. Liebling, C. Palenzuela, O. L. Caballero, E. O’Connor, M. Anderson, and D. Neilsen, Class. Quant. Grav. 33, 184002 (2016), arXiv:1603.00501 [gr-qc] .
  * Dietrich  _et al._ (2017b) T. Dietrich, S. Bernuzzi, M. Ujevic, and W. Tichy, Phys. Rev. D95, 044045 (2017b), arXiv:1611.07367 [gr-qc] .
  * Bernuzzi  _et al._ (2020) S. Bernuzzi _et al._ , Mon. Not. Roy. Astron. Soc. (2020), 10.1093/mnras/staa1860, arXiv:2003.06015 [astro-ph.HE] .
  * Perego  _et al._ (2021) A. Perego, D. Logoteta, D. Radice, S. Bernuzzi, R. Kashyap, A. Das, S. Padamata, and A. Prakash, (2021), arXiv:2112.05864 [astro-ph.HE] .
  * Kastaun  _et al._ (2017) W. Kastaun, R. Ciolfi, A. Endrizzi, and B. Giacomazzo, Phys. Rev. D96, 043019 (2017), arXiv:1612.03671 [astro-ph.HE] .
  * Radice  _et al._ (2020) D. Radice, S. Bernuzzi, and A. Perego, Ann. Rev. Nucl. Part. Sci. 70 (2020), 10.1146/annurev-nucl-013120-114541, arXiv:2002.03863 [astro-ph.HE] .
  * Bernuzzi (2020) S. Bernuzzi, Gen. Rel. Grav. 52, 108 (2020), arXiv:2004.06419 [astro-ph.HE] .
  * Radice  _et al._ (2018a) D. Radice, A. Perego, K. Hotokezaka, S. A. Fromm, S. Bernuzzi, and L. F. Roberts, Astrophys. J. 869, 130 (2018a), arXiv:1809.11161 [astro-ph.HE] .
  * Hotokezaka  _et al._ (2011) K. Hotokezaka, K. Kyutoku, H. Okawa, M. Shibata, and K. Kiuchi, Phys.Rev. D83, 124008 (2011), arXiv:1105.4370 [astro-ph.HE] .
  * Bauswein  _et al._ (2013) A. Bauswein, T. Baumgarte, and H. T. Janka, Phys.Rev.Lett. 111, 131101 (2013), arXiv:1307.5191 [astro-ph.SR] .
  * Bauswein  _et al._ (2020) A. Bauswein, S. Blacker, V. Vijayan, N. Stergioulas, K. Chatziioannou, J. A. Clark, N.-U. F. Bastian, D. B. Blaschke, M. Cierniak, and T. Fischer, Phys. Rev. Lett. 125, 141103 (2020), arXiv:2004.00846 [astro-ph.HE] .
  * Bernuzzi  _et al._ (2016) S. Bernuzzi, D. Radice, C. D. Ott, L. F. Roberts, P. Moesta, and F. Galeazzi, Phys. Rev. D94, 024023 (2016), arXiv:1512.06397 [gr-qc] .
  * Radice  _et al._ (2018b) D. Radice, A. Perego, S. Bernuzzi, and B. Zhang, Mon. Not. Roy. Astron. Soc. 481, 3670 (2018b), arXiv:1803.10865 [astro-ph.HE] .
  * Shibata and Uryu (2000) M. Shibata and K. Uryu, Phys. Rev. D61, 064001 (2000), arXiv:gr-qc/9911058 .
  * Bernuzzi  _et al._ (2014a) S. Bernuzzi, T. Dietrich, W. Tichy, and B. Brügmann, Phys.Rev. D89, 104021 (2014a), arXiv:1311.4443 [gr-qc] .
  * Dimmelmeier  _et al._ (2006) H. Dimmelmeier, N. Stergioulas, and J. A. Font, Mon. Not. Roy. Astron. Soc. 368, 1609 (2006), arXiv:astro-ph/0511394 .
  * Passamonti  _et al._ (2007) A. Passamonti, N. Stergioulas, and A. Nagar, Phys. Rev. D75, 084038 (2007), gr-qc/0702099 .
  * Baiotti  _et al._ (2009) L. Baiotti, S. Bernuzzi, G. Corvino, R. De Pietri, and A. Nagar, Phys. Rev. D79, 024002 (2009), arXiv:0808.4002 [gr-qc] .
  * Nedora  _et al._ (2021) V. Nedora, S. Bernuzzi, D. Radice, B. Daszuta, A. Endrizzi, A. Perego, A. Prakash, M. Safarzadeh, F. Schianchi, and D. Logoteta, Astrophys. J. 906, 98 (2021), arXiv:2008.04333 [astro-ph.HE] .
  * Chandrasekhar (1970) S. Chandrasekhar, Astrophys. J. 161, 561 (1970).
  * Friedman and Schutz (1978) J. L. Friedman and B. F. Schutz, Astrophys. J. 221, 937 (1978).
  * East  _et al._ (2016) W. E. East, V. Paschalidis, F. Pretorius, and S. L. Shapiro, Phys. Rev. D93, 024011 (2016), arXiv:1511.01093 [astro-ph.HE] .
  * Bernuzzi  _et al._ (2008) S. Bernuzzi, A. Nagar, and R. De Pietri, Phys. Rev. D77, 044042 (2008), arXiv:0801.2090 [gr-qc] .
  * Berti  _et al._ (2009) E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Berti and Klein (2014) E. Berti and A. Klein, Phys. Rev. D90, 064012 (2014), arXiv:1408.1860 [gr-qc] .
  * Calderón Bustillo  _et al._ (2021) J. Calderón Bustillo, S. H. W. Leong, T. Dietrich, and P. D. Lasky, Astrophys. J. Lett. 912, L10 (2021), arXiv:2006.11525 [gr-qc] .
  * Bernuzzi  _et al._ (2012) S. Bernuzzi, M. Thierfelder, and B. Brügmann, Phys.Rev. D85, 104030 (2012), arXiv:1109.3611 [gr-qc] .
  * Dietrich  _et al._ (2018a) T. Dietrich, S. Bernuzzi, B. Brügmann, and W. Tichy, in  _2018 26th Euromicro International Conference on Parallel, Distributed and Network-based Processing (PDP)_ (2018) pp. 682–689, arXiv:1803.07965 [gr-qc] .
  * Bernuzzi and Dietrich (2016) S. Bernuzzi and T. Dietrich, Phys. Rev. D94, 064062 (2016), arXiv:1604.07999 [gr-qc] .
  * Fujimoto  _et al._ (2022) Y. Fujimoto, K. Fukushima, K. Hotokezaka, and K. Kyutoku, (2022), arXiv:2205.03882 [astro-ph.HE] .
  * Dietrich  _et al._ (2018b) T. Dietrich, D. Radice, S. Bernuzzi, F. Zappa, A. Perego, B. Brügmann, S. V. Chaurasia, R. Dudi, W. Tichy, and M. Ujevic, Class. Quant. Grav. 35, 24LT01 (2018b), arXiv:1806.01625 [gr-qc] .
  * (86) http://www.computational-relativity.org/, Computational Relativity.
  * Kawaguchi  _et al._ (2018) K. Kawaguchi, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D97, 044044 (2018), arXiv:1802.06518 [gr-qc] .
  * Kiuchi  _et al._ (2017) K. Kiuchi, K. Kawaguchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D96, 084060 (2017), arXiv:1708.08926 [astro-ph.HE] .
  * Kiuchi  _et al._ (2019) K. Kiuchi, K. Kyohei, K. Kyutoku, Y. Sekiguchi, and M. Shibata, (2019), arXiv:1907.03790 [astro-ph.HE] .
  * Camilletti  _et al._ (2022) A. Camilletti, L. Chiesa, G. Ricigliano, A. Perego, L. C. Lippold, S. Padamata, S. Bernuzzi, D. Radice, D. Logoteta, and F. M. Guercilena, (2022), arXiv:2204.05336 [astro-ph.HE] .
  * Brügmann  _et al._ (2008) B. Brügmann, J. A. Gonzalez, M. Hannam, S. Husa, U. Sperhake, _et al._ , Phys.Rev. D77, 024027 (2008), arXiv:gr-qc/0610128 [gr-qc] .
  * Thierfelder  _et al._ (2011) M. Thierfelder, S. Bernuzzi, and B. Brügmann, Phys.Rev. D84, 044012 (2011), arXiv:1104.4751 [gr-qc] .
  * Radice and Rezzolla (2012) D. Radice and L. Rezzolla, Astron. Astrophys. 547, A26 (2012), arXiv:1206.6502 [astro-ph.IM] .
  * Banik  _et al._ (2014) S. Banik, M. Hempel, and D. Bandyopadhyay, Astrophys. J. Suppl. 214, 22 (2014), arXiv:1404.6173 [astro-ph.HE] .
  * Bombaci and Logoteta (2018) I. Bombaci and D. Logoteta, Astron. Astrophys. 609, A128 (2018), arXiv:1805.11846 [astro-ph.HE] .
  * Logoteta  _et al._ (2020) D. Logoteta, A. Perego, and I. Bombaci,  (2020), 10.1051/0004-6361/202039457, arXiv:2012.03599 [nucl-th] .
  * Typel  _et al._ (2010) S. Typel, G. Ropke, T. Klahn, D. Blaschke, and H. H. Wolter, Phys. Rev. C81, 015803 (2010), arXiv:0908.2344 [nucl-th] .
  * Hempel  _et al._ (2012) M. Hempel, T. Fischer, J. Schaffner-Bielich, and M. Liebendorfer, Astrophys. J. 748, 70 (2012), arXiv:1108.0848 [astro-ph.HE] .
  * Lattimer and Swesty (1991) J. M. Lattimer and F. D. Swesty, Nucl. Phys. A535, 331 (1991).
  * Steiner  _et al._ (2013) A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astrophys. J. 765, L5 (2013), arXiv:1205.6871 [nucl-th] .
  * Schneider  _et al._ (2017) A. S. Schneider, L. F. Roberts, and C. D. Ott, Phys. Rev. C96, 065802 (2017), arXiv:1707.01527 [astro-ph.HE] .
  * Alford  _et al._ (2005) M. Alford, M. Braby, M. W. Paris, and S. Reddy, Astrophys. J. 629, 969 (2005), arXiv:nucl-th/0411016 .
  * Engvik  _et al._ (1996) L. Engvik, G. Bao, M. Hjorth-Jensen, E. Osnes, and E. Ostgaard, Astrophys. J. 469, 794 (1996), arXiv:nucl-th/9509016 .
  * Müther  _et al._ (1987) H. Müther, M. Prakash, and T. L. Ainsworth, Phys. Lett. B 199, 469 (1987).
  * Müller and Serot (1996) H. Müller and B. D. Serot, Nucl. Phys. A606, 508 (1996), arXiv:nucl-th/9603037 [nucl-th] .
  * Douchin and Haensel (2001) F. Douchin and P. Haensel, Astron. Astrophys. 380, 151 (2001), astro-ph/0111092 .
  * Lackey  _et al._ (2006) B. D. Lackey, M. Nayyar, and B. J. Owen, Phys. Rev. D 73, 024021 (2006), arXiv:astro-ph/0507312 .
  * Read  _et al._ (2009) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D79, 124032 (2009), arXiv:0812.2163 [astro-ph] .
  * Nagar  _et al._ (2019) A. Nagar, G. Pratten, G. Riemenschneider, and R. Gamba, (2019), arXiv:1904.09550 [gr-qc] .
  * Bernuzzi  _et al._ (2014b) S. Bernuzzi, A. Nagar, S. Balmelli, T. Dietrich, and M. Ujevic, Phys.Rev.Lett. 112, 201101 (2014b), arXiv:1402.6244 [gr-qc] .
  * Zappa  _et al._ (2019) F. Zappa, S. Bernuzzi, F. Pannarale, M. Mapelli, and N. Giacobbo, Phys. Rev. Lett. 123, 041102 (2019), arXiv:1903.11622 [gr-qc] .
  * Radice (2017) D. Radice, Astrophys. J. 838, L2 (2017), arXiv:1703.02046 [astro-ph.HE] .
  * Vitale  _et al._ (2012) S. Vitale, W. Del Pozzo, T. G. Li, C. Van Den Broeck, I. Mandel, _et al._ , Phys.Rev. D85, 064034 (2012), arXiv:1111.3044 [gr-qc] .
  * Storn and Price (1997) R. Storn and K. Price, J. of Global Optimization 11, 341–359 (1997).
  * Damour  _et al._ (2011) T. Damour, A. Nagar, and M. Trias, Phys. Rev. D83, 024006 (2011), arXiv:1009.5998 [gr-qc] .
  * Chatziioannou  _et al._ (2017b) K. Chatziioannou, A. Klein, N. Yunes, and N. Cornish, Phys. Rev. D95, 104004 (2017b), arXiv:1703.03967 [gr-qc] .
  * Gamba  _et al._ (2021b) R. Gamba, M. Breschi, S. Bernuzzi, M. Agathos, and A. Nagar, Phys. Rev. D 103, 124015 (2021b), arXiv:2009.08467 [gr-qc] .
  * Breschi  _et al._ (2021b) M. Breschi, R. Gamba, and S. Bernuzzi, Phys. Rev. D 104, 042001 (2021b), arXiv:2102.00017 [gr-qc] .
  * Sekiguchi  _et al._ (2011) Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Phys.Rev.Lett. 107, 211101 (2011), arXiv:1110.4442 [astro-ph.HE] .
  * Huang  _et al._ (2022) Y.-J. Huang, L. Baiotti, T. Kojo, K. Takami, H. Sotani, H. Togashi, T. Hatsuda, S. Nagataki, and Y.-Z. Fan, (2022), arXiv:2203.04528 [astro-ph.HE] .
  * Chatziioannou and Han (2020) K. Chatziioannou and S. Han, Phys. Rev. D 101, 044019 (2020), arXiv:1911.07091 [gr-qc] .
  * Most  _et al._ (2020) E. R. Most, L. Jens Papenfort, V. Dexheimer, M. Hanauske, H. Stoecker, and L. Rezzolla, Eur. Phys. J. A 56, 59 (2020), arXiv:1910.13893 [astro-ph.HE] .
  * Raithel and Most (2022) C. A. Raithel and E. R. Most, (2022), arXiv:2201.03594 [astro-ph.HE] .
  * Kastaun  _et al._ (2016) W. Kastaun, R. Ciolfi, and B. Giacomazzo, Phys. Rev. D 94, 044060 (2016), arXiv:1607.02186 [astro-ph.HE] .
  * Kiuchi  _et al._ (2018) K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Phys. Rev. D97, 124039 (2018), arXiv:1710.01311 [astro-ph.HE] .
  * Radice  _et al._ (2016b) D. Radice, F. Galeazzi, J. Lippuner, L. F. Roberts, C. D. Ott, and L. Rezzolla, Mon. Not. Roy. Astron. Soc. 460, 3255 (2016b), arXiv:1601.02426 [astro-ph.HE] .