Probing neutron stars with the full pre-merger and post-merger gravitational wave signal from binary coalescences

The gravitational wave signal emitted during the coalescence of two neutron stars carries information about the stars' internal structure. During the long inspiral phase the main matter observable is the tidal interaction between the binary components, an effect that can be parametrically modeled with compact-binary solutions to General Relativity. After the binary merger the main observable is frequency modes of the remnant, most commonly giving rise to a short-duration signal accessible only through numerical simulations. The complicated morphology and the decreasing detector sensitivity in the relevant frequencies currently hinder detection of the post-merger signal and motivate separate analyses for the pre-merger and post-merger data. However, planned and ongoing detector improvements could soon put the post-merger signal within reach. In this study we target the whole pre-merger and post-merger signal without an artificial separation at the binary merger. We construct a hybrid analysis that models the inspiral with templates based on analytical calculations and calibrated to numerical relativity and the post-merger signal with a flexible morphology-independent analysis. Applying this analysis to GW170817 we find, as expected, that the post-merger signal remains undetected. We further study simulated signals and find that we can reconstruct the full signal and simultaneously estimate both the pre-merger tidal deformation and the post-merger signal frequency content. Our analysis allows us to study neutron star physics using all the data available and directly test the pre-merger and post-merger signal for consistency thus probing effects such as the onset of the hadron-quark phase transition.

The gravitational wave signal emitted during the coalescence of two neutron stars carries information about the stars' internal structure. During the long inspiral phase the main matter observable is the tidal interaction between the binary components, an effect that can be parametrically modeled with compact-binary solutions to general relativity. After the binary merger the main observable is frequency modes of the remnant, most commonly giving rise to a short-duration signal accessible only through numerical simulations. The complicated morphology and the decreasing detector sensitivity in the relevant frequencies currently hinder detection of the postmerger signal and motivate separate analyses for the premerger and postmerger data. However, planned and ongoing detector improvements could soon put the postmerger signal within reach. In this study we target the whole premerger and postmerger signal without an artificial separation at the binary merger. We construct a hybrid analysis that models the inspiral with templates based on analytical calculations and calibrated to numerical relativity and the postmerger signal with a flexible morphology-independent analysis. Applying this analysis to GW170817 we find, as expected, that the postmerger signal remains undetected. We further study simulated signals and find that we can reconstruct the full signal and simultaneously estimate both the premerger tidal deformation and the postmerger signal frequency content. Our analysis allows us to study neutron star physics using all the data available and directly test the premerger and postmerger signal for consistency thus probing effects such as the onset of the hadron-quark phase transition.
The GW signal from a binary NS (BNS) coalescence consists of two parts: a premerger and a postmerger. In the premerger phase the binary components inspiral toward each other as they lose orbital energy to GWs and eventually collide [29]. In the case of GW170817, the detectable premerger signal increased in frequency from ∼ 23Hz to many hundreds of Hz as the orbital separation between the two NSs rapidly decayed over ∼ 2 minutes [30] before finally merging at a frequency of ∼ 1500Hz [31]. During the late stages of the inspiral, tidal interactions accelerate the binary evolution and leave an imprint on the signal that depends on the NS matter properties and can thus be used to infer the EoS [32]. The premerger GW signal is typically modeled with compact binary coalescence (CBC) templates that are based on approximate solutions to the general relativity field equations and numerical relativity (NR) simulations, see [33] for a recent review.
After the inevitable merger, the remnant star evolves in a way that depends on the system mass and the EoS, see [34][35][36][37] for reviews. For most EoSs and NS masses, a hypermassive NS supported by differential rotation and thermal effects is expected to be formed and sustained for ∼ 10 − 100ms [38]. During this time, the NS remnant emits a short-duration postmerger GW signal with a characteristic peak at 1500 − 4000Hz [39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57] depending on the mass and EoS, a promising frequency range for upcoming GW detectors 1 . In what follows this is the case we consider. While the exact physics governing the postmerger phase is not fully understood and NR simulations are complicated by factors such as thermal effects, turbulence, magnetohydrodynamical instabilites, neutrinos, and possible phase transitions, simulations make a robust empirical prediction for the frequency of the peak of the postmerger spectrum given the NS mass and EoS.
Information extracted from a postmerger signal about NS properties, is complementary to the premerger. Analysis of a postmerger signal would probe different density [58] and temperature [59,60] regimes of the EoS, as the more massive, hot remnant will have a higher maximum density than the premerger NSs [40]. This complementarity can lead to insights about high-density phenomena such as phase-transitions between the densities probed by the premerger and postmerger signals [61][62][63][64][65][66][67][68][69]. Finally, detection of a postmerger signal would allow us to determine the nature of the merger remnant, namely NS or BH, which has implications for a potential electromagnetic counterpart and its interpretation [70] as well as the determination of the threshold mass for prompt collapse [71].
In the case of GW1708017 and GW190425 only piecewise analyses of the GW signal have been performed focusing either on the premerger or on the postmerger part of the signal [6,30,72,73]. This reduces the computational cost significantly, as the postmerger signal requires a large sampling rate (typically 8194Hz) while the premerger analysis involves long duration data segments (typically 128s). As no postmerger signal is detectable yet [30,72], excluding it from the premerger analysis did not bias the results for the binary parameters [74]. However, looking forward and in the case of a postmerger signal detection, separate analyses will not be able to exploit phase coherence through merger or parameter relations, such as those between the tidal deformability from the premerger signal and the frequency peak of the postmerger spectrum [52,58,75].
Another reason for separate premerger and postmerger analyses is that the morphology and details of the postmerger signal are not well understood, making it difficult to construct a first-principles physical model such as the premerger CBC templates. While approximate analytic models for the postmerger signal have been proposed with the aid of NR simulations [76][77][78][79][80][81][82], these are still limited in accuracy as well as by uncertainties of the simulations they are based on [80]. An alternative are morphology-independent analyses that are not limited to a specific signal type [83][84][85]. Specifically, a modelagnostic approach with the BayesWave algorithm [86][87][88] has been shown to accurately determine the main features of the postmerger signal [31,85], and in some cases do so more accurately than tailored models [81].
BayesWave models a GW signal with a sum of sine-Gaussian wavelets. Both the number of wavelets and their parameters are marginalized over, resulting in a flexible analysis that has been applied to a variety of signals [89][90][91][92][93][94][95][96][97][98]. BayesWave's sine-Gaussian wavelets are particularly adept for postmerger signals that are dom-inated by distinct frequency components, as each component can be approximately modeled by a wavelet. In this context, BayesWave has been used to search for a short-duration, high-frequency signal after both GW170817 [30] and GW190425 [6], returning null results and upper limits on the energy content [30]. A detailed study of BayesWave's sensitivity to postmerger signals is presented in [85].
In this paper we construct a hybrid analysis of the full BNS GW signal that targets both the premerger and the postmerger data. We simultaneously analyze the full signal using 1) a CBC template to describe the well-modeled premerger phase and 2) sine-Gaussian wavelets to capture the less well understood postmerger. The parameters of the CBC template and the wavelets are simultaneously sampled over to obtain the combined multidimensional posterior for all components of the model. Those parameters include the premerger tidal deformability that quantifies the inspiral tidal deformation and the peak frequency of the postmerger spectrum. Our analysis extracts both simultaneously and allows for direct consistency comparisons.
After validating the hybrid modeled/unmodeled analysis on a toy model based on GW150914 data, we apply it to GW170817 and find that the postmerger signal remains undetected, as expected from the upper limit estimates of [30]. We further analyze simulated signals of high signal-to-noise ratio (SNR) for which the postmerger signal is detectable. We show that our analysis can reconstruct the premerger signal as well as the main components of the postmerger signal such as the dominant frequency mode. We simulate signals for which the premerger and postmerger parts are consistent (corresponding to a hadronic EoS) and inconsistent (corresponding to an EoS with a strong phase transition in the relevant density scales) and demonstrate that our analysis can extract either behavior.
The rest of the paper is organized as follows. In Sec. II we describe in detail our analysis and algorithm. In Sec. III we present a proof-of-principle analysis on GW150914. In Sec. IV we reanalyze the GW170817 data. In Sec. V we apply our analysis to simulated data. In Sec. VI we conclude.

II. METHODOLOGY AND MODELS
In order to analyze a GW signal that contains both modeled and unmodeled features we employ a hybrid analysis that makes use of both CBC templates and flexible models for the signal. We base our analysis on the morphology-agnostic data analysis algorithm BayesWave [86,87] by extending it to account for CBC templates similar to [99]. In its core functionality BayesWave simultaneously models a GW signal (referred to as the "signal model" in BayesWave literature), instrumental glitches (the "glitch model"), and the gaussian detector noise (the "PSD model"). BayesWave uses minimal assumptions, i.e., no physical model, to describe these components. Instead, the signal and the glitch are modeled as the sum of a variable number of sine-Gaussian wavelets, whereas the noise model describes the power spectral density (PSD) of the Gaussian noise with a variable number of spline points and Lorentzians. For this work and [99] a fourth component has been added, which models GW signals through CBC templates obtained by solving the two-body problem in general relativity; the CBC model.

A. Overview of models
BayesWave uses sampling methods to characterize the multi-and variable-dimension posterior distribution for independent models that target different features of the GW data. The goal is to draw samples from the posterior distribution function p(M |d), the probability that a model M describes the data d and a prior p(M ), defined as where p(d) is the evidence, and p(d|M ) expresses the likelihood of the data for a given M . We model the data in each interferometric detector d I as a linear combination of multiple components where g I denotes any detector glitches 2 , the full GW signal (h cbc I + h w I ) consists of the CBC waveform model h cbc I and any signal (such as the postmerger) that is not included in the CBC template and is captured by wavelets h w I , and n I is the Gaussian component of the noise. The combination of a signal, glitch, and noise defines the full model M . Given such a model, the likelihood function p(d I |M I ), in short L I , can be computed by noting that the residual r I = d I − M I is consistent with Gaussian noise (see also Eq. (3) of [87]): where C I is a constant that depends on the PSD of the detector noise S n,I (f ), and the noise weighted inner product is defined as where an asterisk denotes a complex conjugate. The remaining ingredients of the analysis are a specific model for each data component, the corresponding prior, and the sampling procedure. We express the GW signal at geocenter and then project it onto each detector in the network to compute the detector response in the frequency domain The intrinsic part of the signal for each GW polarization mode is given by h cbc/w + and h cbc/w × and it is different under the CBC or the wavelet model. The projection then involves a time delay in the arrival of the signal in the different detectors ∆t(α, δ), the detector antenna patterns F + (α, δ, ψ), F × (α, δ, ψ) for each GW polarization mode, parametrized through the sky location α, δ corresponding to the right ascension and declination respectively, and the polarization angle ψ. The projection is independent of whether the signal is expressed through CBC templates or wavelets so the functional form of Eq. (5) and the parameters (α, δ, ψ) are common. The ratio h cbc/w + /h cbc/w × is related to the binary inclination ι and is also shared between the CBC and wavelet models.
By construction BayesWave breaks down the full parameter space defined by all models in different blocks and uses a blocked Gibbs sampler to explore the posterior. Each block of parameters is sampled over with different (reversible jump [100]) Markov Chain Monte Carlo (RJMCMC) samplers. Parallel tempering is employed both for efficient posterior exploration and evidence calculation [101], so all samplers share the same number and temperature of chains and exchanges between the chains proceed for all models simultaneously. In order for this procedure to lead to efficient exploration of the posterior distribution, each block of parameters needs to be selected to be as independent from and uncorrelated to other blocks as possible.
The sampler ensemble consists of (1) a CBC MCMC, Sec. II B, (2) an RJMCMC sampler for the wavelet intrinsic parameters, Sec. II C, (3) an MCMC sampler for the common extrinsic parameters for the CBC and wavelet models, Sec. II D, and (4) an RJMCMC noise model sampler, Sec. II E. The structure is given in Fig. 1 and it is similar to the code used in [99] but instead of nonoverlapping parameter blocks, we now allow for a set of common parameters to be updated in multiple blocks as detailed in the subsequent subsections. Sampling then proceeds iteratively both within each block (O(10 2 ) iterations at a time) and between different blocks (O(10 4 ) iterations). In the rest of this section, we discuss each (RJ)MCMC, which parameters it updates, and how it relates to the other samplers. The discussion is structured to be as self-contained and detailed as possible; readers primarily interested on the results can find them starting in Sec. III.

Gibbs sampler
General code workflow. After some preconditioning, the data are used to obtain a quick fit to the noise PSD Sn(f ) and CBC θ cbc and extrinsic θext parameters. These are then used as a starting point for the blocked Gibbs sampler (red box) which consists of 4 independent samplers. Each iteration of the Gibbs sampler consists of a wavelet RJMCMC (blue box) that updates the wavelet parameters θ w , an extrinsic MCMC (green box) that updates the extrinsic parameters θext, a noise RJMCMC (gray box) that updates Sn(f ), and a CBC MCMC (pink box) that updates θ cbc . Orange rectangles enclose the parameter that is being updated from the preceding box. Each block consists of O(10 2 ) iterations, while the entire red Gibbs box consists of O(10 4 ) iterations.

B. CBC-specific parameters: CBC sampler
A nonprecessing quasicircular inspiral of two compact objects is characterized by a set of 6 intrinsic source parameters and a set of 7 extrinsic parameters that describe the binary's location and orientation relative to the detectors. The CBC sampler updates the parameters that are specific to h cbc + and h cbc × , namely masses, spins, tidal parameters (for BNSs), distance, time, and phase of coalescence. In Fig. 1 the CBC sampler is depicted by the pink box "CBC MCMC" and the relevant parameters are collectively denoted by θ cbc .

Waveform model and Parameters
The CBC sampler is integrated with the LALSimulation [102] waveforms suite and can make use of any non-precessing waveform model; we here focus on PhenomD [103,104] and PhenomDNRT [105][106][107]. The former includes the full inspiral-merger-ringdown (IMR) signal expected from a BBH coalescence. The latter targets the inspiral of BNS systems as it includes the effect of tidal interactions between the NSs, manifesting as a change in the waveform phase. PhenomDNRT models the long inspiral phase and instead of terminating abruptly at merger, it decays rapidly through a window function after the merger frequency [106,107]. PhenomDNRT is well suited for our purposes regarding BNS signals as the waveform remains smooth and continuous, naturally allowing the resulting CBC (inspiral) + wavelet (postmerger) model to be smooth without any forced stitching.
• The masses of the binary components are m 1 and m 2 with the convention m 1 ≥ m 2 .
• The dimensionless spin components aligned with the orbital angular momentum are χ 1 and χ 2 . Current detector sensitivities make individual spin components hard to measure [108,109] so a commonly used spin combination is the mass-weighted effective spin parameter, which is conserved under spin-precession to at least the second post-Newtonian order [110].
• Tidal interactions in BNS signals are encoded by the NS dimensionless tidal deformabilities where k 1,2 is the tidal Love number, R 1,2 is the radius, and m 1,2 is the mass of each binary component. To leading order, the tidal effects are imprinted in the waveform through the binary tidal deformability [111,112] Λ = 16 13 with subleading terms introducing another tidal parameter that is generally not measurable, δΛ [112]. For BBH systems, Λ 1 = Λ 2 = 0 [113,114].
• The CBC sampler further handles 3 extrinsic parameters either because they do not explicitly appear in the wavelet model (time of coalescence t c , luminosity distance D L ) or because they are correlated with intrinsic parameters and need to be simultaneously sampled for efficiency (phase of coalescence φ c ). Despite the well-known correlation between the distance and the inclination [115], we do not update the latter here as the inclination affects the wavelet model as well.

Priors
The priors on the binary intrinsic parameters depend on whether we consider BNSs or BBHs. We sample the detector frame chirp mass M ∈  [116]. We use flat priors for the dimensionless spin components χ 1 , χ 2 ∈ [−0.05, 0.05] for BNSs and χ 1 , χ 2 ∈ [−1, 1] for BBHs. In the BNS case we sample the effective dimensionless tidal parameters uniformΛ ∈ [0, 1000] and δΛ ∈ [−500, 500] with the constraint that the corresponding individual component tidal deformabilities cannot be negative. We can optionally directly sample the individual component tidal parameters uniformly with the default prior range Λ 1,2 ∈ [0, 1000] though all results presented here use the former tidal prior. In all cases, the individual component tidal deformabilities are assumed to be independent of each other and not linked under any assumption about the NS EoS [117][118][119].
We sample the merger time relative to the GPS trigger time (t trig ) with a default window of t c ∈ [t trig − 0.5s, t trig + 1.5s]. Alternatively, it is possible to use a symmetric window of arbitrary size around the trigger time. We also sample in the logarithmic luminosity distance with an adjustment to the Jacobian such that the prior is uniform in volume (i.e., uniform in D 3 L ) with D L ∈ [1, 10000] Mpc. Finally, the coalescence phase prior is uniform in φ c ∈ [0, 2π].
Since the CBC and the wavelet models are active simultaneously, the overall likelihood has to take into account both CBCs and wavelets. The blocked nature of the code means that the wavelet model is not updated in the CBC sampler, so the wavelet-subtracted data d cbc ≡ d − h w are constant as the CBC parameters are being updated. We compute d cbc = d − h w once at the beginning of the CBC sampler block given the current wavelet parameters. The wavelet-subtracted data d cbc then play the role of the data in the sense that d cbc − h cbc is the residual during the likelihood computation.

Proposals
We employ the CBC sampler described in detail in [120] within the BayesWave blocked Gibbs sampler, see Fig. 1. The sampler performs by default 300 iterations that update the CBC parameters described above using a mixture of proposals: Fisher jumps: Gaussian jumps along the eigenvectors of the Fisher information matrix. The Fisher information matrix is recomputed at the start of each CBC sampler block for each parallel chain and using each chain's current location. Differential evolution [121] : Two samples from each chain's history propose a new location. The initial history array is drawn from the prior and then updated every 20 iterations. As each CBC sampler block does 300 iterations by default, the history array contains information about previous blocks. Jiggle: Small Gaussian jumps in individual parameter directions. The size of the proposed jump is scaled with the chain temperature such that high temperature chains have larger jumps and can explore a larger area of the parameter space. The jiggle proposal allows for exploring areas outside those covered by Fisher matrix jumps and are especially useful early on. Mixed : Uniform proposal for the tidal parameters and differential evolution (see above) for the other parameters. This proposal is employed 5% of the time and aids convergence for cases where the tidal parameters cannot be tightly constrained. The proposal facilitates taking larger steps for the tidal parameters whereas the other better constrained parameters evolve with smaller steps, thus leading to a higher acceptance rate than a fully uniform proposal.

Heterodyned likelihood
We integrate the heterodyned likelihood technique [120,122] in the BayesWave CBC (and extrinsic, to be described later) sampler in order to speed up the likelihood evaluation. The technique is based on the fact that all CBC waveforms that contribute to the posterior (i.e. have a likelihood high enough that the MCMC proposal can be accepted) are similar. The likelihood of a model h cbc can then be computed based on a fixed and precomputed reference model h cbc ref that is similar to h cbc , and given the noise PSD and the data.
We can rewrite the residuals r using the reference heterodyne, omitting the detector subscript I for clarity, as: where ∆h cbc is the difference between the reference waveform and the exact waveform at a given set of CBC parameters, d cbc is the wavelet-subtracted data as given to the CBC sampler, The likelihood from Eq. (3) then becomes where we omit the second term of Eq. (3) as the PSD is held fixed during the CBC updates and the term does not affect the relative likelihoods. This form of the likelihood enables several computational advantages as we split the inner product in its components A detailed description of how each components of the heterodyned likelihood above is calculated is available in [120,122]; we here briefly discuss how this is adapted and implemented within the BayesWave blocked Gibbs sampler. The first term of the inner product (r ref |r ref ) is the most straightforward as it is computed once at the beginning of each CBC sampler block using the current wavelet-subtracted data and CBC parameters for each chain.
The key performance improvement of the heterodyne likelihood is achieved by recognizing that the latter two terms in Eq. (11) can be further split up into the product of two terms: a slowly varying term that is computed at each iteration, and a rapidly varying term that is computed once. Both terms are evaluated through a sum over frequencies with bin width ∆f using a Legendre polynomial expansion [120]. The rapidly-varying term requires a small ∆f but this is computed once. The slowly-varying term that needs to be calculated on each likelihood evaluation can instead make use of a coarse frequency grid of ∆f ≤ 4 Hz thus significantly reducing computational cost.
The term (∆h cbc |∆h cbc ) in Eq. (11) introduces a slowly varying phase difference as a function of frequency between the reference waveform and the proposed waveform. However, in the noise weighted inner product, Eq. (4) the denominator might not be slowly varying if the PSD includes sharp spectral lines. Following the approach discussed above, we split the PSD in a smooth broadband component that is slowly varying as a function of frequency, and line features which are only evaluated at discrete intervals. If the PSD is held fixed throughout the analysis, the smooth part is computed once at the beginning. When using a variable PSD model instead, we update the smooth PSD component at the beginning of the CBC sampler block.
The final term of Eq. (11) combining r ref and ∆h cbc is split up into the product of a slowly and rapidly varying part by introducing the smooth component of the PSD in the numerator and denominator of, respectively, the rapidly and slowly varying part. Due to this, the heterodyned difference in the waveforms ∆h cbc is effectively whitened using only the smooth PSD component on a coarse frequency grid. The rapidly varying part consisting of the heterodyned whitened residuals r ref now includes an extra product with the smooth part of the PSD, but since this is only computed at the start of the CBC sampler block, it adds little computational cost.
The heterodyned likelihood speeds up the likelihood calculation by a factor of ∼ T obs ∆f , where T obs is the duration of the analyzed data, thus making it especially useful for long duration signals such as BNSs. The heterodyne implementation in BayesWave is more complex than the one in [120] as it has to account for a potentially changing PSD and the wavelet model. We thus expect more modest overall computational improvements, though the speed-up of each individual likelihood calculation is consistent with [120].

C. Wavelet parameters: Wavelet sampler
BayesWave has a common RJMCMC sampler for its signal and glitch models where the non-Gaussian features (either a GW signal or a detector glitch) are modeled as a sum of a variable number of sine-Gaussian wavelets which can be expressed in the time-domain as where τ ≡ Q/(2πf 0 ) is the decay time, Q is the quality factor, A is the wavelet amplitude, f 0 is the wavelet central frequency, t 0 is the wavelet central time and φ 0 is the phase offset. The glitch model uses independent sums of wavelets for the different detectors, whereas the signal models use a set of wavelets at geocenter and projects them onto detectors through Eq. (5) with where is the ellipticity parameter and N s is the number of wavelets in the signal model. Details about the wavelet sampler and priors/proposals for each wavelet parameter and the number of wavelets are discussed in [86,87] and remain mostly unaltered in our analysis. The prior on the number of wavelets has an upper limit of 100.
We extend the standard settings by adding an option to limit the central frequency range of the wavelets to f 0 ∈ [f w min , f w max ] which defaults to the analysis bandwidth. Since the ellipticity is related to the remaining extrinsic parameters and to the binary inclination, it is not sampled by the wavelet sampler but the extrinsic sampler discussed next.
In Fig. 1 the wavelet sampler is depicted by the blue box termed "wavelet RJMCMC" where θ w corresponds to the relevant parameters: the number of wavelets and the amplitude, quality factor frequency, time, and phase of each wavelet. Similarly to the CBC sampler, the wavelet sampler makes use of the CBC-subtracted data d w ≡ d − h cbc , which are computed once at the beginning of each wavelet sampler block using the current CBC sample.

D. Common parameters: Extrinsic sampler
Since both the CBC and wavelet models target the same GW source, they share a number of extrinsic source parameters that govern the signal projection onto the detector network, Eq. (5). These common parameters θ ext are updated by the extrinsic sampler block, denoted as "extrinsic MCMC" in Fig. 1.
• The source sky location is given by the right ascension α ∈ [0, 2π] and the declination sin δ ∈ [0, 1], and the corresponding prior is flat on the sphere.
• The degree of elliptical polarization in the signal is encoded in ∈ [−1, 1], Eq. (14). For elliptically polarized signals, the ellipticity parameter can be related to the orbital inclination with respect to the line of sight to the binary ι as We use a uniform prior in cos ι ∈ [−1, 1], which is a deviation from previous works where BayesWave samples uniformly in instead.
• An overall phase shift ϕ 0 is sampled uniformly as ϕ 0 ∈ [0, 2π] and it is equivalent to a phase shift in all wavelets simultaneously (see Eq. 12). Despite the degeneracy, introducing this parameter leads to more efficient convergence in the wavelet model.
• We also introduce an overall amplitude scaling parameter A cbc , which is sampled uniformly and is degenerate with the CBC luminosity distance D L . We restrict A cbc to correspond to values for D L within the D L prior range and add a Jacobian factor such that the prior remains uniform in volume, i.e., (D L /A cbc ) 3 .
The extrinsic parameters can be efficiently sampled independently from the CBC intrinsic parameters by noting that a change in the extrinsic parameters induces an amplitude, time, and phase shift in the signal seen in each detector, while leaving its phase evolution unaltered. We therefore compute the geocenter CBC waveform given the fixed intrinsic parameters once and then apply the shifts accordingly to evaluate the waveform at each new proposed set of extrinsic parameters. This has the advantage that we can sample the extrinsic parameters efficiently while not having to recompute the full, expensive CBC waveform phase at each iteration.
In practice we write the CBC response at a given detector I using the geocenter waveform h cbc + and applying the shifts as where ∆t I is the arrival time at the detector relative to the geocenter, the phase shift is and F cbc I is the magnitude of the detector response Even though A cbc is completely degenerate with D L , we find that it increases sampling efficiency in the extrinsic sampler due to the distance-inclination degeneracy [115]. At the end of the block of extrinsic updates, the scaling and shift parameters are reset and propagated to the corresponding CBC parameters before proceeding to the CBC sampler In these equations ∆t a and ∆ϕ 0 correspond to the overall time and phase shift induced on the waveform by the change in extrinsic parameters. The time of arrival of the merger at geocenter t c changes due to the new sky location and this shift is encoded in ∆t a .

Proposals
To update the parameters θ ext the extrinsic sampler uses a mixture of proposals, which are described in detail in [86,87]. In the following we focus on some CBCspecific updates made to the existing proposals.
The most commonly used proposal is the Fisher proposal, where jumps are proposed along the eigenvectors of the Fisher information matrix. The matrix elements, eigenvalues, and eigenvectors are computed once at the start of the extrinsic sampler block by finite differencing using the current chain position. We extend the existing BayesWave Fisher calculation by adding the CBC model such that the final Fisher matrix captures information from the full CBC+wavelet signal. When the sky location is known and fixed (for example for sources with an electromagnetic counterpart), the Fisher matrix is reduced to the parameters that are varied.
The BayesWave extrinsic sampler also employs a twopart proposal that exploits a degeneracy of the extrinsic parameters when multiple detectors are used (also proposed by [123]). The first part consists of a sky-ring proposal that finds a new sky location (α, sin δ) such that the time delays between the detectors are preserved [86,124]. The second part of the proposal updates the remaining parameters (ψ, cos ι, A CBC and ϕ 0 ) either uniformly in their prior range or deterministically such that the waveform is identical at the new sky location compared to the waveform at the previous sky location, see [120] for details. Since the waveform is conserved, so is the likelihood, though this does not guarantee that the proposed sky location will be accepted due the necessary Jacobian factor [120].
Finally, the BayesWave extrinsic sampler includes a uniform proposal that draws from the prior of each parameter. These jumps are particularly useful for higher temperature chains to explore the entire space available as well as early on in the sampling if a signal has not been found.

Heterodyned likelihood
A change in the extrinsic parameters θ ext will affect both the CBC and the wavelet models, therefore both need to be taken into account when computing the likelihood. However, we find that the CBC part of the model dominates the computational cost so we again use the heterodyned likelihood to speed up the calculation. In this case, there are some differences compared to the implementation for the CBC sampler as discussed in Sec. II B 5 because we can no longer assume that the wavelet-subtracted data are fixed. Indeed, a new set of extrinsic parameters changes the full signal projection and the full residual r = d − h w − h cbc has to be recomputed on each iteration.
Despite this, we can achieve a computational improvement by splitting up the heterodyne computation even further to avoid recomputing all components of the heterodyne when the residuals change, i.e., at every iteration. All subcomponents of the heterodyne that do not involve the residuals are computed and stored once at the start of each extrinsic MCMC sampler block, whereas the remaining components are computed at each iteration. If the signal wavelet model is not employed (i.e. an analysis with only the CBC model), the heterodyne likelihood implementation falls back to the same, more efficient, version as for the CBC sampler.
Since the likelihood calculation is more involved and parts of the heterodyne procedure need to be recomputed in each iteration, we expect more modest computational improvements compared to those of Sec. II B 5.

E. Noise parameters: PSD sampler with BayesLine
The PSD of the Gaussian noise is modeled as a combination of splines and Lorenzians targeting the broadband behavior and spectral lines respectively and marginalized over through an RJMCMC in BayesLine [87,125,126]. The analysis presented here makes no updates on this noise marginalization procedure and can optionally include BayesLine as one of the blocks, corresponding to the "Noise RJMCMC" gray box in Fig. 1. Alternatively, the analysis can skip the noise marginalization and use a predetermined fixed PSD. If the PSD is modeled by BayesLine, the heterodyne procedure needs to be repeated at the end of the BayesLine sampling block and before moving to CBC sampling with the latest PSD sample. However, this coincides with a recomputing of the heterodyne elements due to a change of the wavelet model that is subtracted from the data so no further computations are necessary. In the current implementation [125], the wavelet model is subtracted from the data before entering the BayesLine sampling block and we extend this subtraction to include the CBC model.

III. TOY MODEL: GW150914
To illustrate the wavelet+CBC model, we begin with a toy model by analyzing the first detected BBH merger, GW150914 [127]. In reality, the full BBH signal can be modeled with existing CBC templates, something confirmed in the case of GW150914 through comparisons between different models and physical effects [128]. In our case, and in order to test our hybrid analysis on what is probably the best studied signal, we pretend we lacked efficient modeling of a BBH merger and waveform templates terminated shortly before merger.
We analyze 4s of data from LIGO Hanford and LIGO Livingston available through the Gravitational Wave Open Science Center (GWOSC) [129,130] in the frequency band from 16Hz to 2048Hz. The noise PSD is marginalized over as discussed in Sec. II E. For the CBC model, we use either the PhenomD waveform or the PhenomDNRT waveform. PhenomD is appropriate for BBH systems as it models the full IMR signal. On the other hand, we do not expect PhenomDNRT to accurately model the BBH data as it lacks an accurate merger and ringdown portion due to tapering [106]. Regardless, we use it to demonstrate how any leftover signal that is not captured by the CBC waveform can be captured by the wavelets.
We perform three analyses with different models and CBC templates: 1. A full "CBC IMR " analysis where the data are described only by the CBC model with PhenomD.    Fig. 2, the GPS time (ttrig), the models active, the CBC waveform, the wavelet window around the trigger time, and the wavelet bandwidth.

A "CBC
Insp +wavelets" analysis where the data are described by a combination of the CBC model with PhenomDNRT and the wavelet model. We restrict the wavelet model to times t 0 ∈ [t trig − 0.025s, t trig + 0.025s] and frequencies above 150Hz, targeting the merger portion of the signal that is missed by PhenomDNRT.
The run labels and relevant settings are given in Table I. In Fig. 2 we show the signal reconstructions and parameter posteriors for these three analyses. The left panel shows the whitened data and the 90% credible intervals for each reconstruction in LIGO Hanford (top) and LIGO Livingston (bottom). All analyses are able to identify the GW signal and they further lead to consistent signal reconstructions of the early portion of the signal. However, as the signal approaches the merger phase, the CBC Insp analysis deviates from the reference CBC IMR analysis as well as the hybrid CBC Insp +wavelets analysis by underpredicting the strength of the signal. The hybrid CBC Insp +wavelets agrees well with the full CBC IMR , though the uncertainty of the former is larger in the merger phase where the signal is no longer modeled with CBC templates but with wavelets.
The recovered source parameters are given in the right panel of Fig. 2. The CBC IMR analysis leads to consistent results with those reported in previous studies [131][132][133], while the CBC Insp analysis shows a significant bias away from the expected posteriors. The bias is most evident in the chirp mass: PhenomDNRT compensates for the lack of a merger and ringdown by decreasing the chirp mass, thus leading to a longer inspiral phase in an attempt to capture part of the missing merger cycles. Despite this bias, the CBC Insp analysis is still not able to capture the full signal as seen on the left panel. The CBC Insp +wavelets proof-of-concept model results in recovered parameter posteriors that are consistent with the full CBC IMR analysis. The posteriors are not expected to be identical as now the CBC model has access only to the inspiral portion of the signal and therefore lacks information available to the CBC IMR analysis. However, the wavelet model can efficiently capture the missing portion of the signal, allowing the CBC model to recover unbiased system parameters from the inspiral phase only. Figure 3 examines the CBC Insp +wavelets analysis in more detail and the complementary roles of the components in the hybrid model. We plot the signal reconstruction from each submodel as well as their sum. The sum of the individual model components is the same as the CBC Insp +wavelets model in Fig. 2. Figure 3 shows that as expected the BBH inspiral phase is primarily reconstructed by the CBC PhenomDNRT model. As we move from the inspiral toward the merger time, the CBC model starts tapering off and deviating from the full GW signal. At the same time, the wavelet model captures the part of data that is not covered by the CBC model, such that the combination of CBC and wavelets accurately models the full GW signal. This toy model for GW150914 demonstrates the main concept of our hybrid analysis: the full GW signal is modeled with a well-understood component, i.e. the CBC waveform, and a model-agnostic component, i.e. the wavelet model, that covers any parts of the signal that are not included in the CBC model.

IV. CONSTRAINTS ON THE GW170817 POSTMERGER
We turn our attention to GW170817, the first GW detection of a BNS coalescence [5,30,134]. Both the low-frequency inspiral phase and a potential shortduration high-frequency postmerger phase of the signal have been studied separately in detail, with the latter remaining undetected [30,72]. Here we apply our hybrid CBC+wavelet model to study both parts of the GW170817 signal simultaneously and compare our results with separate analyses of the inspiral and postmerger signals.
We analyze 64s of data from LIGO Hanford and LIGO Livingston around GPS time 1187008882.446 [129,130] where the prominent glitch in LIGO Livingston has already been subtracted [94,135]. While Virgo data are available for that time and aided in constraining the source sky location [30], we do not consider them due to the lower sensitivity and the fact that we fix the sky location to the known values of α = 3.446 rad and δ = −0.408 rad [136,137]. We also use a fixed PSD rather than marginalize over the noise model for computational efficiency. We perform three analyses with different models and data: 1. A "CBC Insp " analysis with data in the frequency range (32,2048)Hz that are described only by the CBC model with the PhenomDRT waveform thus focusing on the inspiral signal. For both analyses with the wavelets model, we use a prior on the number of wavelets of D ∈ [2, 100]. The reason for selecting a minimum number of 2 wavelets is different for each analysis. In the case of the hybrid analysis, the CBC part of the model does not terminate at merger, but smoothly tapers into the postmerger phase for a few milliseconds, see Fig 9. We therefore need at least two wavelets, one to undo this effect of the CBC model and the other to capture the true postmerger signal. In the wavelet-only study we again need at least two wavelets such that one wavelet can capture the merger itself which extends into the analysis bandwidth and the other wavelet can capture the contribution from the postmerger signal. The full settings of our three analyses are detailed in Table II. We compare the two analyses that include the inspiral of the signal in Fig. 4 which shows the signal reconstructions focused on the end of the signal (left panel) and the marginalized posterior distribution for selected source parameters (right) of GW170817. Despite the high SNR of the signal, the amplitude is relatively weak compared to the noise level, however both analyses identify the signal and lead to consistent reconstructions for the late  Fig. 4 and Fig. 5, the segment length, the sampling rate, the low frequency cut off, the models active, the CBC waveform, the wavelet window around the trigger time, the wavelet bandwidth, and the minimum number of wavelets. inspiral phase as shown on the left panel. As the binary approaches merger the reconstruction uncertainty increases due to the decreasing detector sensitivity at increasing frequencies. This behavior is evident in both analyses, suggesting that the GW170817 postmerger signal remains undetected as expected from the detector sensitivity at the time of GW170817 [30]. The right panel of Fig. 4 examines selected CBC parameters recovered from the inspiral signal. The recovered parameters are consistent with previous results using the same CBC model [30]. The only difference is that ourΛ posterior has more support for the higher of the two modes present in the results of [30]. This can be attributed to different uses of prior: we use of a prior that is flat inΛ − δΛ with the constraint that Λ 1 > 0 and Λ 2 > 0. This choice results in a marginalized prior forΛ that favors larger values. In contrast, [30] uses a flat marginalized prior forΛ. We have verified that reweighting our posterior to a flat marginalized prior gives consistent results with [30]. Furthermore, the hybrid CBC Insp +wavelet analysis leads to essentially identical results for the CBC parameters as the traditional CBC Insp analysis. This is consistent with [74] that finds that failure to account for the postmerger signal will not lead to biases in the CBC parameters for typical detector sensitivities, but also serves as a sanity check of the CBC Insp +wavelet analysis and the fact that sampling for the joint CBC+wavelet model has converged.

A "CBC
After confirming that the CBC Insp +wavelet analysis matches the CBC Insp analysis as far as the inspiral portion of the signal is concerned, we switch to comparing the two analyses that include a wavelet model for a possible high-frequency postmerger signal. In Fig. 5 we show the reconstructed GW spectrum for CBC Insp +wavelet and wavelet-only. The vertical band gives the 90% credible interval for the merger frequency estimated using Eq. 11 in [138] given the binary component masses and tidal parameters from the CBC model posteriors. This merger frequency is also the location at which the Phe-nomDNRT template starts tapering off, so it is related to the transition between the CBC and the wavelet models.
At frequencies below ∼ 1000Hz the CBC Insp +wavelet reconstruction is dominated by the inspiral signal which is clearly detected as the spectrum lower limit is nonzero. With increasing frequency the detector sensitivity decreases, resulting in larger uncertainties in the recon- structed spectrum. Finally, at around the merger frequency, 1600Hz, the reconstruction uncertainty is large and consistent with no detected signal. A similar picture is drawn by the wavelet-only analysis which again places only upper limits on the signal in frequencies above 1024Hz. The upper limit from the CBC Insp +wavelet and the wavelet-only analyses is comparable, though the former is consistently lower across the frequency band, indicating more stringent constraints on the presence of a postmerger signal. Additionally, the wavelet-only analysis does not lead to a detection of the merger signal in the frequency range (1000,1500)Hz, unlike the CBC Insp +wavelet case. This suggests that this highfrequency portion of the signal is not individually detectable, but only inferred coherently from the preceding inspiral signal.

V. SIMULATED SIGNALS
Going beyond the GW170817 upper limits, in this section we study simulated BNSs where the postmerger emission is detectable by future GW detectors. We assume a network of two detectors at the current location of LIGO Hanford and LIGO Livingston and a zero noise realization. Since all signals have high SNRs we assume that an electromagnetic counterpart has been identified and the sky location of the source is known. Following [31] we do not explicitly select any planned GW detector such as Cosmic Explorer [139][140][141] or the Einstein Telescope [142,143] and their nominal sensitivity. We instead work with the Advanced LIGO design sensitivity [144] and gradually lower the noise PSD, emulating improving detector sensitivity and higher signal SNR.
Since no NR simulation of the full BNS signal as observed by ground-based detectors exists, all simulated signals are constructed in a hybrid fashion: the postmerger signal is obtained through NR simulations, while the premerger signal is computed with PhenomDRT using the same parameters as the NR simulation. We choose extrinsic parameters consistent with GW170817, namely distance D L = 40Mpc, inclination ι = 2.635, polarization angle ψ = 0, and the known sky location. All signals are analyzed with the full CBC+wavelets model with run settings equivalent to the CBC Insp +wavelet analysis from Table II.

A. Waveform from Kastaun and Ohme
The first full BNS waveform we consider was constructed and released by [145]. This hybrid waveform was constructed by combining the PhenomDNRT model for the inspiral and the results of NR simulations for the postmerger GW signal. The hydrodynamical simulations assume a GW170817-like initial system with the hadronic SFHO EoS [146], resulting in a short-lived hypermassive NS. We use the resulting hybrid waveform for a BNS system with mass ratio q = 0.9, detector frame component masses of m 1 = 1.438M and m 2 = 1.294M , and dimensionless tidal parameters Λ 1 = 280 and Λ 2 = 551, giving a binary tidal parameterΛ = 396 and merger frequency f merger = 2110 Hz.
We inject the signal in a detector network with sensitivities 2xDS, 4xDS, 6xDS which respectively denote 2, 4, and 6 time improved strain sensitivity compared to the LIGO design sensitivity (i.e., the design sensitivity divided by 2, 4 and 6). In terms of amplitude spectral density, they correspond to 8.23×10 −24 Hz −1/2 , 4.12×10 −24 Hz −1/2 , and 2.72×10 −24 Hz −1/2 at 3326 Hz, which is the frequency of the main postmerger mode of the signal, f peak . The resulting network SNRs are 142, 384, and 425 for the premerger signal (f < f merger ), and 2.8, 5.6, and 8.4 for the postmerger signal (f ≥ f merger ). Figure 6 shows the whitened time-domain data and reconstruction for the injection at 6xDS, again focusing around the late stages of the signal. The full reconstruction with the combined CBC+wavelets model accurately captures the entire signal. We also plot the individual components of the full model, namely the CBC and the wavelets model separately. As expected, the CBC model captures everything up to merger and then tapers off. The wavelets overlap with the CBC model in the taper region around merger, effectively canceling out the ringdown-like oscillations of the CBC model that do not match the data. At later times, the wavelet model extends to the postmerger part of the signal, capturing its main oscillatory component. Figure 7 shows the reconstructed spectrum from injections on different detectors sensitivities. As the detector sensitivity increases, the signal reconstruction becomes more accurate. The premerger signal is recovered in all analyses given its strength, and increasing sensitivity reduces the reconstruction uncertainty. The postmerger signal, on the other hand, is too weak to be detected in the 2xDS case, and therefore the reconstruction is essentially uninformative and similar to the GW170817 one. At increasing detector sensitivity, the postmerger signal emerges from the noise. At 4xDS the reconstructed spectrum starts identifying the main postmerger peak though only at the ∼ 50% credible level. At 6xDS the reconstructed spectrum not only confidently identifies the main postmerger peak but also has evidence of a subdominant peak at around 2337 Hz.
Posteriors for select premerger and postmerger parameters are given in Fig. 8. The peak frequency of the postmerger signal f peak is defined in the same way as [31,85]: for each reconstruction posterior sample we compute the frequency at the maximum of the spectrum after the merger. If the spectrum possesses no maximum, then a sample is drawn from the f peak prior. All posteriors are consistent with the injected values, and uncertainties decrease with increasing detector sensitivity as expected. The f peak posterior at 2xDS is essentially the prior, consistent with the fact that the postmerger signal was not detected. At 4xDS the f peak posterior starts exhibiting a peak at the injected value, consistent with the partial identification of the postmerger signal in the middle panel of Fig. 7. At 6xDS the peak frequency is accurately measured.

B. Hadronic EoS
We construct further simulated signals based on the postmerger NR simulations from [31], which employ the conformal flatness approximation [147,148]. Compar-isons to fully relativistic studies show a very good agreement with regards to postmerger GW frequencies but an underestimation of the postmerger GW amplitude by some 10% as a result of using the quadrupole formula for GW extraction [47]. Since our analysis does not rely on calibration to any NR simulations, we expect our results below to be unaffected by such uncertainties. Indeed, BayesWave has been shown to produce reliable results on simulated postmerger signals made with different NR codes [81,85].
All simulations in [31] were at the time constructed to be consistent with GW170817, though in light of new data from NICER [13][14][15][16] some of the softest EoSs there are now disfavored. We work with EoS3 from [31], a hadronic EoS that is consistent with radii values inferred in [26] and corresponds to f peak = 2880Hz and R 1.4 = 12.6 km. The inspiral portion of the signal is again described by PhenomDNRT for a GW170817-like system with detector frame component masses m 1 = m 2 = 1.362M , dimensionless tidal parameters Λ 1 = Λ 2 = Λ = 587, and zero spin.
The full waveform is constructed by aligning the projected premerger and postmerger waveforms in an overlap interval in the time domain. We use a transition window [t 1 , t 2 ] where the full waveform transitions from the Phe-nomDNRT template to the NR simulation data. The full waveform is then (22) where h inspiral denotes the PhenomDNRT inspiral model, h NR denotes the postmerger NR simulation, A scale is a scale parameter that can be varied to control the strength of the postmerger signal, and x ensures a smooth transition through a Planck taper function We fix the detector sensitivity to 2xDS and vary the amplitude of the postmerger signal through the scale parameter A scale . The main motivation for this is that it allows us to vary the postmerger signal strength, while keeping the premerger SNR (and thus the computation cost) manageable. This scaling further allows us to address the underestimation of the signal amplitude within the conformal flatness approximation mentioned above. Larger sets of simulations do not find a very tight relation between the postmerger amplitude and binary parameters [79], and simulations in general may over-or underestimate the postmerger amplitude to some extent.
In what follows, we keep the premerger SNR constant at 112, while the postmerger SNR varies and results are presented as a function of the ratio of the premerger to the postmerger peak amplitudes. We explore three values for the ratio of the peak postmerger to the peak premerger amplitudes: 0.7, 0. 9   A scale = 1.5, 2.0 and 2.5 and postmerger SNRs of 6.2, 8.1 and 10.0, respectively). We show an example of the above construction process in Fig. 9.
In Figure 10 we compare the reconstructed spectra from our analyses for the three cases of post/premerger amplitude scaling. The inspiral portion of the signal and the corresponding reconstruction are similar in the three panels. As the amplitude of the postmerger signal in- creases from left to right the reconstructed spectrum includes more detailed features, progressing from hints of a postmerger peak on the left panel to increasingly more confident identification in the middle and right panel. ically we use the fit of [85] that holds for hadronic EoSs and the binary mass and tidal deformability as extracted from the premerger signal to compute the expected f peak . The expected and recovered posterior values for f peak agree, showing that our analysis can correctly conclude that the premerger and postmerger signals are consistent  [26] using data from heavy pulsars, GWs, and the Miller et al. [13,15] radius results using X-ray data. The adopted EoS is inconsistent with this posterior at the 90% level for masses above ∼ 1.8M , but it is consistent to the same level with the less constraining Riley et al. [14,16] results.
with each other given expectations from NR simulations of hadronic EoSs.

C. EoS with phase transitions
Since the hypermassive NS that gives rise to the postmerger signal is characterized by higher core densities and temperatures than the premerger NSs, the postmerger signal has the potential to reveal new high-density physics. One possibility is a strong phase transition in the EoS toward quark degrees of freedom 3 . Such phase transitions could occur at lower densities and thus be detectable with premerger data only [151][152][153][154][155][156][157][158][159][160]. However, as already mentioned above, the tidal deformability is a steeply decreasing function of the NS compactness, making premerger data less constraining about high-mass NSs that could contain quark cores. postmerger data, on the other hand, can be used to probe higher densities and NR simulation-based studies have explored the potential signature of a high-density phase transition on the GW signal [62][63][64][65][66][67][68].
The most prominent signature of a sufficiently strong phase transition would be an increase in the postmerger peak frequency as the softening of the EoS would lead to a more compact merger remnant [63]. It has been proposed that such a frequency increase could be identified if one compares the premerger and postmerger data: for hadronic EoSs the tidal deformability and peak frequency follow the approximately EoS-insensitive relation shown in Fig. 17. This relation could be violated for EoSs with high-density phase transitions as the tidal deformability is determined by the hadronic part of the EoS alone, while the postmerger frequency is affected by the phase transition [63,66] 4 . The degree of deviation from the EoS-insensitive relation, and thus how detectable it is, depends on the strength of the transition [58,63].
As a first example of such a signal, we work again with EoS3 from [31], but the premerger data are constructed with tidal parameters that are systematically shifted compared to their hadronic EoS value. This effectively results in a stiffer hadronic EoS that undergoes a phase transition which softens the postmerger signature to the level of EoS3. All intrinsic parameters remain to be the same as the simulations of Sec. V B with the exception ofΛ = Λ 1 = Λ 2 = 800.The postmerger SNRs remain the same as Sec. V B.
Reconstructed spectra and parameter posteriors are shown in Figs. 13 and 14 respectively. We obtain qualitatively similar results to Figs. 10 and 11 with the main difference being that now the premerger and postmerger are now inconsistent as expected. Figure 14 shows the expected f peak value given the premerger constraints on the binary properties under the assumption of a hadronic EoS. The recovered f peak is inconsistent with this expectation to within its measurement uncertainty, signaling 3 We use the term "strong" phase transition to emphasize that the exact characteristics of the hadron-quark phase transition are currently unclear, and only a sufficiently strong transition, for instance a first-order phase transition with large latent heat, may impact the merger dynamics such that the occurrence of quark matter leads to an unambiguous signature. See, e.g., [58] for a detailed discussion. In fact, quark matter may also resemble the behavior of purely hadronic matter (commonly referred to as the masquerade problem [149,150]), which may not leave a very prominent imprint on the GW signal. 4 In fact, a postmerger frequency deviation may also occur if quark matter is present before the merger [66].
the presence of additional high-density effects in the EoS that affect the postmerger signal. We further assess how well our hybrid analysis could detect a strong phase transition by simulating a signal with the DD2-SF-4 EoS [161,162] and corresponding postmerger simulation from [63], also performed with the conformal flatness approximation. The mass-radius relation for this EoS is given in Fig. 12 and it exhibits the characteristic radius reduction due to strong phase transitions starting at ∼ 1.5M . Given the simulated binary masses of 1.35M for both components, the inspiral signal is emitted by hadronic NSs, while the postmerger signal is affected by the onset of the phase transition [63]. For reference, we also show the 90% symmetric credible intervals of the mass-radius posterior for the EoS derived in [26] using heavy pulsar, GW, and X-ray data. The DD2-SF-4 EoS is inconsistent with the posterior at the 90% level for large masses as it underpredicts the radius of a 2M NS, though it is consistent with current data at the 95% level. If we instead used the less constraining radius results from Riley el at. [14,16], the EoS would be consistent with the posterior at the 90% level. Results are presented in Figs. 15 and 16, where again we find that the postmerger signal can be reconstructed for sufficiently loud signals, and the inconsistency betweeñ Λ and f peak can be identified.
Finally, we elaborate on the premerger and postmerger consistency in Fig. 17. The gray band shows the expected relation betweenΛ and f peak assuming hadronic EoSs. We compute this relation by converting the R − f peak EoS-independent 5 fit by [85] to aΛ − f peak relation. We convert R to Λ using the scaling relation calibrated for GW170817-like systems for our source mass parameters M tot and M given by Eq. 17 in [118]. Overplotted are the two-dimensional posteriors forΛ − f peak for the case of a hadronic EoS (EoS3), the phase-transition version of the EoS3-based signal, and the DD2-SF-4 EoS again with a strong phase transition. The injected values for the hadronic case is consistent with expectations and this is confirmed by the recovered values to within the statistical error. On the other hand, a strong phase transition leads to aΛ − f peak combination that is inconsistent with hadronic expectations and the extracted posteriors are able to identify this behavior.

VI. CONCLUSIONS
We present a hybrid approach to study the full GW signal emitted during a BNS coalescence, including a possible postmerger component. Our method models the inspiral part of the coalescence with waveform templates as implemented in LALSimulation and the postmerger  FIG. 14. One-and two-dimensional marginalized posterior distributions for selected pre and postmerger parameters for the analyses of Fig. 13. The parameters shown are the sourceframe chirp mass M and total mass M , the tidal deformation parameterΛ and the postmerger peak frequency f peak . Black vertical black lines or crosses denote the true values of each parameter. The gray region is the expected value for f peak given the premerger inferred parameters and the fit of [85] that assumes hadronic EoSs. The premerger and postmerger results are now inconsistent with expectations for hadronic EoSs.
part of the signal with a superposition of wavelets. We do not impose phase coherence in the transition between the template and the wavelets, however, the full signal is smooth and coherent as the template model used already extends coherently past merger and into the postmerger part of the signal.
Applying our method to GW170817, we demonstrate that the inspiral parameters are consistent with previous results. We do not detect a postmerger signal for GW170817 but find that the high-frequency portion of the signal between 1000-1500 Hz is only detected using the full hybrid analysis in contrast with the traditional postmerger-only analyses. We apply our method to simulated signals with a detectable postmerger component and show that our full analysis simultaneously reconstructs both the inspiral tidal deformation and the postmerger dominant frequency peak when the SNR is sufficiently high.
In this work, we sample the inspiral tidal parameters and postmerger peak frequency independently, imposing no relation between them and no assumption on the nature of the EoS. This allows us to detect a possible signature of the hadron-quark phase transition in NS mergers, as it lead to a characteristic frequency shift of the postmerger signal. A possible extension of this analysis could use information extracted from a premerger signal to guide the analysis of the postmerger signal. If one restricts to a hadronic EoS, relations between the premerger tidal parameters and postmerger peak frequency could be used to predict the approximate location of the dominant postmerger peak from the inspiral tidal parameters. This could improve the prospects of detecting the postmerger signal in weaker signals, at the expense of assuming the EoS is similar to a hadronic one. A further possible improvement concerns the use of "chirplets", sine-Gaussian wavelets with an evolving frequency [95], which might be better suited for a time-evolving postmerger frequency mode. of postmerger signals show subdominant peaks in the GW spectrum, which can also be used to characterize the properties of the merger remnant star [49,53,163]. Our analysis can identify some of these peaks for sufficiently loud signals, however they are typically less loud than the main peak and can be lost in the noise. Numerical simulations suggest that the main and the secondary peaks in the spectrum are related to each other [53,82,84]. In future work, this feature can be utilized in the form of a prior that links the two peaks and is parametrized in terms of the remnant compactness. Such a prior would enhance the sensitivity of the analysis to secondary peaks as it effectively adds modeling information about the signal similar to the approach of analytic models [76][77][78][79][80][81][82]. Our analysis concerns premerger SNRs in the hundreds, where systematic biases in the waveform models could be important and would have to be mitigated in the lead-up to such improved detectors [74,[164][165][166][167][168]. However, as we make minimal assumptions about the postmerger waveform and use sine-gaussian wavelets, we expect the reconstructed signal to be less susceptible to systematics than analytic models [76][77][78][79][80][81][82] that are limited by the accuracy of the NR simulations on which they are based (see for example Fig. 4 of [80]). Our analysis instead is not limited by NR accuracy and has yielded unbiased results on simulated signals form different codes [81,85]. On the other hand, analytic models can be more straightforwardly attached to premerger waveform templates and enforce phase coherence through merger [80].
Next generation detectors are expected to detect thousands of BNSs [169][170][171], the majority of which will be weak with individually undetectable postmerger signals.
While we study single loud sources here, our method can also be applied to the expected numerous weaker BNS sources which might dominate the overall constraints when combined [172]. Though our analysis does not hinge on the existence of accurate models for the postmerger signal, it is possible that biases in the premerger waveform will lead to biased inferences about the postmerger when making use of EoS-independent relations connecting them. The extremely sensitive observations possible with next-generation detectors indeed require control over a wide range of potential systematic biases, and the flexible analysis presented here helps mitigate such biases from the postmerger signal.

ACKNOWLEDGMENTS
We thank Will Farr and Wynn Ho for many useful discussions. We also thank Sophie Hourihane and Tyson Littenberg for discussions and assistance about BayesWave. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes. This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation. The authors are grateful for computational resources provided by the LIGO Laboratory and  Fig. 15. The parameters shown are the detector frame chirp mass M and total mass M , the tidal deformation parameterΛ and the postmerger peak frequency f peak . Vertical black lines or crosses denote the true values of each parameter. The gray region is the expected value for f peak given the premerger inferred parameters and the fit of [85] that assumes hadronic EoSs. The premerger and postmerger data are now as expected inconsistent with expectations for hadronic EoSs. Comparison of the recoveredΛ-f peak posterior distribution from the EoS3 hadronic and PT data as well as the data with the DD2-SF-4 EoS. The gray bands corresponds to the EoS-independent fit relating the expected values ofΛ-f peak for hadronic NSs. We show 50% (dashed) and 90% (solid) contours, while colors correspond to those of Figs. 11, 14, 16. In all cases, we find the correct agreement or disagreement with the expected hadronic relation for sufficiently loud signals.