Improved Upper Limits on Gravitational-wave Emission from NS 1987A in SNR 1987A

We report on a new search for continuous gravitational waves from NS 1987A, the neutron star born in SN 1987A, using open data from Advanced LIGO and Virgo ’ s third observing run ( O3 ) . The search covered frequencies from 35 – 1050 Hz, more than 5 times the band of the only previous gravitational-wave search to constrain NS 1987A. Our search used an improved code and coherently integrated from 5.10 to 14.85 days depending on frequency. No astrophysical signals were detected. By expanding the frequency range and using O3 data, this search improved on strain upper limits from the previous search and was sensitive at the highest frequencies to ellipticities of 1.6 × 10 − 5 and r -mode amplitudes of 4.4 × 10 − 4 , both an order of magnitude improvement over the previous search and both well within the range of theoretical predictions.

1. INTRODUCTION Piran & Nakamura (1988) first suggested shortly after SN 1987A that the neutron star (NS 1987A) probably born in the supernova could be emitting detectable continuous gravitational waves (GWs).Yet the first search to constrain the behavior of the neutron star through GW upper limits was not performed until recently (Owen et al. 2022).That search used open data from Advanced LIGO and Virgo's second observing run (O2), covered the frequency band 75-275 Hz, and was sensitive in that band to GW signals just below an analog of the pulsar spin-down limit based on the age of the neutron star (Wette et al. 2008).Here we describe a search of open data (Abbott et al. 2023) from Advanced LIGO and Virgo's third observing run (O3) (Tse et al. 2019;Acernese et al. 2019) using a new and improved code covering a wider frequency band (35-1050 Hz).
The detection of neutrinos from SN 1987A (Bionta et al. 1987;Hirata et al. 1987) suggests that a neutron star, rather than a black hole, was the most likely product of this event.Located 51.4 kpc away in the Large Magellanic Cloud (Panagia 1999), NS 1987A is the youngest neutron star in our galactic neighborhood.Electromagnetic searches for a pulsar or non-pulsing neutron star in the remnant SNR 1987A are made difficult by its dust-filled surroundings.Far infrared observations of SNR 1987A by Cigan et al. (2019), however, have detected a relatively warm, compact region of dust that could be powered by a very young cooling neutron star (Page et al. 2020;Dohi et al. 2023).Greco et al. (2021) and Greco et al. (2022) argue that hard X-ray emission suggests the presence of a pulsar wind nebula.
The first searches for continuous GW emission from NS 1987A used stochastic background methods to analyze data from Advanced LIGO and Virgo's observing runs-see Abbott et al. (2021a) and references therein.These searches did not cover a reasonable parameter space.They assumed spin-down rates for NS 1987A of order 10 −9 Hz/s, which is large by the standards of known pulsars but is at least an order of magnitude smaller than the spin-down that would be caused by GW emission at detectable levels (Owen et al. 2022).Stochastic background searches also did not achieve the needed sensitivity to detect GWs at the level of the indirect upper limit on GW strain h age 0 , analogous to the spindown limit for pulsars but based on the age of the object when pulses are not observed.This age-based indirect limit was defined by Wette et al. (2008), who described the basic continuous-wave method for searches for persistent GWs from supernova remnants where there is evidence for a neutron star but where pulses are not observed (such as SNR 1987A).Such searches require that wide bands of frequencies and spin-down parameters (time derivatives of the frequency) be explored.These continuous GW searches use longer signal coherence times than stochastic back-ground searches, and therefore generally require searching over spin-down parameters as well as GW frequencies.The limit h age 0 is a useful figure of merit for search sensitivity.The method of Wette et al. (2008) has been used to search for many objects, starting with the central compact object in supernova remnant Cas A (Abadie et al. 2010) and more recently targeting NS 1987A (Owen et al. 2022).
Until Owen et al. (2022), searches for NS 1987A did not cover the full parameter space or reach the sensitivity of h age 0 .Narrow parameter space searches for continuous GWs from NS 1987A were performed by Sun et al. (2016) using the method described by Chung et al. (2011).Since the Wette et al. (2008) parameter space required a fourth spin-down parameter for NS 1987A (then 19 years old), Chung et al. (2011) narrowed the search by introducing a detailed spin-down model.This is less robust than models making fewer assumptions.Even with a narrow parameter space, Sun et al. (2016)  for mass-quadrupole GW emission and Owen (2010) extended it to currentquadrupole GW emission from r-mode oscillations [see Glampedakis & Gualtieri (2018) for a summary of emission mechanisms].The derivations of these limits assume that GWs dominate the spindown of the star from birth and that the initial spin frequency was much higher than present.The Wette et al. (2008) mass-quadrupole age limit on the GW amplitude h 0 [a measure called the intrinsic strain (Jaranowski et al. 1998)] can be written as a frequency-independent expression that depends on the age of the neutron star a, its distance D, and moment of inertia I: (1) The analogous indirect limit for GW emission from rmodes is given by (Owen 2010) where n = f f / ḟ 2 is the braking index.Here we have modified the expression of Owen (2010) to leave free the parameter A, which is the ratio of r-mode frequency (in an inertial frame) to stellar spin frequency.We convert the above expressions to numerical ranges as follows.Owen (2010) used A = 4/3, appropriate for slow rotation and Newtonian gravity, while general relativistic slow rotation estimates are about 1.39-1.64(Idrisy et al. 2015;Ghosh et al. 2023).The moment of inertia depends on the neutron star equation of state and mass.We take the mass range for neutron stars to be 1.2-2.1 M ⊙ (Martinez et al. 2015;Cromartie et al. 2020).For this mass range and the equations of state used by Abbott et al. (2021b), the moment of inertia range is about 9.1 × 10 44 -4.7 × 10 45 g cm 2 .Finally we use 6-7 for the physically plausible range of r-mode braking indices (Lindblom et al. 1998;Ho & Lai 2000).The resulting range of h age 0 for NS 1987A, using an age of 33 years (applicable for late O3), is about (3) for mass-quadrupole GW emission using Eq. ( 1) and for r-mode GW emission using Eq. ( 2).Inserting these parameters into the results of Wette et al. (2008) and Wette (2012) indicates that a coherent search of O3 data using only two spin-down parameters can surpass the sensitivity of h age 0 for a computing budget of order a million core-hours.This paper describes such a search, which detected no astrophysical signals but placed direct upper limits on the GW strain from NS 1987A.These limits beat the indirect limit h age 0 over a physically consistent parameter space that is considerably larger than the range of frequencies explored in the O2 data search by Owen et al. (2022).

SEARCH METHODS
The search methods used in this paper are similar to those used by Owen et al. (2022).Highlights and changes are summarized here.Readers are directed to Owen et al. (2022) and references therein for details.
The Drill pipeline (Owen et al. 2023) version 1.0.0 was used for this search.It will be described more fully elsewhere.Here we summarize the differences between Drill and previous codes.Drill is a completely new code with functionality similar to that used in Owen et al. (2022).It consists of Python scripts running C codes from LALSuite (LIGO Scientific Collaboration 2020) (v6.25.1 of lalapps and concurrent versions of other packages) that implement the multi-detector Fstatistic (Jaranowski et al. 1998;Cutler & Schutz 2005).Drill is more efficient than the code used in Owen et al. (2022), vetoes signal candidates based on nonparametric statistics, and handles upper limits more consistently.
This search used open data (Vallisneri et al. 2015;Abbott et al. 2023)  Short Fourier Transforms (SFTs) generated from Frameformat time-domain data sampled at 4 kHz.All SFT logic and data selection were handled by the lalapps MakeSFTs program.O3 and our search included data from the Hanford, WA (H1) and Livingston, LA (L1) 4 km LIGO interferometers and the Cascina, Italy (V1) 3 km Virgo interferometer.Although V1 was generally less sensitive, we verified that including its data improved search sensitivity at fixed computational cost.The observation time spans for our search bands were also set to achieve a fixed computational cost.Following Jaranowski et al. (1998) the start time of each span was chosen to maximize the data time divided by the joint power spectral density (psd) of strain noise, which is approximately equivalent to maximizing the search sensitivity.
For NS 1987A we used the (J2000) right ascension and declination α = 05 h 35 m 27 s .998,δ = −69 • 16 ′ 11 ′′ .107(5) from Cigan et al. (2019).We used an age of 33 yr (suitable for the late O3 data used) and the distance 51.4 kpc (Panagia 1999) to determine the parameter space and infer source properties.Observation spans and other search parameters derived from the age are listed in Table 1.The signal parameter space was chosen similarly to Owen et al. (2022).The GW frequency in the solar system barycenter frame was modeled by where t 0 is the time at the beginning of the span and the parameters (f, ḟ , f ) are evaluated at epoch t 0 .The ranges of ( ḟ , f ) for a given value of f were with the braking index n ranging from n min = 3 to n max = 7 and a being the neutron star's age.Unlike in previous searches, a minimum braking index of n min = 3 was used to keep the highest values of − ḟ below about 5 × 10 −7 Hz/s, where Owen et al. (2022) found that the SFT length of 1800 s can become problematic.We verified that the other consistency checks described in Owen et al. (2022) were satisfied.This range of braking indices is consistent with the minimum spin-down for a given h 0 (Owen 2010), and is greater than the maximum value covered by the all-sky surveys such as Dergachev & Papa (2023).It does not include some of the more extreme range recently proposed by Morales & Horowitz (2023), which is more appropriate for older stars.
The frequency band was split into a low frequency band from 35-125 Hz, a medium frequency band from 125-450 Hz, and a high frequency band from 450-1050 Hz.The boundary at 125 Hz was chosen so that even the fastest spinning young pulsar (Marshall et al. 1998) emits GWs in the low band.The boundary at 450 Hz was chosen so that the middle band would avoid noise artifacts due to violin modes of the LIGO test mass suspensions (Covas et al. 2018;Davis et al. 2021).The overall lower bound of 35 Hz was chosen so that an a priori estimate of search sensitivity (Wette 2012) indicated that Eq. ( 1) would be achieved.The overall upper bound of 1050 Hz was chosen for the same reason and because going much higher in frequency causes difficulties with the analysis such as the spin-down range mentioned above.
This search ran a total of roughly 1.2 × 10 6 core-hours on the Texas Tech "Nocona" computing cluster, split into batch jobs of about 8 core-hours each.Integration times and other parameters are shown in Table 1.Each search job covered the full range of spin-down parameters appropriate for its frequency band.Widths of these bands ranged from about 1-37 mHz depending on frequency.Parameter spacings in ḟ and f were of order 10 −11 Hz/s and 10 −18 Hz/s 2 respectively.Search jobs contained about 4-9 × 10 8 templates each, for totals of about 2, 4, and 6 × 10 13 templates for the low, medium, and high frequency bands respectively.
The approach of Drill to vetoing signal candidates does not rely on time-frequency behavior or known instrumental lines for a priori vetoes, as in Owen et al. (2022).Instead it outputs a 2F histogram for each search job without recording any specific candidates on the first pass.This ameliorates major storage and input/output issues in noisy bands.From each histogram Drill takes the loudest 2F (to within the binning resolution of 0.1) and computes (approximately due to the binning) the discrete Cramér-von Mises statistic (Choulakian et al. 1994).The continuous version of the Cramér-von Mises statistic can be written where C and C * are the observed and expected cumulative distributions respectively.For the F-statistic, the latter is a χ 2 with four degrees of freedom.Thus ω emphasizes the middle of the distribution (where non-Gaussian noise is concentrated) more than the tail (where a detectable signal will be) and is a good way of checking for "bad" noise bands without discarding loud signals.
We determined the threshold of ω empirically.Even for Gaussian stationary noise, the output of the LALSuite F-statistic code is not precisely χ 2 (4) distributed due to approximations used to increase computational speed.We checked the behavior of ω in real noise and simulated Gaussian noise, with and without injected signals.In Gaussian or nearly Gaussian noise, we found that ω has a mean of 0.0035 and a standard deviation of 0.0012.Therefore ω = 0.02 corresponds to about thirteen standard deviations and should not veto clean noise bands even with our large number of trials.Loud enough signals should cause a high ω due to many templates triggering at high 2F with slightly wrong parameters.Through injection studies we found that signals are not spuriously vetoed until 2F ∼ 10 6 , which is orders of magnitude above the physical limit h age 0 .

SEARCH RESULTS
We checked the search results for candidate signals as follows.First we vetoed the entirety of each search job which produced ω ≥ 0.02, totaling about 20 Hz or 2% of the total search band.We calculated the 95% confidence threshold in idealized Gaussian noise for each wide band (2F about 75.6, 76.5, and 76.4 respectively) and recorded which jobs exceeded it (4, 6, and 36 jobs respectively).Then we visually inspected histograms of the surviving jobs by the same criteria as in Aasi et al. (2015), looking for fat tails rather than the thin tails indicative of injected signals.All histograms but one were rejected at a glance, and the remaining one (with 2F = 79.8 and ω = 0.001) hinted at abnormality on closer inspection.Also, that job searched around 998.4 Hz, a frequency known to be contaminated (in LIGO data) by violin modes of the test mass suspensions (Covas et al. 2018;Davis et al. 2021).Nevertheless, we followed it up.That search job was rerun, keeping detailed information on all templates with 2F ≥ 40.Plotting 2F versus frequency showed multiple peaks indicative of noise lines.We then searched the entire frequency band of that job with double the integration time.A real signal would produce double the value of 2F, or about 160.The double-length followup search produced several 2F peaks only a little over 100, and every search job dramatically failed the Cramer-von Mises test (ω ∼ 0.05).Therefore we concluded that we found no astrophysical signals.
We then set upper limits on h 0 in 1 Hz bands, similar to Owen et al. (2022).We considered a population of signals with fixed h 0 and random values of intrinsic parameters (f, ḟ , f ) as well as extrinsic parameters described in Jaranowski et al. (1998).We estimated what h 0 would be detected at a rate of 90% with 2F larger than the largest non-vetoed search result in that upper limit band.A semi-analytic estimate of h 0 was checked by software-injecting 1000 signals per 1 Hz upper limit band.The Drill pipeline does this more efficiently by cutting down on disk input/output bottlenecks, corrects a minor inconsistency in the way vetoed bands were incorporated into the upper limits in previous analyses, and computes upper limits for all bands that were not too heavily vetoed.A heavily vetoed band is one where eliminating the search jobs exceeding the ω threshold vetoes more than 10% of the upper limit band, and therefore a 90% upper limit is not meaningful.These amounted to 60 upper limit bands out of the 1015 covered by the search, or about 6%.
The left panel of Fig. 1 displays our 90% confidence upper limits on h 0 as a function of frequency, except in the 60 heavily vetoed bands.The discontinuities at 125 and 450 Hz are caused by differences in integration times in the three search bands.The (red) horizontal solid lines in the left panel of Fig. 1 show the range of h age 0 from mass-quadrupole GW emission, Eq. (3), and the (green) horizontal dashed lines show the range for rmode GW emission, Eq. ( 4).Note that the lower red and green lines coincide.The observed upper limits on h 0 are less (better) than the average values of the indirect limits  3), and the (green) horizontal dashed lines show the range for r-mode GW emission, Eq. ( 4).(Note that the h age 0 lower limits coincide.)The right panel shows corresponding upper limits on the dimensionless neutron-star ellipticity ϵ and r-mode amplitude α.
h age 0 over the full frequency band.Our search places limits on GW emission from NS 1987A that are better than the strictest h age 0 estimates over the astrophysically most interesting part of the frequency band, 50-600 Hz.
The efficiency of our search can be expressed in terms of the sensitivity depth (Behnke et al. 2015) where S h is the harmonic mean strain noise power spectral density.For our search D is about 29 Hz −1/2 , 22 Hz −1/2 , and 18 Hz −1/2 in the low, middle, and high frequency bands respectively.This is comparable to or somewhat worse than Owen et al. (2022), as one would expect from the short integration times (Wette 2023).Upper limits on h 0 imply upper limits on the fiducial neutron star ellipticity ϵ (Jaranowski et al. 1998)  We have performed a search for continuous GWs from NS 1987A over a much wider range of frequencies than the only previous physically consistent search (Owen et al. 2022) using an improved code on improved data.While we did not detect any astrophysical signal, we set upper limits which constrained the behavior of NS 1987A over a much broader parameter space than ever before.
When translated into fiducial neutron star ellipticity or r-mode amplitude, our upper limits (at the highest frequencies) are an order of magnitude better than the highest frequency limits set by Owen et al. (2022).Our best ellipticities are just above 10 −5 .Such elastic deformations are possible for quark stars and quark-baryon hybrid stars (Owen 2005;Johnson-McDaniel & Owen 2013), and are skirting the maximum predicted for normal neutron stars (Morales & Horowitz 2022).For a magnetic deformation, our strain upper limits imply a limit on the internal magnetic field of a few times 10 15 or 10 14 G if the protons in the core are or are not superconducting respectively (Lander 2014;Ciolfi & Rezzolla 2013).Our best upper limits on r-mode amplitude are starting to enter the range of theoretical predictions (Bondarescu et al. 2009).
These comparisons show that our search had a sensitivity to GW emission from NS 1987A compatible with a variety of theoretical predictions, not just the most extreme ones as in Owen et al. (2022).We used simple coherent integrations of relatively short spans of O3 data.Searches of better data from the current Advanced LIGO and Virgo's fourth observing run (O4) to Cosmic Explorer (Evans et al. 2021;Gupta et al. 2023) will further improve on this sensitivity, especially with more sophisticated data analysis techniques (Wette 2023).This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA.This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation, 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 research was supported in part by NSF grant 2012857 to the University of California at San Diego, NSF grants 1912625 and 2309305 to Texas Tech University, and NSF grant 2110460 to the Rochester Institute of Technology.The authors acknowledge computational resources provided by the High Performance Computing Center (HPCC) of Texas Tech University at Lubbock (http://www.depts.ttu.edu/hpcc/).This paper has LIGO Document number P2300319.

Figure 1 .
Figure 1.The left panel displays observational 90% confidence upper limits on h0 from NS 1987A in 1 Hz bands as a function of frequency.The (red) horizontal solid lines show the range of h age 0 from mass-quadrupole GW emission, Eq. (3), and the (green) horizontal dashed lines show the range for r-mode GW emission, Eq. (4).(Note that the h age on the r-mode amplitude α defined byLindblom et al. (1998) usingOwen (2010), in these expressions are uncertain by roughly a factor of three due to uncertainties in the neutron star mass and equation of state.Moreover, general relativity complicates these expressions in ways that have not yet been calculated.These upper limit estimates on ϵ and α for NS 1987A from this search are shown in the right panel of Fig 1.

Table 1 .
from O3 in the form of 1800 s Data parameters used in this search.Times are UTC.