Post-merger Gravitational-wave Signal from Neutron-star Binaries: A New Look at an Old Problem

The spectral properties of the post-merger gravitational-wave signal from a binary of neutron stars encodes a variety of information about the features of the system and of the equation of state describing matter around and above nuclear saturation density. Characterizing the properties of such a signal is an “old” problem, which first emerged when a number of frequencies were shown to be related to the properties of the binary through “quasi-universal” relations. Here we take a new look at this old problem by computing the properties of the signal in terms of the Weyl scalar ψ 4. In this way, and using a database of more than 100 simulations, we provide the first evidence for a new instantaneous frequency, f0ψ4 , associated with the instant of quasi-time-symmetry in the dynamics, and which also follows a quasi-universal relation. We also derive a new quasi-universal relation for the merger frequency fmerh , which provides a description of the data that is 4 times more accurate than previous expressions while requiring fewer fitting coefficients. Finally, consistent with the findings of numerous studies before ours, and using an enlarged ensemble of binary systems, we point out that the ℓ = 2, m = 1 gravitational-wave mode could become comparable with the traditional ℓ = 2, m = 2 mode on sufficiently long timescales, with strain amplitudes in a ratio ∣h 21∣/∣h 22∣ ∼ 0.1–1 under generic orientations of the binary, which could be measured by present detectors for signals with a large signal-to-noise ratio or by third-generation detectors for generic signals should no collapse occur.


INTRODUCTION
The observation of a gravitational-wave (GW) signal from the binary neutron-star (BNS) merger event GW170817 (The LIGO Scientific Collaboration & The Virgo Collaboration 2017), and the detection of an electromagnetic (EM) counterpart, has testified the enormous potential of GW astronomy.Starting from early works with simplified equations of state (EOSs) (see, e.g., Shibata et al. 2005;Anderson et al. 2008;Liu et al. 2008;Baiotti et al. 2008;Hotokezaka et al. 2011), increasingly more comprehensive simulations of these events, which involve an ever more detailed description of the microphysics (Bauswein et al. 2019;De Pietri et al. 2019;Gieg et al. 2019;Tootle et al. 2022;Most et al. 2022;Camilletti et al. 2022;Ujevic et al. 2023), of the magnetic-field evolution (Rezzolla et al. 2011;Dionysopoulou et al. 2013;Ciolfi et al. 2019;Sun et al. 2022;Zappa et al. 2023), and its amplification (Kiuchi et al. 2015;Palenzuela et al. 2022;Chabanov et al. 2023), and of transport of neutrinos (Foucart et al. 2022;Zappa et al. 2023), allow one to make predictions from the early inspiral up to the long-term evolution of the postmerger remnant (De Pietri et al. 2020;Kiuchi et al. 2022).During each stage in the evolution of the binary, the features of the GW and EM signals change in a characteristic manner, encoding information on the properties of the constituent neutron stars and of the hypermassive neutron star (HMNS) produced after the merger and, hence, on the governing EOS.
Characterising the properties of the post-merger GW signal is a rather "old" problem, which has first emerged when a number of peculiar frequencies were shown to be related with the properties of the binary through quasi-universal relations, i.e., relations that are almost independent of the specific EOS.These relations have been suggested for the GW frequency at merger f mer (Read et al. 2013;Bernuzzi et al. 2014;Takami et al. 2015;Rezzolla & Takami 2016;Most et al. 2019;Bauswein et al. 2019;Weih et al. 2020;Gonzalez et al. 2022), the dominant frequency in the postmerger spectrum f 2 (see, e.g., Oechslin & Janka 2007;Bauswein & Janka 2012;Read et al. 2013;Rezzolla & Takami 2016;Gonzalez et al. 2022), and other frequencies identifiable in the transient period right after the merger (Bauswein & Stergioulas 2015;Takami et al. 2015;Rezzolla & Takami 2016).Fits to these quasi-universal relations have been employed in a number of studies (see, e.g., (Bauswein et al. 2016;Baiotti & Rezzolla 2017) for some reviews).These EOSinsensitive relations can help enormously in constraining the EOS of matter at nuclear densities, marking the possible appearance of phase transitions (Most et al. 2019;Weih et al. 2020;Liebling et al. 2021;Prakash et al. 2021;Fujimoto et al. 2023;Tootle et al. 2022;Espino et al. 2023), inform waveform models (Bose et al. 2018;Breschi et al. 2019); however, see Raithel & Most (2022) for possible violations of these relations.
Another relatively "old" problem in the characterisation of the GW signal from BNS is the one about the relative weight of the lower-order multipole ℓ = 2, m = 1.Numerical simulations have highlighted that the HMNS can be subject to a nonaxisymmetric instability that powers the growth of ℓ = 2, m = 1 mode of the rest-mass density distribution and, hence, of the corresponding GW signal (see, e.g., East et al. 2015;Lehner et al. 2016a,b;Radice et al. 2016;East et al. 2019;Papenfort et al. 2022).This mode can already be seeded by the initial asymmetry of the system in the unequalmass case, or develop by the shearing of the contact layers of the binary constituents upon merger.While the ℓ = 2, m = 2 GW mode is the primary contributor to the GW signal, it is damped faster than the other modes leading to interesting secular behaviours.
We here take a new look at both of these old problems by considering the spectral properties of the GW signal when computed in terms of the Weyl scalar ψ 4 .In this way, we are able to find three novel features that can enrich our understanding of the GW signal from BNS mergers.In particular, we first highlight the presence of a new instantaneous frequency, which we dub f ψ4 0 , that can be associated with the instant of quasi time-symmetry in the postmerger dynamics.Interestingly, we find that a quasi-universal relation exists for f ψ4 0 as a function of the tidal deformability κ T 2 and of the binary mass ratio q.Second, by employing a large number of BNS simulations, some of which are taken from the CoRe database (Gonzalez et al. 2022), we obtain a new quasi-universal relation for f mer as a function of κ T 2 and q that not only requires a smaller number of coefficients, but also provides a more accurate description of the data.Finally, as already suggested in (Papenfort et al. 2022), we provide evidence that the ℓ = 2, m = 1 GW mode could become the most powerful mode on secular timescales after the merger.

NUMERICAL AND PHYSICAL FRAMEWORK
Our analysis is based on the GW signal computed via numerical simulations of BNS mergers in full general relativity computed with the codes described in (Radice et al. 2014a,b;Most et al. 2019a,b;Papenfort et al. 2021;Tootle et al. 2021) and using a number of different EOSs (see below).In addition, we employ part of the data contained in the CoRe database (Gonzalez et al. 2022), from where we select only simulations with the highest-resolution.The combined data of 118 irrotational binaries covers the range q : [2.4, 3.33] M ⊙ in the total ADM mass at infinite separation, and κ T 2 ∈ [33,458] in the tidal deformability.The dataset comprises a variety of EOSs including some with quark matter (Prakash et al. 2021;Logoteta, Domenico et al. 2021;Alford et al. 2005;Demircik et al. 2022;Tootle et al. 2022).
A crucial role in our analysis is played by the use of the Weyl scalar ψ 4 in place of the standard dimensionless strain polarisations h +,× .The two quantities are mathematically equivalent and related by two time derivatives (i.e., ψ 4 = ∂ 2 t (h + − ih × ); see (Bishop & Rezzolla 2016) for a review).However, while ψ 4 is computed from the simulations, h +,× are obtained after a nontrivial double time integration (the transformation from h +,× to is ψ 4 trivial as it involves derivatives and not integrals; see (Calderon Bustillo et al. 2022) for a data-analysis framework based on ψ 4 , which can obviously be employed for all types of compact-object binaries).More importantly, the evolution of the GW frequency from ψ 4 is less rapid than from the strain, i.e., ∂ t ln f ψ4 GW (t) ≪ ∂ t ln f h GW (t), thus making it easier and more robust to characterise the features of the ψ 4 GW signal.In this sense, while ψ 4 and h +,× are related by simple time derivatives, the analysis carried out with the former does provide additional information as it allows for the determination of properties that are harder to capture with the latter.

OLD AND NEW FREQUENCIES
Figure 1 reports the complete information of the GW signal from a representative binary in our sample.Using a 3D representation, we report on the left the ℓ = 2, m = 2 mode of the GW signal ψ 4 (t) (light red) and its amplitude |ψ 4 (t)| (dark red), the instantaneous frequency f ψ4 GW (t) (black), and the power spectral density (PSD) √ 2f ψ4 (blue) as a function of the frequency f (see Rezzolla & Takami 2016, for details on the definition).Also indicated are the three main frequencies in our analysis: the frequency at merger f ψ4 mer , i.e., the GW frequency at the first maximum of |ψ 4 |, the frequency at quasi time-symmetry f ψ4 0 , i.e., the GW frequency at the first minimum of |ψ 4 |, and the dominant frequency of the HMNS emission f ψ4 2 .To help the eye, we also mark with lines the corresponding times t ψ4 mer (dashed), t ψ4 0 (dotted), and frequencies (dashed, dotted and dot-dashed respectively).The right panel of Fig. 1 shows the same quantities but when computed from the strain.By comparing the black lines in the left and right panels it is straightforward to realise that the variation of f h GW (t) is much larger than that in f ψ4 GW (t) over the same interval of ∼ 1 ms after the merger.It is this very rapid change in f h GW (t) that makes the identification of f h 0 extremely difficult, if not impossible.Note also that while 3D representation of the complete information in the GW signal from a representative binary (Tootle et al. 2022).The left panel shows the ℓ = 2, m = 2 mode of the GW signal ψ4(t) (light red) and its amplitude |ψ4(t)| (dark red), the instantaneous frequency f GW (t) (black), and the power spectral density (PSD) √ 2f ψ4 (blue) as a function of the frequency f .Also indicated are the frequency at merger f ψ 4 mer , the frequency at quasi time-symmetry f ψ 4 0 , the dominant frequency of the HMNS emission f2, and the corresponding times where these frequencies appear.The right panel shows the same quantities but when computed in terms of the GW strain.
in both representations f mer < f 2 < f 0 , the numerical values of the various quantities are similar but not identical.However, f ψ4 2 ≃ f h 2 to very good precision (the largest differences are ≲ 4%) simply because this frequency is relative to a mostly monochromatic GW signal; hence, hereafter we simply assume f ψ4 2 = f h 2 =: f 2 .Finally, because in all representations f 0 is largest frequency measured, even a crude measure of largest frequency in the signal will serve as a first estimate of the f 0 frequency. 1 Besides marking the time of the first amplitude minimum, from a physical point of view t ψ4 0 corresponds to the time when the two stellar cores have reached the minimum separation and are about to bounce-off each other.At this instant, the corresponding amplitude of ψ 4 shows a clear minimum, while the instantaneous GW frequency a local maximum [the discussion in the Supplemental Material (SM) illustrates this behaviour very clearly by employing the toy model introduced in (Takami et al. 2015)].

QUASI-UNIVERSAL RELATIONS
We next proceed to the derivation of quasi-universal relations that can be employed to deduce the physical properties of the binary.Following the approach started already in (Takami et al. 2014(Takami et al. , 2015;;Rezzolla & Takami 2016), which captures the logarithmic variation of a properly rescaled mass and frequency, we express the relevant frequencies in terms of a power expansion of the mass ratio q, i.e., log 10 n , where f is any of the frequencies we consider 1 From a numerical point of view, we note that the f 0 frequency is always below ∼ 4 kHz and is therefore much smaller than the typical sampling frequency of the ψ 4 scalar, that is ≃ 80−100 kHz. (i.e., Hereafter, we will refer to this generic fitting functions as F 1 . Figure 2 provides a 3D representation of the measured GW frequencies f ψ4 mer and f ψ4 0 as a function of κ T 2 and q (see also the SM for fits to f h mer and f 2 ).Also reported is the fitting surface described by F 1 , with the best-fit parameters listed in Table 2 of the SM for all the frequencies considered.Furthermore, for each frequency we report below the relative error of the fit in the two principal directions of the fit, κ T 2 and q.Despite their simple form, our fits for f ψ4 mer and f ψ4 0 capture the data very well, showing average relative errors that are ≲ 1% and maximal relative errors ≲ 2% for other than equal-mass binaries.
It is interesting to compare our functional fitting form F 1 for f h mer , which needs only five fitting coefficients, with the one proposed in (Breschi et al. 2022) for irrotational binaries, which we will refer to as F 2 , and that requires twice as many coefficients.In order to compare F 1 and F 2 it is first necessary to distinguish the "pipeline", that is, the technical procedure employed to extract the frequencies from the data.We thus indicate with P 1 the pipeline discussed above and with P 2 that released in (Gonzalez et al. 2022).Naturally, each fitting function can be applied to either pipeline, so that F 1 (P 1 ) indicates the use of our fitting form to data computed with our pipeline.In Fig. 3, we present the relative differences between the measured frequencies for the 118 binaries considered and the corresponding values from the fit, with different rows referring to the four possibilities.
Overall, the comparison in Fig. 3 shows that F 1 leads to smaller relative errors with a maximum residual error of ∼ 2% and an an average residual error that is between two and four times smaller than for F 2 .As a cautionary note we as measured from the data (coloured circles) and presented as a function of κ T 2 and q.Also reported are the best-fit surfaces, while shown below are the relative errors of the fit in the two principal directions.Stars mark binaries modelled with the V-QCD EOS and thus having a strong first-order phase transition (Demircik et al. 2022;Tootle et al. 2022).should remark that we have specialised the fitting F 2 , which is more general and can include spinning and eccentric binaries, to the case relevant for this comparison, namely, irrotational binaries.Hence, our conclusions apply only to such binaries.

SECULAR GW EMISSION
The last point we cover regards the relative strengths of the ℓ = 2, m = 2 and ℓ = 2, m = 1 GW modes.The importance of the latter was first pointed out in (Paschalidis et al. 2015;Lehner et al. 2016a) and is produced by corresponding asymmetries in the rest-mass density.The emergence of an m = 1 deformation is well-known to occur in isolated stars (Chandrasekhar 1969;Shibata et al. 2000;Baiotti et al. 2007;Franci et al. 2013;Löffler et al. 2015) that have a sufficiently large amount of rotational kinetic energy T and emerges when the ratio T /|W |, where W is the gravitational binding energy, exceeds a certain threshold (in a systematic analysis, Baiotti et al. 2007, have shown that this happens for T /|W | ≳ 0.25).In such stars, the m = 1 mode in the restmass density would grow exponentially reaching equipartition with the m = 2 mode, and, subsequently, represent the largest deformation.A similar phenomenology seems to be present also for the postmerger remnant, as already hinted in (Papenfort et al. 2022), but as shown more clearly by the numerous binaries considered here.Figure 4  | = 1 within the simulated time, the large majority exhibits a trend that we try to capture by extrapolating linearly in time after averaging the last 5 ms of the evolution.Using a colour code to distinguish binaries with different q, it becomes clear that the initial strength of the m = 1 mode is inversely proportional to the mass ratio, so that for an extremely asymmetric binary, i.e., q ≲ 0.6, |ψ 21 4 |/|ψ 22 4 | can be more than two orders of magnitude larger than for equal-mass binaries.At the same time, the initial mode-amplitude ratio does not depend on the merger dimensionless spin (Papenfort et al. 2022).Unsurprisingly, a similar behaviour can be observed when  |, with the left panel showing a selected set of binaries to highlight the behaviour for different mass ratios (colour code) and the right one reporting all of the binaries.Black dotted horizontal lines are used to mark the position when the two modes have equal amplitudes, while straight coloured lines report a linear extrapolation after averaging the last 5 ms of the evolution; note that the large majority of binaries exhibits a growing trend that is maintained throughout the computed evolution.
computing the signal-to-noise (SNR) ratio in the m = 2 and m = 1 modes.Specifically, by computing a time-windowed SNR ratio we can quantify the growing contribution of the subdominant mode in a similar fashion to estimates given by (Lehner et al. 2016b) and find that its increases to ≃ O(1) if the m = 1 is not suppressed (see SM for details).
If confirmed by systematic, long-term evolutions, this finding would change the standard picture in which the largest signal in the BNS postmerger is to be expected at the f 2 frequency.Rather, the trend reported here suggests that the most powerful feature in the PSD for long-lived HMNSs may actually appear at a frequency ≃ 1 2 f 2 .Because this falls in a more favourable region of the detectors sensitivities, and assuming a detection angle not favouring either of the modes, the corresponding signal-to-noise ratio will grow proportionally to the ratio between detectors noise at 1 2 f 2 and f 2 , hence by a factor of ∼ 2 for LIGO or Virgo.
While promising, this prospect should be accompanied by some caveats.First, it is possible that the growth rate may be weaker than the one estimated here.Second, the extrapolation assumes that the HMNS will not collapse to a black hole before reaching |ψ 21 4 |/|ψ 22 4 | ∼ 1 and while this is likely for soft EOSs and low-mass binaries, it may not happen if the EOS is stiff and the binary massive.Third, all binaries in our sample have zero deformation.A robust conclusion that can be inferred from the results shown in Fig. 4 is that remnants with a long lifetime, as it was likely the case for GW170817 (Rezzolla et al. 2018;Gill et al. 2019;Murguia-Berthier et al. 2021), will reasonably have the ℓ = 2, m = 1 as the least-damped mode.Hence, considerable spectral power should be present at frequencies 1 2 f 2 and f 2 , with the main strain amplitudes in a ratio |h 21 |/|h 22 | ∼ 0.1 − 1 for generic orientations (e.g., for an inclination of 2 arctan(1/2) ∼ 53 • two modes have the same spin-weighted spherical-harmonics coefficients).

CONCLUSION
Leveraging on a rich literature developed over the last ten years on this subject, we have considered again the spectral properties of the signal when computed in terms of the Weyl scalar ψ 4 rather than in terms of the GW strain h +,× .Exploiting the better behaviour of ψ 4 , we were able to highlight three novel features that can be used to better infer physical information from the detected signal.
First, by employing a large number of simulations spanning a considerable set of EOSs and mass ratios, we have shown the existence of a new instantaneous frequency, f ψ4 0 , that can be associated with the instant of quasi timesymmetry in the postmerger dynamics.This corresponds to when the stellar cores in the merger remnant have reached their minimum separation and are about to bounce-off each other.Just like other spectral frequencies of the BNS GW signal, f ψ4 0 also follows a quasi-universal behaviour as a function of the tidal deformability κ T 2 and of the binary mass ratio q, for which we provide a simple and yet accurate analytical expression.Second, we have obtained a new quasiuniversal relation for the merger frequency f h mer as a function of κ T 2 and q.The new expression not only requires a smaller number of fitting coefficients than alternative expressions in the literature, but it also provides a more accurate description of the data, with a residual error that is four times smaller on average.Finally, we have pointed out the evidence that the ℓ = 2, m = 1 could become the most powerful GW mode on sufficiently long timescales, with strain amplitudes for the dominant modes that are in a ratio |h 21 |/|h 22 | ∼ 0.1 − 1.Should this mode not be suppressed by the collapse of the HMNS to a black hole or by other dissipative effects such as magnetic fields, considerable spectral power should be present at frequencies 1 2 f 2 , where it could be detected in conditions of smaller signal-to-noise ratios or by third-generation detectors.
The results presented here can be improved by enlarging the number of BNS simulations considered, by increasing the variance in the microphysical description (e.g., including simulations with magnetic fields and neutrino transport), by performing additional long-term evolutions, and by extending the fitting approach to binaries with spins and eccentricity.We will explore these extensions in future work.Data policy.The relevant data that supports the findings of this paper is available from the first author and can be shared upon a reasonable request.

TOY-MODEL ANALOGY
Although it is possible to characterise the new spectral feature f 0 simply in terms of the instantaneous GW frequency at the time t ψ4 0 at which the Weyl scalar has its first minimum, it is more interesting to associate with this definition also a physical interpretation.As mentioned in the main text, t ψ4 0 effectively corresponds to when the two stellar cores have reached the minimum separation and are about to bounce-off each other.At this instant, the radial velocity of the two stellar cores and the angular velocity of the HMNS has a minimum, so that the corresponding amplitude of ψ 4 has a minimum and instantaneous GW frequency a local maximum.
To illustrate this behaviour, we employ the toy model developed in Takami et al. (2015) (and also employed in other studies, e.g., Ellis et al. (2018); Lucca et al. (2021)) to describe the dynamics of the postmerger remnant and which consists of an axisymmetric disk rotating rapidly at a given angular frequency, say Ω 2 , to which two spheres reproducing the two stellar cores are connected but are also free to oscillate via a spring that connects them (see Fig. 17 Takami et al. (2015)).In such a system, the two spheres will either approach each other, decreasing the moment of inertia of the system, or move away from each other, increasing the moment of inertia.Because the total angular momentum is essentially conserved, the system's angular frequency will vary between a minimum value Ω 1 (corresponding to the time when the two spheres are at the largest separation) and a maximum value Ω 3 (corresponding to the time when the two spheres are at the smallest separation).Overall, the mechanical toy model will rotate with an angular frequency that is a function of time and bounded by Ω 1 and Ω 3 , where more time is spent and hence more spectral power is accumulated.If a damping is introduced, the excursion of the oscillations between the spheres will decrease over time and eventually stop; when this happens, the toy model will simply rotate at the frequency Ω 2 , as does the HMNS after the transient period and before other deformations (e.g., the ℓ = 2, m = 1 mode) affect it.
The Lagrangian of the toy model can be found in Takami et al. (2014) and from it it is possible to derive the differential equation for the radial displacement r(t) where M and m are the masses of the disc and of the spheres, R is the radius of the disc, c 1 an integration constant related to the total angular momentum, r 0 the natural displacement of mass and b is responsible for dissipative effects.The dissipation due to the emission of GWs and which results in the two spheres getting closer to one another with time and with a decreasing radial displacement is introduced by varying the natural displacement r 0 via via an exponential function exp(−t/τ ) with a suitably chosen relaxation time τ .Note that this modification carries straight from the Lagrangian to the final ODE without additional time derivatives, as it is present in the derivative of the Lagrangian with respect to the canonical coordinate r(t) and not its conjugate.
The system is sufficiently simple that it is possible to compute the GW amplitude and frequency derived from the quadrupole formula so that the strain components read with d the distance from the source to the detector.The ψ 4 polarisations are obtained by differentiating twice with respect to time.Introducing the abbreviations A := ṙ2 + r(r − 2rΩ 2 ) and B := r(4 ṙΩ + r Ω), as well as recalling that Ω(t) = φ(t) we write them now explicitly where ψ 4,+ := ∂ 2 t h + and ψ 4,× := ∂ 2 t h × .The amplitude and the instantaneous frequency of the signal are calculated using the usual definitions The rather lengthy final expressions read: Figure 5 reports the solution of the system when considering the following set of representative parameters: M = m = 10, k = 0.2, b = 0.25, c 1 = 25, R = 22, r 0 = 5, ṙ0 = 0.01, τ = 1000.More specifically, the left panel reports the GW frequency (top row) as computed from the toy model and from an actual simulation (marked by the index "NR"), the corresponding GW amplitudes (middle row) and the radial displacement r(t) (bottom row).Note how the GW frequencies and amplitudes are anticorrelated and the minimum in the radial displacement at t ψ4 0 corresponds to the maximum frequency, i.e., f ψ4 0 , and minimum amplitude, as expected in a condition of quasi time-symmetry (in the toy model f h 0 is actually a local minimum of f h GW (t), a case often found in NR data).The right panel of Fig. 5, on the other hand, is used to show the relation between the GW quantities of the toy model and other dynamical quantities such as the time derivatives of the angular velocity, ∂ t Ω(t) and ∂ 2 t Ω(t), and of the radial displacement, ∂ t r(t) and ∂ 2 t r(t) (the top row of the right panel contains the same information as in the top and middle row .The same as in Fig. 2 in the main text but for the GW frequencies f h mer and f2.Note how F1 is equally accurate in modelling other frequencies, also when obtained from the strain as in the case for f h mer and f2. of the left panel).It is clear that for the first few oscillations the minima of the amplitude (maxima of the frequency) correspond to moments in time when the angular frequency Ω attains a local maximum, and the radial displacement a local minimum.We find it rather remarkable that the model originally developed to provide justification for the dominant postmerger frequency can be used to mimic the presence of the f 0 frequency with no additional adjustments, appearing as a result of chosen parameters and initial conditions.Finally, while we do not consider it here, it would be useful to generalise the model to admit unequal sphere masses, i.e., m/M < 1, and assess the impact it has on the f 0 frequency and on the evolution of the m = 2 and m = 1 modes.

ON THE EFFECTIVENESS OF THE FIT
To illustrate effectiveness of the fitting functional form F 1 in capturing the spectral properties of the signal in a quasi-universal manner, we report in Fig. 2 the same 3D representation as of the GW frequencies shown in Fig. 2 but for the GW frequencies f h mer and f 2 ; note that f h mer is referred to as f max in Takami et al. (2014); Rezzolla & Takami (2016) and that the f 2 frequency corresponds to the peak of the PSD computed using the definition in Tootle et al. (2022), with the time integration encompassing the full GW signal.Clearly, the fitting function F 1 is equally accurate in modelling other frequencies, quite independently of whether they have been computed from the Weyl scalar or from the GW strain, as it is as it is the case for f h mer and f 2 .We also note that the equal-mass limit of the fitting function F 1 for the frequency f h mer , recovers quite accurately the simpler fit examined in Takami et al. (2014) which expressed the merger frequency as a first-order expression of the tidal deformability f h mer = a 0 + b(κ T 2 ) 1/5 .When setting q = 1 and n = 1/5 in F 1 , we obtain b = b 0 + b 1 + b 2 = −0.199(cf., Table 1), thus resulting in a relative difference that is ≲ 2% with respect to the value b = −0.195found in Takami et al. (2014).At the same time, we also note that the fit of f h mer is overall slightly better than that for f ψ4 mer (see Table 1); while we do not believe this to be statistically very significant, it may be due to the slightly larger range in which f ψ4 mer is measured.The only exception to the remarkably good fit of the frequencies is found in the f ψ4 0 frequency, where a noticeable deviation from the trend is visible for binary systems with equal or very-unequal masses, regardless of the EOS employed (see lower part of the bottom-right panel of Fig. 2).We attribute this decrease in accuracy to the already described difficulties of measuring this frequency using the ψ 4 GW signal, which are even more severe when using |h|.From a physical point of view, the moment of time symmetry corresponds to the time when the non-axisymmetric quadrupolar deformations of the HMNS are at a minimum, which leads to a severely suppressed GW signal (indeed the amplitude is at a minimum).This is particularly severe for binaries with q ≃ 1, where the ℓ = 2, m = 2 deformation is the largest and is significantly suppressed at the moment of time-symmetry,

Figure 3 .
Figure3.Relative differences between the measured f h mer frequencies and the corresponding values from the fit.Different rows refer to the four different possibilities of applying the fitting function Fi to the data pipeline Pi with i = 1, 2; arrows indicate relative differences above 10%.Note how the new fitting function and pipeline, i.e., F1(P1), provide errors that are about four times smaller.
reports the evolution of the ratio of the GW amplitudes in the two modes |ψ 21 4 |/|ψ 22 4 |, with the left panel showing a selected set of binaries and with the right panel reporting all binaries (the raw timeseries are smoothed over a window ∆t = 0.5 ms).While only a few binaries in the sample reach |ψ 21 4 |/|ψ 22 4

Figure 4 .
Figure 4. Evolution of the ratio of the two main GW modes |ψ 21 4 |/|ψ 22 4|, with the left panel showing a selected set of binaries to highlight the behaviour for different mass ratios (colour code) and the right one reporting all of the binaries.Black dotted horizontal lines are used to mark the position when the two modes have equal amplitudes, while straight coloured lines report a linear extrapolation after averaging the last 5 ms of the evolution; note that the large majority of binaries exhibits a growing trend that is maintained throughout the computed evolution.

Figure 5 .
Figure5.Solution of the toy model for a representative set of parameters.The left panel reports the GW frequency (top row) from the toy model and from an actual simulation (marked by the index "NR"), the GW amplitudes (middle row) and the radial displacement r(t) (bottom row).The right panel shows instead the relation between the GW quantities and other dynamical quantities of the toy model, such as the time derivatives of the angular velocity and of the radial displacement.Note that the GW frequencies and amplitudes are anti-correlated and the minimum in radial displacement at t ψ 4 0 corresponds to the maximum frequency f ψ 4 0 and minimum amplitude.Also, the minima of the amplitude (maxima of the frequency) correspond to moments in time when the angular frequency Ω attains a local maximum, and the radial displacement a local minimum.