Modeling the postmerger gravitational wave signal and extracting binary properties from future binary neutron star detections

Gravitational wave astronomy has established its role in measuring the equation of state governing cold supranuclear matter. To date and in the near future, gravitational wave measurements from neutron star binaries are likely to be restricted to the inspiral. However, future upgrades and the next generation of gravitational wave detectors will enable us to detect the gravitational wave signatures emitted after the merger of two stars, at times when densities beyond those in single neutron stars are reached. Therefore, the postmerger gravitational wave signal enables studies of supranuclear matter at its extreme limit. To support this line of research, we present new and updated phenomenological relations between the binary properties and characteristic features of the postmerger evolution. Most notably, we derive an updated relation connecting the mass-weighted tidal deformability and the maximum neutron star mass to the dominant emission frequency of the postmerger spectrum. With the help of a configuration-independent Bayesian analysis using simplified Lorentzian model functions, we find that the main emission frequency of the postmerger remnant, for signal-to-noise ratios of $8$ and above, can be extracted within a 1-sigma uncertainty of about 100 Hz for Advanced LIGO and Advanced Virgo's design sensitivities. In some cases, even a postmerger signal-to-noise ratio of $4$ can be sufficient to determine the main emission frequency. This will enable us to measure binary and equation-of-state properties from the postmerger, to perform a consistency check between different parts of the binary neutron star coalescence, and to put our physical interpretation of neutron star mergers to the test.


I. INTRODUCTION
The extreme densities and conditions inside neutron stars (NSs) cannot be reached in existing experiments.This makes NSs a unique laboratory to study the equation of state (EOS) governing cold-supranuclear dense material.Following the first detection of a gravitational wave (GW) signal originating from the coalescence of a binary neutron star (BNS) system, GW170817, by the Advanced LIGO [1] and Advanced Virgo detectors [2], it became possible to constrain the NS EOS by analyzing the measured GWs [3][4][5][6][7].Because of the increasing sensitivity of GW interferometers, multiple detections of merging BNSs are expected in the near future [8].This will make GW astronomy an inevitable tool within nuclear physics community.
In general, there are two ways to extract information about the EOS governing the NS's interior from a GW detection.The first method relies on the modeling of the BNS inspiral [9][10][11][12][13] and on waveform approximants that include tidal effects, represent accurately the system's properties, and are of sufficiently low computational cost that they can be used in parameter estimation pipelines, e.g., [10,14,15].The zero-temperature EOS is then constrained by measuring a mass-weighted combination of the quadrupolar tidal deformability Λ or similar parameters that characterize tidal interactions, e.g., [16,17].
To date, the advanced GW detectors have only been able to observe the inspiral of the two NSs [39,40] and no postmerger signal has been observed.This nonobservation is caused by the higher emission frequency at which current GW detectors are less sensitive.But, the increasing sensitivity of the 2nd generation of GW detectors (Advanced LIGO and Advanced Virgo) will not only increase the detection rate of BNS inspiral signals, there will also be the chance of observing the postmerger signal for a few 'loud' events.Ref. [41] finds that for sources similar to GW170817 but observed with Advanced LIGO and Advanced Virgo's design sensitivities, the postmerger part of the BNS coalescence might have an SNR of ∼ 2-3.The planned third generation of GW interferometers, e.g., the Einstein Telescope [42][43][44] or the Cosmic Explorer [45], have the capability to detect the postmerger signal of upcoming BNS mergers with SNRs up to ∼ 10.
Unfortunately, the postmerger spectrum is influenced in a complicated way by thermal effects, magnetohydrodynamical instabilities, neutrino emissions, phase transitions, and dissipative processes, e.g., [24,25,[46][47][48][49][50].Currently, any postmerger study relies heavily on expensive numerical relativity (NR) simulations and there is to date no possibility to perform simulations incorporating all necessary microphysical processes.Therefore, our current theoretical understanding of this part of the BNS coalescence is overall limited.In addition, there has been no NR simulation yet which has been able to show convergence of the GW phase in the postmerger.While this observation can be generally explained by the presence of shocks or discontinuities formed during the collision of the two stars, it also increases our uncertainty on any quantitative result.
Nevertheless, the community tried to construct postmerger approximants focusing on characteristic (robust) features present in NR simulations.The discovery of quasi-universal relations is a building block for most descriptions of the postmerger GW spectrum.Clark et al. [51] showed that a principle component analysis can be used to reduce the dimensionality of the spectrum for equal mass binaries once the different spectra are normalized and aligned such that the main emission frequencies coincide.Effort has also been put to model the plus polarization in time domain using a superposition of damped sinusoids incorporating quasi-universal relations [35].Relying on a very accurate f 2 estimate for an accurate rescaling of the waveforms Easter et al. [37] created a hierarchical model to estimate the postmerger spectra.
Here, we follow a similar path and try to describe the GW spectrum with a set of a three-and a six-parameter model function with a Lorentzian-like shape.Comparing our ansatz with a set of 54 NR simulations, we find average mismatches of 0.18 for the three-parameter and 0.15 for the six-parameter model; cf.Tab.I. Our approximants do not incorporate directly quasi-universal relations, but are constructed to describe generic postmerger waveforms.Thus, our analysis is flexible and allows to describe almost arbitrary configurations.Employing our model in standard parameter estimation pipelines [52,53] of the LIGO and Virgo Collaborations, we find that we can extract the dominant emission frequency in the postmerger for a number of tests.To our knowledge, this is the first time a model-based (but configuration independent) method is employed within a Bayesian analysis of the postmerger signal.
Once the individual parameters describing the postmerger spectra are extracted, we use fits for the peak frequency to connect the measured signal to the properties of the supranuclear EOS and the merging binary.This way, one can combine measurements from the inspiral and postmerger phase to provide a consistency test for our supranuclear matter description.
Although not used here, we want to mention an alternative approach, which employs the morphology-independent burst search algorithm called BayesWave [54,55].Ref. [36] showed that this approach is capable of reconstructing the postmerger signal and allows to extract properties from the measured GW signal.Even for a measured postmerger SNR of ∼ 5, the main emission frequency of the remnant could be determined within a few dozens of Hz.Compared to BayesWave, our simple model functions might have the advantage that without any modifications of the current code for statistical inferences, in particular the LALInference module [53] available in the LSC Algorithm Library (LAL) Suite, they can be added to existing frequency domain inspiral-merger waveforms describing the first part of the BNS coalescence, e.g., [10,14,15,[56][57][58], to construct a full inspiral-merger-postmerger (IMP) waveform directly employable for GW analysis.Such an IMP study can also be carried out within the BayesWave approach, but seems technically harder since one has to combine model-based and non-model-based algorithms.
Our paper is structured as follows.In Sec.II, we discuss the general time domain and frequency domain morphologies of the postmerger signal as obtained from NR simulations.Based on this discussion, we derive new quasi-universal relations for the time between the merger and the first time domain amplitude minimum, and the first time domain amplitude maximum.We also extend existing quasi-universal relation for the main emission frequency f 2 of the GW postmerger spectrum and its amplitude in the frequency domain.In Sec.III A we discuss two different Lorentzian model functions and their performance to model NR simulations.In Sec.III B a full Bayesian analysis of a set of NR model waveforms is performed.In Sec.III C we show how our analysis can be used to constrain the EOS and how to test consistency between the inspiral and postmerger.We conclude in Sec.IV.We list in the appendix the NR data employed for the construction of the quasi-universal relations presented in the main text.
Unless otherwise stated, this paper uses geometric units by setting G = c = 1.Throughout the work, we employ the NR simulations published in the Computational Relativity (CoRe) database [59].In addition, where explicitly mentioned, we increase our dataset by adding results published in [31,33].We refer the reader to Tab.I for further details about the individual data.

A. Time Domain
While the inspiral GW signal is characterized by a chirp, i.e., a monotonic increase of the GW amplitude and frequency, the postmerger emission shows a nonmonotonic amplitude and frequency evolution.Figure 1 presents one example of a possible postmerger waveform.In the following, we highlight some of the important features characterizing the signal.
First postmerger minimum: By definition, the inspiral ends at the peak of the GW amplitude (merger) marked with a red circle in Fig. 1.After the merger, the amplitude decreases showing a clear minimum (red squared marker) shortly afterwards, see [61,62] for further discussions.Around this intermediate and highly non-linear regime, different frequencies are excited for a few milliseconds, see e.g., [21,34] for further details.While it was already known that the merger frequency can be expressed by a quasi-universal relation, e.g., [20,63], we find that the time between merger and this amplitude minimum also follows a similar relation.
In Fig. 2 we show the time between the merger and the first amplitude minimum, ∆t min /M , as a function of the mass-weighted tidal deformability Λ = 16 13 5  . ( with the individual dimensionless tidal deformabilities Λ = 2 3 k 2 ( R M ) 5 , where k 2 labels the dimensionless = 2 Love number and R labels the radius of the isolated NSs.We show with different colors the mass ratio of each setup defined as q = M A /M B 1, cf. colorbar of Fig. 2.
We find a clear correlation between ∆t min /M and the mass-weighted tidal deformability Λ.A good phenomenological representation is given by with the parameters α = 2.4681 × 10 1 , β = 2.8477 × 10 −3 , γ = 6.6798 × 10 −4 obtained by a least-square fit for which the root-mean-square (RMS) error is 2.4608.Interestingly, the two highest mass ratio simulations do not follow Eq. ( 2).This is caused by the different postmerger evolution for these high-mass ratio setups.While the amplitude minimum is produced when the two NS cores approach each other and potentially get repelled, configurations with very high mass ratio show almost a disruption during the merger, i.e., the lower massive NS deforms significantly under the strong external gravitational field of its companion.
One possible application for the quasi-universal relation for ∆t min /M is the improvement of BNS waveform approximants, i.e., it might help to determine the amplitude evolution after the merger of the two NSs.In particular, incorporating an amplitude tapering after the merger with a width of ∆t min provides a natural ending condition for inspiral-only approximants, e.g., NRTidal [14,58,64] or tidal effective-one-body models [65][66][67].Therefore, Eq. ( 2) might become a central criterion to connect inspiral and postmerger models.
First postmerger maximum: After the minimum of the GW amplitude, the amplitude grows and reaches a maximum, marked with a red diamond in Fig. 1.One finds that the main binary property determining the amplitude of this first postmerger GW amplitude maximum is the mass ratio of the binary q, cf.Fig. 3 with 1 − (5.3149 × 10 −1 )q 1 − (2.3420 × 10 −1 )q .
(3) The qualitative behavior is again related to the possible tidal disruption of the binary close to the merger for unequal mass systems.We note that even if the secondary star does not get disrupted, the maximum density in the remnant shows one peak rather than two independent cores [68,69] which leads to a smaller first postmerger peak and overall on average a smaller GW amplitude.

B. Frequency domain
We obtain the frequency domain waveform by fast Fourier transformation (FFT) of the time domain GW strain h(t): As before, we consider only the dominant 22-mode of the GW signal.
In Fig. 4 we show the frequency domain GW signal (blue solid line) of THC:0001.The merger frequency of this particular configuration is 1638 Hz and it is marked with a red circle.The main feature of the postmerger spectrum is the dominant peak characterizing the main emission frequency f 2 , which for the setup shown in the Figs. 1 and 4 is about 2354 Hz.For a better interpretation, we also present the frequency domain postmerger spectrum in green.Such postmerger-only waveforms are obtained by FFT after applying a Tukey window [70] with a shape parameter 0.05 at t min (where the shape parameter represents the fraction of the window inside the cosine tapered region) and will be used for our injections to test our parameter estimation infrastructure.
f 2 -frequency: The dominant feature in the postmerger frequency spectrum is the dominant emission mode of the merger remnant at a frequency f 2 .As mentioned before, a number of works, e.g., [18][19][20][21][22] have discussed possible EOS-insensitive, quasi-universal relation for the f 2 -frequency.
Building mostly on the work of Bernuzzi et al. [22], we derive a new relation for the f 2 frequency.First, we extend the dataset of 99 NR simulations employed in Ref. [22] and use a set of 121 data by incorporating additional setups published as a part of the CoRe database [59]; cf.Tab.I. Second, we are switching from to which yields a tiny improvement (∼ 0.1%) in the RMS error against the NR data, but more notably relates directly to the mass-weighted tidal deformability Λ measured most accurately from the inspiral part of the signal.In addition to the dependence of the mass-weighted tidal deformability Λ, the postmerger evolution depends also on the stability of the formed remnant and how close it is to the black hole formation.This information is in part encoded in the ratio between the total mass M and the maximum allowed mass of a single non-rotating NS M TOV .By incorporating an additional M/M TOV dependence we are able to reduce the RMS error by ≈ 28%.Therefore, we define a parameter ζ by a linear combination of κ T eff and M MTOV (see also [71][72][73] for a similar approach), The free parameter a = −131.7010 is determined by minimizing the RMS error.Finally, the dimensionless frequency M f 2 is fitted against ζ using a Padé approximant: with α = 3.4285 × 10 −2 , A = 2.0796 × 10 −3 and B = 3.9588 × 10 −3 .We present Eq. ( 8) together with our NR dataset and a one sigma uncertainty region (shaded area) in Fig. 5.
Frequency domain amplitude of f 2 : Finally, we want to briefly discuss the dependence of the f 2 -peak amplitude on the binary properties.While the f 2 -frequency correlates clearly to ζ, we have not been able to find a similar tight relation between any combination of the binary parameters and the amplitude | h22 (f 2 )|.The only noticeably imprint which we have been able to extract comes from the mass ratio q, where generally higher mass ratios lead to a smaller amplitude |r h22 (f 2 )/M | as shown in Fig. 6 with We note that because of the large uncertainty, we see Eq. ( 9) more as a qualitative rather than a quantitative statement about the postmerger spectrum.However, the overall amplitude decreases for an increasing mass ratio seems to be a robust feature and might help to interpret future GW observations.

III. MODEL FUNCTIONS, f2 MEASUREMENT, AND INSPIRAL-POSTMERGER CONSISTENCY A. Lorentzian Approximants
Based on our previous discussion and the dominance of the characteristic f 2 frequency, we start our consideration with a simple damped sinusoidal time domain waveform to model the postmerger waveform.The Fourier transform of a damped sinusoidal function is a Lorentzian function, Eq. ( 10).In the simplest case which we consider, we use 3 unknown coefficients (c 0 , c 1 , c 2 ) corresponding to the amplitude, the dominant emission frequency and the inverse of the damping time, respectively, and write the frequency-domain signal as: Equation (10) suggests that the amplitude peak of the GW postmerger spectrum and also the main postmerger phase evolution are connected to the same frequency characterized by c 1 .
Maximizing over (c 0 , c 1 , c 2 ), we compute the mismatches between the used NR data from the CoRedatabase (Tab.I) and the model function, Eq. (10). Figure 7 shows all mismatches, which on average are ∼ 0.18.
The mismatches can be further decreased by adding three additional coefficients: For Eq. ( 11) the amplitude and phase evolution are independent from each other and we obtain average mismatches of 0.15, i.e., about 17% better than for the threeparameter model.While one might argue that the additional introduced degrees of freedom hinder the extraction of individual parameters in a full Bayesian analysis, it might also be possible that the more flexible 6-parameter model recovers signals with smaller SNRs.Thus, we continue our study with both model functions Eqs. ( 10) and ( 11).Finally, one obtains the plus and cross polarizations from Eq. ( 10) and Eq. ( 11) by incorporating the inclination (ι) dependence: hc , hp can be employed directly to infer information from the postmerger part of a GW signal or to construct a full IMP-waveform for BNSs.

B. Validating the Parameter Estimation Pipeline
In this section, we present for four selected cases the performance of the three-and six-parameter models.We inject the NR waveforms immersed in the same simulated Gaussian noise with total network SNRs ranging from SNR 0 to SNR 10 assuming that Advanced LIGO and Advanced Virgo detectors run at design sensitivity [8] 1 .Fig. 8 shows the injection of THC:0021 with SNR 8 in both the time (top panel) and frequency domain (bottom panel).For each injected waveform, a Tukey window with shape parameter 0.05 is applied at t min to isolate the postmerger signal and avoid Gibbs phenomenon.All simulated signals are injected with zero inclination angle, zero polarization angle ψ and sky location (α, δ) to be (0, 0).We estimate parameters using Bayesian inference with the LALInference module [53] available in the LALSuite package.Sampling is done on 9 (12) parameters 1 The corresponding power spectral density (PSD) files, LIGO-P1200087-v18-aLIGO DESIGN.txt and LIGO-P1200087-v18-AdV DESIGN.txt are available under the LALSimulation module in the LALSuite package.with nested sampling algorithm lalinferencenest [74], where i runs from 0 to 2 (5) for the three-(six-) parameter model, and t c and φ c are the reference time and phase, respectively.The priors are chosen to be uniform in [0, 10 −20 ]s −1 on c 0 , uniform in [1500, 4096]Hz on c 1 and c 5 , uniform in [1,400]Hz on c 2 and c 4 , uniform in [0, 6] on c 3 , uniform in [0, 2π] on α, ψ and φ c , uniform in [−1, 1] on cos(ι) and sin(δ), and uniform in [trigger time − 0.05s, trigger time + 0.05s] on t c where the trigger time is the signal arrival time at the geocentric frame.
Fig. 9 shows the posterior for c 1 , i.e., our best estimate of the f 2 -frequency for SNRs up to 8 for our 4 examples, which we mark in Tab.I. We present the recovery with the 3-and the 6-parameter model in the top and bottom panels, respectively.The solid vertical line represents the injected f 2 -frequency and the dashed line represents the FIG. 9. Posteriors for the parameter c1 in Eq. ( 10) (top panels) and Eq. ( 11) (bottom panels) for a variety of SNRs.c1 can be directly related to the peak in the frequency domain spectrum and therefore relates to the f2 frequency.The IDs of the four injected waveforms are # 18, 20, 32, 34, i.e., THC:0021, THC:0031, BAM:0048, BAM:0057 of [59].The chosen set covers various EOSs, mass ratios, and masses and is therefore used as a testbed for our new algorithm.FIG. 10.Posteriors for the parameter c1 and c5 in Eq. ( 11) for 2 SNRs.c5 peaks at a frequency close to f2 but it is significantly less constrained than c1.
estimate according to the quasi-universal relation Eq. ( 8) together with a one-sigma uncertainty (gray shaded region).
We summarize the main findings below: (i) The three-parameter and six-parameter approxi-mants perform similarly.
(ii) Depending on the exact setting (e.g., intrinsic source properties, noise realization, sky location) one can recover the f 2 frequency with an SNR of ∼ 4 for the best and ∼ 8 for worst considered scenarios.
(iii) Interestingly, one finds that also c 5 relates to a frequency which is close to the f 2 frequency, however, c 5 is significantly less constrained than c 1 (Fig. 10).
(iv) Once 3rd generation detectors are available and so the postmerger SNRs of ∼ 10 are obtained, the systematic uncertainties of the quasi-universal relations become larger than the statistical uncertainties; cf.dashed and solid, vertical black lines.

C. Inspiral and Postmerger consistency
Finally, we want to illustrate how a future detection of a postmerger GW signal will help to constrain the source properties and the internal composition of NSs.As shown before, the f 2 -frequency can be extracted through a simple waveform model (or alternatively by using BayesWave, e.g., [36]).To connect the f 2 -frequency with the source parameters, one needs to employ quasi-universal relations as presented in Sec.II and some information obtained from the analysis of the inspiral GW signal.In particular, the total mass M can be measured precisely using state-of-art BNS inspiral waveforms, e.g., [10,14,15 77].For GW170817 the uncertainty of M could be reduced to ±0.04M once EM information had been included [71,78].Thus, we will use an uncertainty of ∆M = ±0.04Mas a conservative estimate.addition, we have to know the maximum TOV-mass M TOV .
Current estimates for M TOV are based on the observation of J0740+6620 [79] with M = 2.17 +0.11 −0.10 M and the assumption that GW170817's endstate was a black hole [80][81][82][83] such that M TOV 2.17-2.35M .Due to the increasing number of BNS detections in the future, we expect that the uncertainty of M TOV can be considerably reduced so that we will use an uncertainty of ±0.04M .From this information, we can compute the ζinterval consistent with the observed inspiral signal from the Λ posteriors, cf.vertical red shaded region in Fig. 11.We then connect the ζ estimate obtained from the inspiral with Eq. ( 8) and the f 2 -measurement of the postmerger signal.This consistency analysis is somewhat connected to the inspiral-merger-ringdown consistency test for BBH [84], but not only has to assume the correctness of general relativity but also that our understanding of supranuclear matter and the EOS-insensitive quasiuniversal relations are valid.
Figure 11 summarizes our main results.Generally the GW measurements can be considered consistent between the inspiral and postmerger observations as well as with the quasi-universal relation relating the main postmerger frequency with the binary properties.We find that for all cases, the quasi-universal relation and its 1-sigma uncertainty region lie within the intersection of the red shaded region (inspiral) and the blue/orange/green horizontal regions (postmerger).Thus, (as expected) all simulations are consistent with (i) general relativity and (ii) the nuclear physics descriptions used as a basis for NR simulations to derive the quasi-universal relations.For future events this approach will allow us to probe our understanding of physical processes under extreme conditions, and in cases where the quasi-universal relation seems violated to even derive new relations based on GW measurements (under the assumption that general relativity is correct).We note that even in the case where either general relativity or the quasi-universal relations would be violated, we might not be able to determine reliably the violation based on one individual event, but stronger constraints can be obtained by combining multiple BNS events.

IV. SUMMARY
In this work, we discussed the general morphology of a BNS postmerger in both the time and frequency domain.We presented quasi-universal relations for the time at which the first postmerger amplitude minimum happens and the strength of the first postmerger amplitude maximum.In general, the time between the merger and the amplitude minimum increases with an increasing Λ, while the amplitude of the first postmerger maximum decreases with an increasing mass ratio.In the frequency domain, we improved the existing quasi-universal relations of M f 2 by extending the employed NR dataset (121 simulations in total) and adding an extra dependence of M/M TOV .The extra term M/M TOV characterizes how close the setup is to the black hole formation.
We find that a three-(six-) parameter Lorentzian can model the postmerger waveform with average mismatch of 0.18 (0.15).To test these model functions, we performed an injection study, in which we simulated the detector strains with four different BNS configurations immersed in the same simulated Gaussian noise assuming Advanced LIGO and Advanced Virgo at design sensitivities.We find that in the best cases the Lorentzian models could measure the dominant emission frequency f 2 once the signal has an SNR of 4 or above; however, for most scenarios higher SNRs ∼ 8 were required.
Employing the new quasi-universal relation for M f 2 described in this work, we could present consistency tests between the inspiral and the postmerger signal; cf.Fig. 11.

FIG. 3 .
FIG.3.A scatter plot of the first peak in the postmerger spectra after merger versus q.The color shows the massweighted tidal deformability Λ.The shaded region indicates the 1-sigma (±1.7915 × 10 −2 ) uncertainty.

FIG. 5 .
FIG. 5. M f2 as a function of ζ.The color shows the massweighted tidal deformability Λ.The shaded region indicates the 1-sigma (±1.1025 × 10 −3) uncertainty.In addition to the CoRe-dataset employed to derive the previously shown quasi-universal relations, we include here the published results of[31,33].

FIG. 7 .
FIG. 7. Mismatch between a subset of the NR data listed in Tab.I and the three-and six-parameter models.
FIG. 11.Schematic plot showing how one can constrain ζ from f2 measurements where the spread f2 are measured by one standard deviation.Top panel is for three-parameter model and the bottom panel is for six-parameter model.The vertical red shaded region corresponds to the ζ-interval consistent with a hypothetical inspiral signal assuming an uncertainty of ±0.04M for M and MTOV, and ±30 for κ T eff .The exact value is marked as vertical red-dashed line.