1 Introduction↩︎

The detection of GW150914 [1], which was the first observed gravitational wave (GW) event, initiated the beginning of GW astronomy. With the completion of three observing runs (O1, O2 and O3) by the ground-based detectors advanced LIGO [2] and advanced Virgo [3], recently joined by KAGRA [4], [5] more than 90 GW events are included in the Gravitational Wave Transient Catalog (GWTC) [6]–[11] (see [12] for additional candidates), two of which, GW170817 [13] and GW190425 [13], [14], were identified as Binary Neutron Star (BNS) mergers. This has opened new avenues for probing the physics at the extreme densities encountered at the center of compact stars (for comprehensive reviews, see, e.g. [15]–[19] and references therein). The detection of new BNS merger signals in the coming observational runs (such as the ongoing O4 run) is expected to further refine our understanding of compact stars and their Equation of State (EoS) (e.g., see [20]–[35] for how we could benefit from further observations).

GW170817 was the first GW event that allowed us to constrain the EoS of neutron stars (NS) and rule out some very soft models of high-density matter [13], [16], [36]–[43]. Moreover, this event marked the beginning of the era of multi-messenger astronomy, since it was detected with both Gravitational and Electromagnetic (EM) radiation [44], [45], spanning from radio, through optical and x-rays [46] to gamma rays [47], [48]. More such multi-messenger events are expected to be observed in the future observing runs (e.g. [49]) and provide us with the opportunity to combine different information for the same event[50]–[59].

Despite the improvement in the sensitivity of current detectors to the O4 level, there are no high hopes of detecting the post-merger part of the GW waveform of BNS mergers during O4 [60]–[64]. However, such signal could be detected with future detector upgrades, beyond the fifth observing run (O5) and more likely with the planned third-generation detectors [65]–[70]. The importance of such post-merger waveform signals acts as a motivation for dedicated observatories, such as the one suggested in [71], which we will henceforth call the High Frequency (HF) detector design.

The detection of GWs in the post-merger phase of BNS mergers will allow us to place constraints on the NS EoS in a density regime that cannot be probed directly by extracting the tidal deformability of NS in the inspiral phase. Moreover, the post-merger phase is rich in additional physics that could be probed (high temperature, shock waves, magneto-rotational instabilities, unstable oscillation modes, etc.).

Simulations of BNS mergers have shown [72]–[74], that the remnant can either promptly collapse to a black hole (BH), survive for a short duration (of order several or tens of milliseconds) before a delayed collapse takes place, or even avoid collapse all together and remain as stable NS, if the component masses are sufficiently small and the EoS sufficiently stiff. In the post-merger phase, the highly deformed remnant continues for some time to radiate GWs, if it is in a quasi-stable state or a stable state. There are different studies, such as Ref. [75], [76], that implement techniques to predict whether the merger will form a remnant or will promptly collapse into a black hole based on the gravitational-wave signal.

The strongly differentially rotating remnant can avoid prompt collapse if the total mass of the binary is above a threshold mass. The remnant then becomes a transient Hypermassive NS (HMNS) [77]–[81]. The merger process excites linear nonaxisymmetric oscillation modes, nonlinear combination tones, and other transient effects in the remnant; see e.g. [78], [80], [82]–[91] and references therein. The main oscillation mode is the \(l=m=2\) \(f\)-mode, the frequency of which is typically denoted as \(f_{\rm{peak}}\) in the literature on post-merger GW spectra¹.

The excitation of these modes leaves a unique signature in the GW spectrum, typically between \(\mathrm{2~kHz}\) and \(\mathrm{4~kHz}\), which depends on the masses of the two binary components and the EoS, see, e.g.  [16], [73], [78], [83]–[85], [92]–[97] and references therein. Thus, the oscillation modes are closely related to the nature of the resulting NS remnant, and by analyzing the information of the GW spectrum, one can infer various characteristics of the NSs, e.g.  [61], [92], [98]–[101] and [85] for a review and references therein.

For example, there are empirical relations that correlate \(f_{\rm peak}\) with the radius of the non-rotating NS in the inspiral phase [78], [84], [85], [95], [102], [102]–[105]. A recent study  [106] demonstrates how to employ such empirical relations to extract \(f_{\rm peak}\), when combining information from several subthreshold events. The hierarchical Bayesian method proposed in [106] overcomes the limitations posed by a low signal-to-noise ratio (SNR) in retrieving \(f_{\rm peak}\).

Apart from stable oscillations excited during merger, rotational instabilities could also develop in the remnant (see [107] and references therein). Of particular interest are the dynamical shear instabilities, often manifested as the low\(-|T/W|\) rotational instabilities  [108]–[110] where \(T\) is the rotational kinetic energy and \(W\) the gravitational binding energy of the star. Differentially rotating stars can develop such instabilities for relatively low values of the rotational parameter \(\beta = |T/W|\), [108], [111]–[113]. In particular, one condition for the low\(-|T/W|\) instability to set in, is that the pattern speed of the \(l=m=2\) \(f\)-mode matches the rotational speed of the star [109], [110], [114]–[117].

Due to the instability’s nature, depending on the rotational parameter, the impact of different rotation laws of HMNSs on the low\(-|T/W|\) rotational instability was studied by [109] and they attempted to give a qualitative explanation of its features observed in numerical simulation of either BNS remnants (e.g., the Ref. [118], [119]) or rapidly rotating cold NSs (e.g., the Ref. [115]). It was found that in a differentially rotating BNS merger remnant, an oscillation mode could co-rotate with the matter at two different points. The parameter \(\beta\) was found to be as low as 0.02, and as this parameter increases, the oscillation mode stabilizes to a very specific value. Some numerical simulations of BNS mergers (see, e.g. [118], [120]) include evidence for the re-excitation of the \(l=m=2\) \(f\)-mode in the remnant.

In this study, we focus on two different cases of post-merger waveforms from the study of [120], where the MPA1 EoS [121] was used. We explore whether rotational instabilities could be detected with future upgraded or 3rd-generation networks of detectors. For parts of our analysis, we use the BayesWave pipeline [122]–[124], a Bayesian data-analysis algorithm useful to recover a GW signal via sine-Gaussian wavelets, which has already been extensively used for BNS post-merger studies (see, e.g., Ref. [62], [63], [125]–[128]). We inject the simulated signals into colored Gaussian noise corresponding to different detectors. Our sources are considered to be optimally oriented with respect to the detector located in Livingston or with a non-optimal inclination and sky position at distances of at least² 40 Mpc. In order to infer the capabilities of detecting rotational instabilities in BNS merger remnants, we calculate the overlap for the instability-part only.

In the next sections, we demonstrate the effects of the low\(-|T/W|\) rotational instability (Sec. 2) and define the different network configurations (Sec. 3), with which we aim to detect the instability. In Sec. 4, the injections of our investigation are presented, and in Sec. 5 a description of our analysis method is given, which is based on the BayesWave algorithm. The simulation setup and evaluation process are presented in Sec. 6. Finally, our results are provided in Sec. 7 and discussed in Sect. 8.

2 Post-merger waveform↩︎

In order to constrain the detectability of the low-\(|T/W|\) rotational instability, we employ numerically generated post-merger waveforms, which reveal the re-excitation of the \(m=2\) \(f-\)mode. More specifically, we study cases of NS with equal masses, which were simulated for the MPA1 EoS [121] using the Einstein Toolkit software [129] in [120]. We focus on the cases of 1.50+1.50 \(M_\odot\) and 1.55+1.55 \(M_\odot\) equal-mass mergers.

Fig. 1 shows the GW signal³ (gray) of a BNS merger at 40 Mpc with MPA1 EoS and a total mass (\(M_{\rm{tot}}\)) of 3.1 \(M_\odot\) from [120]. For the first \(O(10)\) ms after merger, the GW signal is decaying. However, after this initial period, the rotational instability sets in and causes a re-excitation, which maintains appreciable amplitude, showing only a slower damping rate with time. In the upper panel of Fig. 1, the waveform is split into the signal from the merger up to the onset of the instability (dashed blue curve) and the subsequent instability part (green), when the instability develops. In the lower panel of Fig. 1, we observe that the Amplitude Spectral Density (ASD) of the full signal has an additional narrow peak, reaching about twice the ASD of the signal truncated at the time of the onset of the instability (dashed blue curve). The main focus of the present work is to investigate whether the presence of such an additional narrow peak could become detectable with different GW detector networks.

3 Network Configuration↩︎

Our study focuses on the detectability of the low-\(|T/W|\) rotational instability in post-merger NSs by future observatories. For that reason, we consider three cases of network configuration, for which a summarized description is given in Table 1. The first network configuration corresponds to twice the projected sensitivity of the fifth planned observing run (O5) Advanced Virgo [3], [130], and \(A+\) (LIGO) [2], [130]).

The second network is more sensitive for post-merger NS signals and consists of the planned Cosmic Explorer (CE) [131], [132] and the planned Einstein Telescope (ET) [133], comprising the 3G detectors, to become operational within the next decade. For both detectors, there have been different estimates for their sensitivities, but in this study we choose to work with the “CE 20 km" for the post-merger curve [131], [134] and the”ET-D" [133], [135], because these designs have greater sensitivity to post-merger NS signals. Regarding their locations, the CE detector is assumed to be located in Livingston and ET in Cascina [136].

The third configuration refers to the recent proposal for a 25 km L-resonator interferometer [71], denoted as HF (High Frequency), which we assume the LIGO-Livingston location. A significant limitation in future detectors, which are based on dual-recycled-Fabry-Pérot-Michelson Interferometers, is the loss in the signal extraction cavity (SEC). The HF design aims to suppress this loss-limit at high frequencies, resulting in a larger than 1/yr detection rate of the post-merger signal [71], significantly better than other 3G detector designs. The HF design includes an “L" shape optical resonator of 25 km arm length and laser wavelength at 1064 nm (for a detailed description, see Sec. IV in Ref. [71]). Their proposed configuration targets post-merger signals of BNSs (sensitive between 2-4kHz) with a peak sensitivity at 3kHz. If realized, the HF detector would be an ideal interferometer for investigations focusing on the NS post-merger signal, since the peak frequency of such signals is between 2kHz and 4kHz. For our chosen models, the peak frequency of the post-merger signal is close to 3 kHz [120].

We investigate source cases that are either optimally oriented with respect to the Livingston detector or randomly oriented. The antenna pattern response functions \(F_+, \,F_\times\) of the interferometer as defined in the sky plane, are \[\begin{align} F_{+,\times} =& -\frac{1}{2}\big(1+\cos^2\theta\big)\cos2\phi \cos2\psi \\ & \mp \cos\theta\sin2\phi \sin2\psi, \end{align} \label{eq:responses}\tag{1}\] with (\(\theta, \, \phi\)) denoting the sky-location of a source, relative to the axes of the detector arms. The two polarization components \(h_+\) and \(h_\times\) are with respect to axes in the plane of the sky that are rotated by an angle \(\psi\) (see Fig. 3 of Ref. [137]). Then, one can define the right ascension of the source, \(\alpha=\phi +\)GMST, and its declination, \(\delta=\pi/2 - \theta\), where GMST is the Greenwich mean sidereal time of arrival of the signal. Therefore, the quantities \(\alpha,~\delta\) maximize the responses \(F_+,F_{\times}\), when \(\phi\) is equal to the detector longitude (\(-1.58430937078\) rad at the location of the LIGO Livingston detector) and \(\theta\) the detector latitude (\(0.53342313506\) rad). Furthermore, we assume that the binary is seen face-on, i.e. the inclination is \(\iota=0\) (see [138]). With these assumptions, the source is considered to be optimally oriented with respect to the detector.

4 Injections↩︎

For each BNS merger case, we construct two different source cases to perform the injection, differing only in their sky localization and orientation. One case has a non-optimal orientation and inclination (using the inferred values for the GW170817 event [13], [44], [139], [140]). The second case is an optimally oriented source with respect to the LIGO-Livingston detector (denoted as L). The parameters of the injected model are given in Table 2. The distance from the source depends on the detectability of the network configuration (see Table 1) and ranges from 40 \(\mathrm{Mpc}\) to 200 \(\mathrm{Mpc}\). The lower limit of 40 \(\mathrm{Mpc}\) is set due to the inferred distance of GW170817.

We assume two different cases of total mass \(M_{\rm{tot}}\) 3.0 and 3.1 \(M_\odot\). In Ref. [120], a re-excitation of the GW signal is demonstrated, which is a characteristic of the presence of rotational instabilities. Although the inclusion of magnetic fields would affect the rotational profile, as a first step, we consider simulations without magnetic fields or other viscous effects.

The duration of the simulated signal is 74 \(\mathrm{ms}\), which is then injected into 1 second of colored Gaussian noise. The frequency band of the analysis is (1800, 4096)⁴ Hz for the case of \(M_{\rm{tot}}=3.0M_\odot\) and (2000, 4096) Hz when \(M_{\rm{tot}}=3.1M_\odot\). The lower cutoff frequency is chosen with the aim of leaving out parts of the spectrum that are influenced by the inspiral signal.

5 Data Analysis methodology↩︎

In order to assess the detectability of the various signals that appear with different characteristic spectra, we adopt a template-agnostic method in a Bayesian framework. In this framework one begins with the Bayes’ Theorem, which is expressed as \[p(\vec{\lambda} |s, \mathcal{M})=\frac{p(s|\vec{\lambda} , \mathcal{M})p(\vec{\lambda} , \mathcal{M})}{p(s, \mathcal{M})}, \label{Bayes32theorem}\tag{2}\] where \(p(\vec{\lambda} |s, \mathcal{M})\) is the posterior for the parameters \(\vec{\lambda}\) of the signal \(h\) given the data \(s\) and model \(\mathcal{M}\); \(p(\vec{\lambda} , \mathcal{M})\) is the adopted prior of the parameters of interest; and \(p(s|\vec{\lambda} , \mathcal{M})\) is the likelihood function of the data, given the particular signal. \(p(s, \mathcal{M})\) is the marginalized likelihood, or evidence (integral of the likelihood function over the parameter space) \[\begin{align} \mathcal{Z}&=\int p(\vec{\lambda} ,s|\mathcal{M})d\vec{\lambda} \\ &=\int p(s|\vec{\lambda} , \mathcal{M})p(\vec{\lambda} | \mathcal{M})d\vec{\lambda} , \end{align} \label{evidence}\tag{3}\] which, in parameter estimation procedures serves solely as a normalizing constant [141].

However, the evidence is crucial for model selection purposes, because it essentially characterizes the capabilities of our adopted model to describe the measurements [62], [122], [142]. Having the evidences of two competing models, one can compute the Bayes Factor as the ratio of the evidences and infer which model is supported best by the measured data. In practice, the integral of Eq. (3 ) is quite challenging to compute directly, but one can adopt numerical methodologies to approximate it (e.g. via Thermodynamic Integration [143] or for different approaches, see [144]). However, there are situations where the set of models to be tested is relatively large, and computing the evidence for each combination of models can be quite inefficient or even computationally prohibiting. Instead, a trans-dimensional Markov Chain Monte Carlo (MCMC) method can be adopted, which would yield an estimate of the Bayes factor directly [122], [145], [146]).

For the reason mentioned above, in this work we adopt the BayesWave pipeline [122], [123], which utilizes Bayesian inference to sample a dynamical parameter space. In practice, signal and noise are modeled by an ensemble of Morlet–Gabor wavelet functions by employing a Reversible-Jump (RJ) MCMC [147] algorithm, where the optimal number of wavelets and their corresponding parameters are estimated from the data. In fact, the more complex the data, the larger the number of wavelets is, as in Fig. 2. By “optimal”, we refer to the statistically most probable model that can sufficiently describe the observations. Essentially, Occam’s razor is applied on the model complexity, given the data [62], [148]. The same approach is used to separate the noise (together with possible non-stationarities, such as glitches [149]) and the gravitational-wave signal components.

As already mentioned, BayesWave uses a more flexible parameterization for the signal \(h\). The signal is modeled as a sum of functions \(\Psi_i(\vec{\lambda} )\), each depending on a set of parameters \(\vec{\lambda}\). Then, \(\Psi_i\) can be written as [122], \[\begin{align} \Psi(t;A, f_0, \tau, t_0, \phi_0)=&Ae^{-(t-t_0)^2/\tau^2}\\ &\times\cos[2\pi f_0(t-t_0)+\phi_0], \end{align} \label{eq:wavelet32basis}\tag{4}\] where \(t_0,\, f_0\) are the central time and frequency of the wavelet, respectively. In Eq. (4 ), \(A\) is the wavelet amplitude, \(\phi_0\) a phase offset and \[\tau = \frac{Q}{2\pi f_0}, \label{eq:tau}\tag{5}\] with Q being the quality factor⁵. Thus, the GW signal at the geocenter is represented as \[h_+(t) = \sum_{i=0}^D\Psi_i(t;A, f_0, \tau, t_0, \phi_0), \label{eq:hp}\tag{6}\] where \(D\) is the number of wavelets. An example of such a representation is given in Fig. 3, where the ensemble of wavelets \(\Psi_i\) is to be used in order to obtain the \(h_{+}(t)\) reconstruction.

Continuing in the frequency domain, we write \(\tilde{h}_\times(f)\) as \[\tilde{h}_\times(f) = \epsilon \tilde{h}_+(f) e^{i \pi/2}, \label{eq:hc}\tag{7}\] where \(\epsilon\) is the ellipticity⁶ parameter [122]. For a linearly-polarized wave, \(\epsilon=0\), while \(\epsilon=1, -1\) yields circular polarization. The expression in Eq. (7 ) is a consequence of modeling GW signals as a superposition of the elliptical state, which allows us to encapsulate all possible morphologies of a fully polarized monochromatic wave. By decomposing the signal in terms of spherical harmonics \(_{-2}Y_{lm}\), one can obtain the ellipticity parameter relative to the inclination \(\iota\) of the source as in [150] \[\epsilon = \frac{2\cos\iota}{1+\cos\iota^2}. \label{eq:epsilon32relative32to32iota}\tag{8}\]

Then, the strain in the frequency domain measured by the detector is expressed in the form [151] \[\tilde{h}(f) = \Bigg[\tilde{h}_+(f)F_+ +\tilde{h}_\times(f)F_\times\Bigg] e^{i 2\pi\Delta t(\alpha, \delta)+i\phi}, \label{eq:hf32reconstruction}\tag{9}\] where \(\Delta t(\alpha, \delta)\) is the time-delay relative to the arrival time at the geocenter and \(\phi\) the phase at a fiducial reference time.

To calculate the reconstructed signal using Eq. (9 ), we need to estimate both extrinsic and intrinsic parameters. These are the sky position, orientation, and number of wavelets \((\alpha, \, \sin\delta, \, \psi, \, \cos\iota, \, \phi, \, D)\), and the parameters of each wavelet \((A, \, f_0, \, \tau, \, t_0, \, \phi_0)\), respectively. As mentioned above, these are dynamically sampled with the RJMCMC algorithm. We employ uninformative priors, assuming that \(t_{0}\in \mathcal{U}[t_{c}-0.5,t_{c}+0.5]\), where \(t_{c}\) is the GW trigger time, and \(\phi_{0}\in \mathcal{U}[0,2\pi]\). For the amplitude, we choose the default prior suggested in [122]. Regarding the quality factor \(Q\), its maximum value is set to \(300\) for all network configurations, except the HF detector. For this detector, we set it to 800 due to its improved sensitivity, which results in a relatively high post-merger SNR and consequently in a higher quality factor. Finally, we employ a prior for the dimensionality analogous to [62], i.e. \(D_{\rm{min}}=2\). Furthermore, we set the number of iterations to \(4\times10^6\), of which half are considered a burn-in period and later discarded. The final chain is additionally thinned, by keeping every 100th sample. Thus, we calculate \(N\)=20000 reconstructions.

As a measure to quantify the quality of the reconstructions, we use the overlap with respect to the injected signal. The overlap between the true injected signal in the frequency domain \(\tilde{h}_s(f)\) and the reconstructed \(\tilde{h}_r(f)\) is obtained via \[\mathcal{O}=\frac{\langle \tilde{h}_s,\tilde{h}_r \rangle}{\sqrt{\langle \tilde{h}_s,\tilde{h}_s \rangle}\sqrt{\langle \tilde{h}_r,\tilde{h}_r \rangle}}, \label{eq:overlap}\tag{10}\] where the \(\boldsymbol{\langle \cdot \, , \cdot \rangle}\) represents the inner product between two real time series. This is defined as [152] \[\langle a,b \rangle=2 \int_{f_\mathrm{low}}^{f_\mathrm{high}}\frac{\tilde{a}(f)\tilde{b}^*(f)+\tilde{a}^*(f)\tilde{b}(f)}{S_n(f)} d f , \label{eq:inner95product}\tag{11}\] with \(S_n(f)\) being the detector’s one-sided noise Power Spectral Density (PSD) and \((f_\mathrm{low} , f_\mathrm{high})\) the bandwidth of the analysis, while the asterisk \((^\ast)\) denotes the complex conjugation. The overlap takes values between -1 and 1  [122], [150]; the closer the quantity to 1, the better the reconstruction, while \(\mathcal{O}\) = -1 means an anti-correlation. Note that the optimal signal-to-noise ratio (SNR) is defined as \[\mathrm{SNR}=\sqrt{\langle \tilde{h}_s,\tilde{h}_s \rangle}. \label{eq:SNR32equation}\tag{12}\]

6 Evaluation method↩︎

The purpose of this work is to assess the detectability of the rotational instabilities in the post-merger GW spectrum. We anticipate that all networks (Table 1) will accurately detect and reconstruct the signal that corresponds to the “merger" and”early" post-merger phases (up to \(O(10)\)ms after merger, where the main post-merger GW emission takes place), but not to the “late" post-merger signal (i.e. to the emission after \(O(10)\)ms from merger), where rotational instabilities could trigger a re-excitation of GW emission. It is essential to note that the overlap of the entire post-merger signal will be rather close to one in all cases considered here, since the numerator in Eq. (10 ) is dominated by the strong merger and early post-merger part of the waveform. To avoid biased results (in terms of the ability to infer the presence of rotational instabilities), we specifically constrain the overlap calculation to only the instability part of the signal.

First, we used the BayesWave pipeline to analyze the entire post-merger signal, and obtain its extrinsic and intrisic parameters (described in Sec. 5). As an example, Fig. 4a shows the resulting reconstruction of the entire post-merger signal in the frequency domain, for the case of a BNS source with \(M_{\rm{tot}}=3.1~M_\odot\) at 40 Mpc, optimally oriented with respect to the HF detector. The corresponding reconstruction in the time domain is shown in Fig. 4b. Since this is the most sensitive detector we consider and the source is placed at the nearest realistic distance, the SNR of the post-merger phase was 133, resulting in a very accurate reconstruction up to several tens of milliseconds after merger, including the instability part. In the inset of Fig. 4b, we show in more detail a magnification of the beginning of the instability part, where the accuracy of the reconstruction can be clearly seen.

7 Detectability of the rotational instabilities↩︎

Next, we isolated and evaluated the detectability of the instability part of the signal. To do this, we cropped the time series (around 0.4950s as shown in Fig. 4b) and retained only the part of the signal at \(t>0.4950s\), which contains the instability part. We applied a very weak Tukey window (\(a_{\rm{wt}}=0.01\)) and resized the signal onto the initial duration (1 second) and re-sampled them to the desired sampling rate (8192 Hz). The same was applied to the 20000 samples of the posterior reconstructed signal and a new median of the reconstructed instability part was computed and used in evaluating the overlap restricted to the instability part. The detailed results of this investigation are presented in Sect. 7. In Appendix 10, we present additional cases (for the same model, but different detector networks), where we analyzed the entire post-merger signal.

7.1 Reconstruction of the instability part↩︎

Following the procedure described in Sec. 6, we reconstructed the instability part of the post-merger signals for the different cases we considered.

The first network configuration in Table 1, corresponding to possible future upgrades of the HLV network at twice the sensitivity of O5, did not have the required sensitivity for the instability part to be reconstructed, for any of the injections we considered. At this sensitiviy, only the initial post-merger signal can be reconstructed (see Appendix 10).

Next, we consider the CE+ET network configuration in Table 1, for the case of an optimally oriented source at a distance of 40Mpc, with \(M_{tot}=3.0 M_\odot\). The top panel of Fig. 5 shows, in the time domain, the injected signal for the entire post-merger phase, the instability part and its corresponding median reconstructed signal. The corresponding frequency-domain representation is shown in the bottom panel of the same figure. For this case, the instability part has a small amplitude and the median reconstructed signal has a very small overlap with the injected data. The case of non-optimal orientation is even less promising. Hence, for this configuration, one would not be able to discern the presence of a rotational instability.

For the same CE+ET network configuration, we also consider the case of a merger with \(M_{\rm{tot}}=3.1 M_\odot\). In this case, we do not assume optimal orientation and the distance is 40 Mpc. The top panel of Fig. 6 shows that the instability part is reconstructed in the time domain, although the median waveform has a larger amplitude than the injected data. In the bottom panel of Fig. 6, which displays the corresponding picture in the frequency domain, one can observe that the 90% confidence interval (CI) region is rather broad, but the instability part is still well within this region and with a somewhat lower amplitude than the median reconstructed Fourier transform \(\tilde{h}(f)\). We conclude that, for this case, one would be able to infer the presence of the rotational instability, but the parameter estimation would not be very accurate.

Finally, we consider the HF detector design in Table 1. The two panels of Fig. 7a show that for an optimally-oriented source with \(M_{\rm{tot}}=3.0 M_\odot\) at a distance of 40Mpc, the instability part can be reconstructed fairly accurately in the time and frequency domains. Increasing the distance to 80Mpc (two panels of Fig. 7b) results in a similar median reconstructed signal, but with a broader 90% CI that would not allow an accurate parameter estimation. For both distances of 40Mpc and 80Mpc, the reconstruction reveals the presence of two frequencies with a separation of about 150Hz, which explains the apparent “beating” in the time domain representation. One frequency has a value comparable to the dominant \(f_{\rm peak}\) frequency in the early post-merger phase, whereas the other frequency is somewhat smaller. The two distinct frequencies have similar amplitudes in the frequency domain, indicating that they both persist throughout the duration of the signal in the instability part. Hence, they could be due to two different unstable modes operating at the same time. This remains to be confirmed with a detailed mode analysis.

For the case of \(M_{\rm{tot}}=3.1 M_\odot\) the HF detector achieves a very accurate reconstruction, when the source is optimally-oriented at 40Mpc, as shown in the two panels of Fig. 8a. For this case, only a single frequency (comparable to the dominant \(f_{\rm peak}\) frequency in the early post-merger phase) is active in the instability part and it is recovered with a very narrow 90% CI in the frequency domain, which would allow for accurate parameter estimation. When the same source is set at a significantly larger distance of 200Mpc, the reconstruction is still fairly accurate, but with somewhat broader 90% CI.

7.2 Overlap calculations for the instability part↩︎

Fig. 9a shows the calculated overlap between the median reconstructed signal at the CE detector and the injected signal for the instability part as a function of the SNR of this part. The model in this case is the \(M_{\rm }=3.0M_\odot\) model and we varied the distance to the source to achieve different SNR values. For this model, the instability part is very weak and with the CE+ET network only a very small overlap of \({\cal O} \simeq 0.36\) is achieved. For the same model, the HF detector design performs significantly better, as one can reach \({\rm SNR}\sim 21\) for the instability part of the signal, with an overlap of \({\cal O} \simeq 0.977\).

In 9b the corresponding overlap for the case of \(M_{\rm }=3.1M_\odot\) is shown. This model has a significantly stronger instability signal, allowing the CE+ET network to reach an \({\rm SNR}\sim 8\) with an overlap of \({\cal O} \simeq 0.96\). On the other hand, the HF detector reaches an \({\rm SNR}\sim 44\) with an overlap of \({\cal O} \simeq 0.99\).

To allow a direct comparison between the \(M_{\rm }=3.0M_\odot\) case and the \(M_{\rm }=3.1M_\odot\) case, we display the overlap obtained at the CE detector and with the HF detector design in Figs. 10, respectively.

8 Summary and Conclusions↩︎

Several studies on the dynamical evolution of differentially rotating neutron star configurations suggest a prospective excitation of a low-\(|T/W|\) rotational instability, which could play an important role in the long-term evolution of a BNS post-merger remnant [108]–[110]. This phenomenon causes a dynamical shear instability [109], [110], while stars with a large entropy could be stabilized against it [153]. Evidence for a re-excitation of the of the \(m=2\) \(f\)-mode can be seen in the numerical simulations of e.g. [118]–[120].

We explored the potential for detecting such rotational instabilities if they occur in BNS merger remnants. To do that, we inject the numerically simulated signals of equal-mass BNS mergers produced by [120], with MPA1 EoS and total mass of 3.0 and 3.1 \(M_\odot\), into different network configurations (see Table 1) assuming colored Gaussian noise. We considered that the merge occurs at realistic distances, that is, 40-200 Mpc, and examined different source orientations (see Table 2). Finally, we employed BayesWave to compute N=20000 reconstructed signals, and then we calculated the overlap between the injected waveform and the median reconstructed of the instability part of the signal.

Our study indicates that a network of advanced detectors upgraded to twice the sensitivity of the A + design would only be able to detect the “early" post-merger signal at the closest distance of 40 Mpc, but it would miss the presence of rotational instabilities. We find that for a BNS with \(\rm{M_{tot}}=3.1\) \(M_\odot\) using the MPA1 EoS, the network comprised of the proposed CE and ET detectors will be able to detect rotational instabilities if the event takes place at distances less than 80 Mpc. The HF design would detect rotational instabilities for sources within \(\sim O\)(200Mpc).

For the case of a BNS merger with MPA1 EoS and a lower mass of \(\rm{M_{tot}}=3.0\) \(M_\odot\), all considerd detector designs reconstruct the early post-merger signal of the BNS merger, but not necessarily the signal regarding rotational instabilities. In this case, the spectrum of the waveform reveals the presence of two frequencies, the dominant mode and a second, nearby frequency, which have comparable amplitudes in the frequency domain. For this model, the CE+ET network would not allow for an accurate parameter estimation of the instability part, but the HF would have acceptable accuracy, especially at a distance of 40Mpc.

In this first study, we have considered a simple setup for the BNS simulations, where only hydrodynamics is taken into account and the effect of magnetic fields, bulk viscosity or other viscous effect are ignored. It will be important to perform detailed studies taking into account all physical effects, to determine which ones allow or prohibit the development of rotational instabilities. Future observational studies will either confirm the presence of rotational instabilities or set upper limits and constrain physical properties that contribute to their suppression.

We are grateful to Miquel Miravet for carefully reviewing this manuscript and for the fruitful conversation about BayesWave injections. We also thank Katerina Chatziioannou for useful comments and clarifications on the paper in Ref. [62] and Meg Millhouse, Sophie Hourihane for more information about the BayesWave pipeline. We acknowledge Teng Zhang, Huan Yang, Denis Martynov, Patricia Schmidt and Haixing Miao for providing us a Python notebook to reproduce the sensitivity curves in Ref. [71]. AS acknowledges the Bodossaki Foundation for support in the form of a Ph.D. scholarship. Figures in this manuscript were produced using Matplotlib [154]. The package PyCBC [155] were used for benchmarking of different waveform models, for transformations between the frequency and time domain, and for other results. We are grateful for the computational resources provided by the LIGO Laboratory and supported by the U.S. National Science Foundation Awards PHY-0757058 and PHY-0823459. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan.

9 Injected waveform SNR for different cases↩︎

The spectrum of the waveforms (see Table 2) and the Amplitude Spectral Density (ASD) of each detector sensitivity (Table 1), plus the sensitivity curves for O4 and O5) are presented in Fig. 11. The waveforms correspond to a BNS (MPA1 EoS) with a total mass of 3.1\(M_\odot\) (top panel) and 3.0\(M_\odot\) (bottom panel). The event is considered at 40 Mpc and is optimally oriented with respect to the detector. Labeled SNRs correspond to the optimal SNR of the entire post-merger phase.

For instance, the optimal post-merger SNR for the MPA1 EoS with \(M_{\rm{tot}}\) = 3.1 \(M_\odot\) is computed via \[\left(\frac{S}{N}\right)^2=4 \int_{2000}^{4096} d f \frac{|\tilde{h}(f)|^2}{S_n(f)} \label{eq:optimal95post95merger95snr95example95155}\tag{13}\] where the limits in the integral are in Hz.

From Fig. 11, one can indicate that there is a poor possibility of obtaining valuable information about shear instability during O4 (SNR = 2) and O5 (SNR = 4). Only with a sensitivity twice the sensitivity of A + or better could an SNR of at least 8 be obtained, which is necessary for the detection of sigle events (see, e.g.,  [64]).

10 Analysis for the full signal↩︎

In this Appendix, we consider the same analysis method mentioned in Sect. 5, but without cropping the signal and extracting only the instability part; instead, we perform the analysis in the entire post-merger phase (full injected signal). Figs. 12 and 13 represent the case of a BNS merger with the MPA1 EoS and \(M_{\rm{tot}}=3.1~ M_\odot\) at 40 Mpc. The network cases correspond to (a) “\(2~\times~\rm{O}5\)" and (b)”CE+ET" (see Table 1) and the figures refer to the L1 and CE channels, respectively. Figures show the injected waveform (gray), the median reconstructed waveform (pink), and its 90% CI in the frequency domain (Fig. 12) and time domain (Fig. 13). Additionally, in Fig. 12, the analyzed data (light gray), which is labeled as “Noise+injection", are shown. The noise is Gaussian and is generated with respect to the ASD sensitivity (black dashed line). Both figures indicate the capability of all detectors to recover the post-merger GW signal.

For the case of the “\(2~\times~\rm{O}5\)" network configuration, Fig. 12 reveals the capacity of the network to capture the peak frequency, while the "early" post-merger is well reconstructed. The performance of the”\(\rm{CE+ET}\)" network seems very good.

For a quantitative evaluation, we calculate the overlap \(\mathcal{O}\) between the entire injected signal and the median reconstructed signal as a function of the optimal post-merger SNR. The overlaps corresponding to the “\(2~\times~\rm{O}5\)" LIGO detector are good, despite the fact that we note difficulties in recovering the shear instability. The overlap is larger due to the good reconstruction of the early post-merger signal. In both mass models, the CE detector reaches \({\rm SNR}\sim 48\) with an overlap of \({\cal O} \simeq 0.99\), while the HF design reaches \({\rm SNR}\sim 135\) with an overlap of \({\cal O} \simeq 0.999\).

The median recovered peak frequency \(f_{\rm{peak}}\) (in Hz) for the MPA1 EoS and for both mass cases, as a function of the optimal post-merger SNR is shown in Fig. 15. For the case of a BNS merger with a total mass of 3.0 \(M_\odot\) (Fig. 15a), the recovered \(f_{\rm{peak}}\) reaches an uncertainty of \(\Delta f_{\rm{peak}}\sim 35\) Hz for the HF design. The “CE+ET" network configuration seems to be able to give similar results for an optimal post-merger network SNR 37. For lower SNRs, the data become more uninformative and the error in the peak frequency becomes higher.

When the total mass of the BNS is 3.1 \(M_\odot\) (Fig. 15 b), the HF detector captures almost perfectly the peak frequency (\(\Delta f_{\rm{peak}}\lesssim 3.4\)Hz), resulting in good estimates up to 200 Mpc (post-merger SNR \(\sim 34\)). Regarding the performance of the “CE+ET" network, the peak frequency is well estimated when the source is optimally oriented at 40 Mpc (post-merger network SNR \(\sim37\)).

References↩︎