Tidal reconstruction of neutron star mergers from their late inspiral
Souradeep Pal K Rajesh Nayak


0000-0003-2172-8589]Souradeep Pal
Indian Institute of Science Education And Research Kolkata 

Mohanpur, Nadia - 741 246 

West Bengal, India
0000-0002-6814-7792]K Rajesh Nayak
Indian Institute of Science Education And Research Kolkata 

Mohanpur, Nadia - 741 246 

West Bengal, India

We investigate the measurement correlation between the effective spin and the effective tidal deformability in gravitational wave signals from binary neutron star mergers. We exploit the fact that the tidal effects in a binary system are prominent when the components are closer during the late-inspiral. Thus, we indicate to a computationally efficient strategy of extracting the tidal information compressed within seconds before the merger. We report our observations for  and explore the suitability of our approach for the upcoming observation scenarios. Fast and accurate measurements of the tidal deformability parameters can be used to inform astronomers in prioritizing the electromagnetic follow-up efforts for such sources.




§ INTRODUCTION
 
Even decades after the first observational evidence for the existence of  neutron stars, the pursuit for a consistent model for their equation of state continues in the broad astrophysics community. The equation of state (EoS), which describes how the matter inside the stars resists collapse against their own gravitational pull, can be constrained by their masses and the radii which are inferred through observations. Recent detections of the neutron star mergers by the current generation of gravitational wave observatories have opened up a new window of investigations. 

The first three observing runs (O1-O3) of the advanced LIGO and the Advanced Virgo have reported two significant detections of binary neutron star (BNS) mergers <cit.>. The joint detection and the association of gravitational waves (GW) and electromagnetic (EM) emissions for the event  have put multi-messenger astronomy into the limelight <cit.>. However, till now   is the only reported multi-messenger GW event observed by multiple observatories. A couple of neutron star - black hole (NSBH) systems have also been detected <cit.> but no EM counterparts were reported. 

Accurate measurement of the tidal deformability in a neutron star caused by its companion object is singularly helpful in constraining the radius and the EoS of the star. Additionally, these measurements can be used to probe the black-hole-nature of the components and provide valuable insights on the EM-properties, crucial for decision making by the astronomers during EM follow-ups. 

Several techniques, with or without involving EM observations, have been used to measure the tidal effects primarily from the inspiral part of the GW signal <cit.>. While estimating the binary's parameters, Bayesian methods are generally robust, however, in the case of tidal deformability parameters, the inferences can be prone to variations from the choice of the prior model <cit.>. Nevertheless, full Bayesian analyses on fairly large populations of simulated BNSs are computationally challenging, given the long duration of the signals. Packages like  offer fast schemes for inferring the binaries' properties, but can still take hours to process a single BNS signal even on a medium-sized computing cluster <cit.>. Typically, a BNS signal spends hundreds of cycles in a detector's sensitive frequency band, however, the last several cycles to the merger would carry the best imprints of the tidal deformation effects. These interactions in a binary system are expected to be prominent when the components are in close proximity of each other before the merger. This is indicated in the literature <cit.>. However, the standard approach is to match signals starting at lowest possible frequencies (flower∼ 20 Hz) by using the data containing the complete inspiral <cit.>. 

In this work, we show that in the full frequency band of the detectors, spins can interfere with the accurate measurement of the tidal deformability parameters. Thus, we focus on the late-inspiral part of the signals by starting the analyses at a higher frequency (flower∼ 150 Hz), using the data available just seconds before the merger. Using simulations, we show that this approach can efficiently provide accurate estimates of the effective tidal parameter, λeff in the presence of the component spins. Here we consider the component spins that are aligned to the orbital angular momentum of the binary system. Note that λeff is better estimated than the individual tidal deformability parameters in a similar sense that the mass-weighted effective spin parameter, χeff is more accurately measured than the individual aligned-spins. These two quantities represent effects that are not explicitly correleated and can be independently calculated using

    χeff=m_1χ_1z+m_2χ_2z/m_2+m_2 ,


    λeff=16/13m_1^4(m_1+12m_2)λ_1+m_2^4(m_2+12m_1)λ_2/(m_2+m_2)^5
where m_i, χ_i and λ_i represent the individual masses, the dimensionless spin parameters and the dimensionless tidal deformability parameters respectively.

There are several techniques to blend these effects together into waveform models. The classical approach using the post-Newtonian (PN) approximation is mostly valid in the low frequencies of the signal when the separation of the components is large [ terminates at the frequency corresponding to the innermost stable circular orbit <cit.>. All the waveform models used in this work are publicly available and are implemented in  <cit.>.]. Using the effective-one-body (EOB) formalism, it is possible to include the tidal corrections till the initial contact of the two bodies <cit.>. In general, the coalescence of a neutron star merger becomes increasingly difficult to model as the components reach the merger, a regime where the numerical relativity (NR) informed approaches are expected to be more accurate <cit.>. It is also possible to incorporate the NR tidal corrections to precessing spin waveform models. However, parameter estimations using these waveform models have inferred only nominal spins for the GW events reported so far <cit.>.

The work is organized as follows. In Section <ref>, we discuss an overview of the GW searches for BNS systems involving spin as well as the tidal effects. We use particle swarm optimization (PSO) for recovering signals from the data. PSO is a stochastic algorithm that iteratively optimizes a given fitness function in an arbitrary dimensional parameter space. We demonstrate a novel method to efficiently decouple the effective spin and the effective tidal deformability. We primarily investigate the match computed between the data and model waveforms starting at different low frequency cutoffs. The inner product between two vectors a and b is given by

    ⟨a|b⟩=4∫_flower^fupperã(f)b̃^*(f)/Sn(f)df .
Here Sn(f) represents the power spectral density (PSD) of the data which downweights the contribution of noise in the desired frequency domain [flower, fupper]. In Section <ref>, we report our observations for  using the new approach and discuss a prospective strategy for such analyses in the upcoming observing scenarios. We make a brief remark on the computational cost of the analyses conducted with the new approach.


§ METHODS




 §.§ Searches with the tidal parameters

Gravitational wave searches constitute the first steps of astrophysical data processing whenever the analyzable data arrive from the observatories after calibration. The broad goal of the searches is to essentially identify events of astrophyical origin from the strain timeseries data to enable detailed follow-up investigations on them. The candidate events need to be consistent with several astrophysical assumptions such as that the triggers at different detectors caused by astrophysical events cannot be separated in time by more than the light's travel time between the detectors. Triggers are computed from the matched-filter SNR timeseries, ρ(t), realizing Eq. <ref> with the strain data h and a template waveform q

    ρ^2(t;θ)=4∫_flower^fupperh̃(f)q̃^*(f;θ)/Sn(f)e2πiftdf .

The data from the observatories are cross-correlated with banks of template waveforms that include various physical parameters (θ) describing the binaries. The number of the such parameters is often regarded as the search space dimensionality. Specific search configurations are adequately described in the literature <cit.>. Searches for binary neutron stars in the Advanced LIGO - Advanced Virgo era have generally included small aligned-spins in addition to the component masses for the templates <cit.>. The first template bank search for neutron star mergers with the tidal parameters was only conducted recently by <cit.>. However, their templates did not include any spin and resulted into no additional significant detection. Instead the search reported even smaller SNRs for some of the events already detected by aligned-spin searches. In general, including larger parts of the existing parameter space or extending the dimensionality of the existing parameter space is computationally challenging and can also result into higher levels of background events. We test the detection algorithm based on PSO to include the tidal parameters in the presence of aligned-spins. PSO allows us to explore a higher-dimensional search space without the need of designing and validating a template bank while requiring a relatively smaller number of template point realizations. The algorithm finally returns the optimal detection statistic along with a set of the point estimates for the source parameters in a frequenctist way <cit.> [Earlier, a PSO based approach has been reported that explores the orbital eccentricity in the presence of moderate aligned-spins]. However, here we do not explicitly conduct a search targeting real events over the archived LIGO-Virgo-KAGRA (LVK) data and conclude based on simulations only. For this, we inject about 1000 simulated signals starting at 20 Hz into simulated data segments having O4 model noise PSDs. The component masses and the tidal parameters are randomly chosen from [1M_⊙,3M_⊙] and [0, 4000] respectively having a uniform distribution. These injections carry precessing spins with each spin component having a normal distribution around zero with the spin magnitudes restricted upto 0.95.

The signals are recovered with different combinations and/or ranges of template parameters using the PSO-based search algorithm. We use the  model to generate templates choosing from a subset of parameters given in Eq. <ref>. A combination includes at least the individual aligned-spins or the tidal parameters as the search parameters in addition to the component masses. The search space ranges are also allowed to vary for the aligned-spins. The maximum search dimensionality explored is 6 which includes both the aligned-spin and the tidal parameters.

    θ ={m_1,m_2,χ_1z,χ_2z,λ_1,λ_2} .

Within the PSO framework, the following two factors independently affect the recovery of the optimal match between the data and a template- (a) the waveform model, which ensures that the template is equipped with the physics to adequately model the incoming signal and (b) the template sampling, which ensures that a sufficient matched-filter computations are performed throughout the search space in order to locate the peak value of the match function. For computational efficiency, a faster waveform model is preferred. On the other hand, the sampling is preset in a way that just recovers signals at their optimal statistic as demonstrated in Fig. <ref>. In this work, we have tuned the search so that for a given stretch of data, around 25500 template points are computed. These many points in the search parameter space are used to generate the template waveforms which are then used in the matched-filtering of the data from each detector. The resulting triggers from the individual detectors are processed with several consistency tests to obtain the coincident events. The events are assigned a detection statistic from their reweighted SNRs and a false alarm rate based on time-slides <cit.>. The methods adapted for the PSO algorithm are described in <cit.>.

A summary of the simulation results is the following. Between spins and tides, the former play the dominant role in obtaining the optimal detection statistic. We find that the templates used in searches should have at least large aligned-spins to recover sources that may have generic spins. On the other hand, the tidal parameters appears to have very small contribution to the build-up of the detection statistic. A search with templates having tidal parameters but no spin can lose a significant fraction of the retrievable statistic. We also exercise these searches on real datasets from O3b that do not contain any reported GW signal or known hardware injection <cit.>. These datasets are 4096 seconds in duration and the searches use a similar setup as described earlier. We observe that all the template configurations respond to the noise events from the detector data in a similar fashion, resulting into equivalent rates of non-astrophysical background events, as summarized in Fig. <ref>. Thus we conclude that the search parameter space should include the tidal parameters in the presence of large generic aligned-spins to approach the optimal sensitivity towards the observable binary neutron stars.








 §.§ Estimation of the effective tidal parameter

In this section, we discuss the effect of the measurement correlation between the effective spin and the effective tide using the likelihood function. The likelihood (ℒ) of a signal with the parameters θ in the data d assuming a Gaussian noise model is given as

    ℒ=e^-⟨d-h(θ)|d-h(θ)⟩/2 ,

where ⟨a|b⟩ represents the inner product of two vectors a and b defined in Eq. <ref>. The likelihood function is extensively used in the Bayesian approach to generate the posterior probability distributions of a desirable subset of the parameters. Several of these parameters can be degenerate with one another- the measurement of one can interfere with the measurement of another. The log-likelihood ratio (LLR) is often used to quantify a competing signal hypothesis against the noise hypothesis which simplifies to,

    LLR=⟨d|h(θ)⟩-1/2⟨h(θ)|h(θ)⟩ .

To understand the measurement correlation between the effective spin χeff and the effective tidal parameter λeff, we sample the LLR as a function of λeff and χeff for a random BNS injection while keeping other signal parameters fixed at their injected values. The injected signal is for a simulated source with component masses (1.4M_⊙, 1.4M_⊙) at 200 Mpc generated using the  waveform model. The signal is projected onto a simulated HL detector network (O3). For sampling the likelihood (or the LLR), we choose moderate values of χi upto ±0.5 and λi upto 4000 to generate an equally spaced grid of points in the χeff-λeff plane. These points can be easily obtained by varying the values of χi and λi which are preserved thoughout the simulation. Finally the LLR is computed at all the grid points and the contours containing the top likelihood points are plotted. The likelihood computations are repeated by changing the value of flower, as shown in Fig. <ref>. We observe that the contours include a broad range of values for λeff when the waveforms are computed from a small flower∼ 23 Hz whereas it peaks sharply near the injected χeff value (blue). This trend keeps reversing on increasing the flower with the likelihood function becoming more contrained towards λeff (pink), in spite of the fact that the likelihood value diminishes in this process due to the loss of the signal strength in the relevant frequency band.






In the above simulation, we assumed the masses to be fixed at the injected values for computing the likelihood, however in real cases, the component masses are unknown and thus need to be estimated as free parameters. In fact in the most general case, all the signal parameters need to be treated as free parameters and be searched over their possible domains. However, in the following study we just free the component masses to explore the resulting six-dimensional space, since for the current scope, we are interested in obtaining intrinsic parameters accurately, specifically the tidal parameters, as summarized in Fig. <ref>. We keep other extrinsic parameters like the coalescence time, sky location, distance, etc. fixed at their injected values. In practice, this can also be tackled with a suitable marginalization scheme or by roughly treating some parameters as semi-free parameters that are available from external observations. For example, in the case of , we can combine the information from the EM observations and use the known sky position to sample over a smaller window for optimization. On the other hand, the coalescence time of the event can be sampled around the trigger time of the GW event which is a standard practice in GW parameter estimation. Since computing the likelihood function on such a high-dimensional space with a regular grid is computationally expensive and also not optimal, we use PSO with the LLR as the fitness function. The injection distribution of this simulation is similar to that mentioned in the previous section except that these injections are non-precessing with uniform aligned-spin components generated and recovered using the  model.

Fig. <ref> shows the efficiency in recovering the LLR (or the SNR) as compared to the true LLR (or the true SNR) for multiple random injections. This is obtained by releasing the extrinsic parameters as semi-free parameters as shown in Table. <ref>. The PSO algorithm stochastically locates the optimal solution, realizing the maximum likelihood or match through a point estimate in the search parameter space. The parameter values of the optimally matching template are prone to statistical variations from several sources, for example from a random fluctuation due to a specific noise realization or the choice of the seed for initializing the algorithm. Note that with a higher flower, the computational cost of the analysis becomes reasonably small mainly due to (a) the smaller duration of the data that is analyzed, ∼4 seconds before the merger and (b) the faster generation of template waveforms due to the smaller number of inspiral cycles computed. We find that these analyses take seconds (or minutes) on multi-core (or single-core) machines. The reduced computational cost also makes it possible to conduct multiple rounds of the same analysis, each with a different seed to initialize the PSO algorithm. Thus, begining with a distinct seed, the algorithm chooses a distinct path, eventually converging on the contour of the optimal likelihood, which we call a multi-swarm exploration. This is demonstrated in Fig. <ref> for a -like injection. We infer the source parameters from the resulting distribution of the point estimates thus obtained as described in the next section.




§ OBSERVING SCENARIOS



 §.§ Current observations





In the previous section, we discussed the PSO algorithm in the 6D intrinsic parameter space. Here we are interested in the BNS signals from the Advanced LIGO and the Advanced Virgo observing runs, whose extrinsic parameters like the sky location, the luminosity distance and the inclination angles are approximately known from observations. The BNS events reported so far are  and . The latter was observed only by the LIGO Livingston observatory and no electromagnetic counterpart could be observed. On the other hand, the former event was recored at both the LIGO Hanford and the LIGO Livingston observatories with a coincident network SNR > 32 in the ∼ 23 Hz - 1024 Hz band, making it the loudest GW event till the end of O3. The HL coincident SNR of ∼ 16 in the ∼ 150 Hz - 1024 Hz band apparently makes it a suitable candidate for our analysis. A sub-threshold signal was also observed in Virgo improving the sky localization of the source. The sky location and the distance information for the source are also approximately known from other astronomical (EM) observations <cit.>. We conduct the multi-swarm analyses with the GW data ∼ 4 seconds close to the merger time of the event with the goal of getting more insights on its tidal effects. 

While we expect that the estimation of Mchirp and χeff are more accurate in the full frequency range, this analyses still provides comparable estimates with the previously reported results. The component mass estimates are consistent with the Mchirp-uncertainty contour, which we expect to be slightly broadened due to the smaller signal strength in the chosen analysis band, whereas the estimated χeff is fairly consistent with a nominal value as shown in Fig. <ref>. As demonstrated in the earlier section, λeff which is expected to be better estimated with a higher flower, is significantly different from that of the previously reported results <cit.>. However, note that such an estimate is prone to variations arising from the random noise fluctuations and the inference can be dependent on the specific noise realization of the instruments at the time of recording of the signal. The underlying variations are indicated for scenarios in O3 and beyond in Fig. <ref>. Thus, we cannot establish the true estimate for λeff with certainty based on the above observation, primarily due to the effect of the noise realization. Note that in the full Bayesian setting, a 90% credible region could exceedingly miss the true values at the current or past sensitivities as discussed in <cit.>. Either multiple such observations are needed to reduce the effect of noise realizations from the detectors or exceptionally loud events can possibly inform us about the true tidal deformability parameters in the future. Thus, we do not explore further implications on the EoS of neutron stars based on the above inferences. However, we can conclude that a less stiff EoS may not be strictly ruled out based on the few observed neutron star binaries. This may be in contradiction to the general preference for a stiff EoS in the community.


 §.§ Future observations


In the next generation of ground-based GW detectors, the rate of detections is likely to increase depending on the astrophysical abundancy of the sources and the factors like the sensitivity and the uptime of the instruments. The total SNR of the majority of the detectable events could be significantly large given the number of low-frequency inspiral cycles that the plausible flower would contain. Roughly, a signal from 1 M_⊙ - 1 M_⊙ would last for more than 10,000 seconds in a detector's sensitive frequency band starting at, say 5 Hz. Data analysis, including the standard Bayesian parameter estimation, of such long signals is a general challenge for GW astronomers. However, with the proposed sensitivities of the Cosmic Explorer (CE) and the Eintein Telescope (ET), a GW170817-like event could still register SNRs greater than 100, even with a higher flower as shown in Fig. <ref>. We find that in such scenarios [including LIGO Voyager (VO)], the effects of the random fluctuations from different noise realizations would be negligible and can thus, lead to reliable inferences from the estimated tidal parameters. For these future observations, the tidal properties can be efficiently extracted by exploring only the higher frequencies, despite the improvement in the low frequency sensitivities of the detectors like CE and ET. We expect the spin-tidal measurement correlation to become more pronounced for the full frequency band analyses with the improved low frequency sensitivities. Thus the method presented here can be promising, however, to include the related effects at further higher frequencies, a higher sampling rate would be considered (e.g. to increase fupper from 1024 Hz to 2048 Hz, we just need to increase the sampling rate from 2048 Hz to 4096 Hz [which would roughly be just twice as costly as presented here; the overall cost of the analyses would still remain low since a small duration of data is used]). Note that the Mchirp information for low mass binaries is a fairly accurate byproduct of the detection algorithms. The technique proposed here could possibly lead to data products that provide valuable insights on the EM properties of the sources and can also be used to provide source classification in low-latencies.


§ CONCLUSIONS

This study shows that ignoring large generic aligned-spins in the presence/absence of tides in the detection stage leads to sub-optimal recovery of neutron star signals. However, the measurement of the effective spin parameter of the neutron star components can affect the measurement of the effective tidal parameter in the full frequency bandwith of the detectors. We have discussed a novel technique to extract the tidal information primarily from the late-inspiral part of neutron star mergers signals observable with the ground-based detectors. Though a restricted frequency band reduces the recovered SNR, we find the strategy effective in probing the tidal deformability effects in binary systems. We have also explored an empirical low frequency cutoff, flower for the signals  that minimizes the errors in the point estimates of the effective tidal parameter, λeff using particle swarm optimization (PSO). We observed that a frequency band around 150 Hz - 1024 Hz efficiently provides more reliable estimates for a test population of simulated binary neutron stars than the full frequency band. We also found that this technique benifits from the numerical relativity (NR) informed waveform models which excel in modeling the near-coalescence regime. The study can be extended further to understand the correlation with/among other effects such as eccentricity and spin-precession. 

For -like events, we conclude that with this method, the λeff estimates are expected to become free from uncertainties due to the detectors' noise realizations in the future observing plans. The new approach could be especially effective in the upcoming observing scenarios when the duration of these signals would be longer and hence, computationally challenging with the improved/proposed low frequency sensitivities of the current and the upcoming ground-based detectors. Further investigations with the above strategy are desirable for binaries containing sub-solar mass objects or black holes accompanying the neutron stars. 

In the case of neutron star-black hole (NSBH) systems, we observe that the tidal deformability of the NS tends to depend on the mass of its BH counterpart, as shown in Fig. <ref>. It is believed that the tidal deformability of a BH is zero which can be used as a test of the BH-nature of a component. However, the BH can cause measurable tidal deformability to the accompanying NS under some circumstances. As demonstrated in this work, the tidal component is almost uninforming at flower∼ 20 Hz and starts constraining at a higher flower. Nevertheless, we find that even at a higher flower, the contours tend to be unconstrained for large mass ratios (q) systems, possibly indicating the absence of a measurable tidal deformability. Because a large mass BH would merge quickly with a typical NS, thus allowing almost no time to tidally deform and disrupt the NS. Generally, the tidal interactions can be expected when the BH mass is comparable to the NS mass, while also increasing the chances for an EM emission. We also encourage the readers to refer to an interesting discussion in the following review article <cit.>.


§ ACKNOWLEDGMENTS

This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. 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.

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 article has a LIGO Document number LIGO-P2400130. We have primarily used the PyCBC software in conducting various analyses in this work. We acknowledge the use of the Sarathi computing cluster at IUCAA for computational/numerical work. A significant part of the numerical work was also carried out at the Center of Excellence in Space Sciences, India (CESSI). CESSI is a multi-institutional Center of Excellence hosted by the Indian Institute of Science Education and Research (IISER) Kolkata and has been established through funding from the Ministry of Education, Government of India. SP wants to thank the financial support from the Council for Scientific and Industrial Research (CSIR), India through File No:09/921(0272)/2019-EMR-I.



aasjournal
