First Constraining Upper Limits on Gravitational Wave Emission from NS 1987A in SNR 1987A

We report on a search for continuous gravitational waves (GWs) from NS 1987A, the neutron star born in SN 1987A. The search covered a frequency band of 75-275 Hz, included a wide range of spin-down parameters for the first time, and coherently integrated 12.8 days of data below 125 Hz and 8.7 days of data above 125 Hz from the second Advanced LIGO observing run. We found no astrophysical signal. We set upper limits on GW emission as tight as an intrinsic strain of $2\times10^{-25}$ at 90\% confidence. The large spin-down parameter space makes this search the first astrophysically consistent one for continuous GWs from NS 1987A. Our upper limits are the first consistent ones to beat an analog of the spin-down limit based on the age of the neutron star, and hence are the first GW observations to put new constraints on NS 1987A.


INTRODUCTION
Due to the detection of neutrinos from the supernova (Bionta et al. 1987;Hirata et al. 1987), it has been known since 1987 that SN 1987A probably produced a neutron star. This "NS 1987A" is the youngest neutron star known near our galaxy, 51.4 kpc away in the Large Magellanic Cloud (Panagia 1999). After many years of unsuccessful searches for a pulsar or non-pulsing neutron star, indirect evidence has accumulated in recent years. Cigan et al. (2019) observed infrared emission from a relatively warm, compact blob of dust that Page et al. (2020) showed could be powered by a 30 yr old cooling neutron star. Greco et al. (2021) and Greco et al. (2022) argue that the hard x-ray emission indicates the presence of a pulsar wind nebula.
Continuous gravitational wave (GW) emission from NS 1987A has been suggested since Piran & Nakamura (1988). Such GWs could be produced by a nonaxisymmetric deformation of the neutron star, by free precession, or by long-lived r-mode oscillations, as summarized e.g. by Glampedakis & Gualtieri (2018).
The most recent search for continuous GW emission from NS 1987A used stochastic background methods to analyze data from Advanced LIGO and Virgo's first three observing runs (Abbott et al. 2021a). Like previous similar searches referenced therein, Abbott et al. (2021a) assumed a small spin-down for NS 1987A. Of order 10 −9 Hz/s, this spin-down is large by the standards of known pulsars but small because it implies an enormous unphysical accretion spin-up to counteract the spin-down implied by GW emission at the high level required for detection (see below). Wette et al. (2008) scoped out broad band continuous GW searches for supernova remnants (such as 1987A) where there is evidence for a neutron star but pulses are not observed, and hence a wide range of frequencies and spin-downs (time derivatives of the frequency) must be covered. (Continuous GW searches involve longer coherence times than stochastic background searches, and typically search over spin-down parameters as well as GW frequencies.) The approach of Wette et al. (2008) was first used on the young neutron star in Cas A (Abadie et al. 2010), and similar methods have been used on other likely locations of neutron stars not observed as pulsars (Aasi et al. 2013a(Aasi et al. , 2015Zhu et al. 2016;Sun et al. 2016;Abbott et al. 2017Abbott et al. , 2019Dergachev et al. 2019;Ming et al. 2019;Piccinni et al. 2020;Lindblom & Owen 2020;Millhouse et al. 2020;Papa et al. 2020;Jones & Sun 2021;Abbott et al. 2021bAbbott et al. , 2022aMing et al. 2022). Wette et al. (2008) defined a key figure of merit, an indirect limit on GW emission similar to the spin-down limit for pulsars-i.e., a best case amplitude of GW emission. Even when the spin-down is not known, one can assume that it has been dominated by GWs since the birth of the star and that it has spun the star down significantly from its birth frequency. This results in a frequency-independent limit based on the age of the star, extended by Owen (2010) (1) Here h 0 is a measure of GW amplitude called the intrinsic strain (Jaranowski et al. 1998), D is the distance to the neutron star, a is its age, I is its moment of inertia, and n = ff /ḟ 2 is its braking index. The braking index is about 5 or 7 if the GW emission is due to a corotating nonaxisymmetry (ellipticity) or an r-mode respectively. The fiducial moment of inertia in Eq. (1) is on the low end of the predicted range, and depending on the star's mass and the nuclear matter equation of state h 0 could go up by about a factor 2 (Abbott et al. 2021c), so the range of h age 0 for NS 1987A is about 2-4×10 −25 . Sun et al. (2016) performed the most recent search for NS 1987A that used continuous GW methods, using the setup of Chung et al. (2011), and summarize earlier searches for that star. Continuous wave methods are generally more sensitive than stochastic background methods but more computationally intensive. Since the Wette et al. (2008) wide parameter space was unfeasible for a source as young as NS 1987A (19 years old for the data used and requiring a fourth spin-down parameter), Chung et al. (2011) narrowed the search by introducing a detailed spin-down model. But this is less robust than a model that makes few assumptions like Wette et al. (2008), and even with a narrow parameter space Sun et al. (2016) did not place upper limits beating the indirect limit h age 0 . Recent all-sky surveys for continuous GWs such as Abbott et al. (2022b) and Dergachev & Papa (2022) reach h age 0 in the direction of NS 1987A but cover too small a spin-down range for it.
As of Advanced LIGO's second observing run (O2), NS 1987A was 30 years old. For that age and a 51.4 kpc distance (Panagia 1999), h age 0 is comparable to what recent GW searches of supernova remnants such as Lindblom & Owen (2020) have achieved using only two spindown parameters. The results of Wette et al. (2008) and Wette (2012) can be combined to estimate that a coherent search of O2 data can use only two spin-down parameters and surpass the sensitivity of h age 0 for a computing budget of order a million core-hours on a modern cluster.
Here we describe such a search, which detected no astrophysical signals but placed the first direct upper limits on GW strain from NS 1987A to beat the indirect limit h age 0 over a wide and physically consistent parameter space.

SEARCH METHODS
Since our search methods were much like those of Lindblom & Owen (2020) and similar papers, we only summarize highlights and changes here and direct the reader to Lindblom & Owen (2020) and references therein for details. Our input parameters and some derived parameters are given in Table 1.
We used LIGO open data (Vallisneri et al. 2015;Abbott et al. 2021d) from O2, the most recent data publicly available when we started our computational runs, in the form of 1800 s Short Fourier Transforms (SFTs). O2 included data from the Hanford, WA (H1) and Livingston, LA (L1) 4 km interferometers. Once the integration time spans for our search bands were determined by computational cost (see below), we selected the stretch of data for each span to maximize sensitivity, which is proportional to live time over the power spectral density (psd) of strain noise (Jaranowski et al. 1998).
The integration method was the multi-detector Fstatistic (Jaranowski et al. 1998;Cutler & Schutz 2005), which efficiently accounts for the modulation of longlived signals due to the rotation of the Earth. Because 2F is a quadrature of four matched filters, in stationary Gaussian noise it is drawn from a χ 2 distribution with four degrees of freedom. If a signal is present, the χ 2 is noncentral and the amplitude signal-to-noise ratio is roughly F/2.
We assumed that the demodulated signal frequency evolved in the solar system barycenter frame as where the reference time t 0 is the beginning of the observation and the parameters (f,ḟ ,f ) are evaluated at that time. That is, we assumed no binary companion to NS 1987A, no glitches during the spans of integration, little timing noise, and insufficient frequency drift to require a third derivative.
To choose the parameter space-i.e., ranges of (f,ḟ ,f )-we first split the search into low and high frequency bands divided at 125 Hz, roughly twice the spin frequency of the fastest known young pulsar (Marshall et al. 1998). The latter should be the frequency of most efficient emission of mass quadrupolar GWs. For a given frequency f, the ranges of (ḟ ,f ) were chosen the same way as in Lindblom & Owen (2020). That is, n minḟ with the braking index ranging from n min = 2 to n max = 7. These ranges correspond to a wide range of observed and predicted behaviors, and are consistent with the minimum spin-down required for a given h 0 -see, e.g., Owen (2010). Here κ is the ratio of GW frequency to spin frequency. Note that the minimum value of −ḟ is greater than the maximum value covered by all-sky surveys (Abbott et al. 2022b (1), for r-mode emission from a low mass (and moment of inertia) star or mass quadrupole emission from an intermediate mass star, about 2.5 × 10 −25 . For a computational cost of 10 6 core-hours per band, this resulted in the integration times and other parameters shown in Table 1.
The code was an improved version of that used in Lindblom & Owen (2020), based on the S6SNRSearch tag of the LALSuite software package (LIGO Scientific Collaboration 2020) and its implementation of the F-statistic. The search parameter space was split into roughly 10 5 batch jobs per band, each taking roughly ten hours on the Texas Tech supercomputing cluster "Quanah." We also used these jobs as an ad hoc way of clustering candidate signals (see below).
We did not a priori veto candidates based on timefrequency behavior or lists of known instrumental lines. Due to the rapid spin-down a detectable signal would have, most templates overlap a spectral line for some time and such vetoes would render much of the search band unusable for setting upper limits. The rapid spindown has the advantage, however, of diluting the effect of a narrow line on any given template, as the template relatively rapidly moves out of the disturbed frequency band. This is shown by the relatively small number of candidates (see below). We did use the interferometer consistency veto, a simple check first used in Aasi et al. (2013b) that the joint 2F is greater than the value from either interferometer alone.
We performed consistency checks as in Lindblom & Owen (2020), plus two more necessitated by the unusual youth of the target and hence the high value of spin-downs searched: We confirmed that a third frequency derivative is not needed in Eq. (2) and that the standard SFT length of 1800 s is not too long.
Inspection of the parameter space metric (Wette et al. 2008) shows that omitting the third frequency derivative can result in a substantial mismatch between signal and template for the parameters and integration times used here. However Jaranowski & Krolak (2000) argued that correlations with lower derivatives allow for neglect of the third derivative at much longer integration times than we use, while still keeping a low mismatch and thus a high fraction of the ideal 2F. Essentially, correlations allow a large template bank to pick up the signal efficiently at a shifted position (f,ḟ ,f ). This argument is weak near the edges of parameter space, so we checked against a set of software injections (with the highest third derivatives allowed by our braking index range) and confirmed that omitting the third derivative causes no appreciable loss in 2F for a population of signals. As Jaranowski & Krolak (2000) argue, there is no detectable effect.
One might expect the SFT length to be a problem for frequency derivatives high enough to send a signal through multiple SFTs within the duration of one SFT, i.e. for |ḟ | of order 1/(1800 s) 2 or 3 × 10 −7 Hz/s-which is also the maximum |ḟ | covered by our search. The injection checks of our upper limits (see below) already test this to some extent, but we performed additional injection studies dedicated to this issue. We found no significant losses for values up to 5 × 10 −7 Hz/s, well beyond what we searched.

SEARCH RESULTS
We examined the search results for candidate signals surpassing 95% confidence in Gaussian noise, corresponding to 2F thresholds of 77.1 and 77.5 for the low and high frequency bands respectively. (These values were determined using an effective number of independent templates found by a Kolmogorov-Smirnov distance minimizer between the observed distributions of loudest events per search job and the Gaussian noise prediction.) The high band produced no candidates. The low band produced 29 search jobs with candidates mostly clustered around 83.32 Hz and 100.0 Hz, with sin-gle jobs at 107.1 Hz and 108.5 Hz. The highest 2F was 95. We examined the search jobs as in Lindblom & Owen (2020) and found that all 2F histograms and frequency plots showed contamination by broad noise lines. This was sufficient to rule out the candidates as astrophysical signals, but we followed up by checking against lists of known instrumental artifacts (Covas et al. 2018). The 83.32 Hz and 100.0 Hz clusters are due to known lines in L1 and H1 respectively, and the narrower cluster at 108.5 Hz is due to a known line in H1. The narrower cluster at 107.1 Hz has 2F improbably dominated by H1. The large clusters produced up to 0.1% candidates triggering the interferometer consistency veto, several orders of magnitude above typical search jobs. Therefore we do not claim any astrophysical signals in our search.
In the absence of signals, we set new astrophysically meaningful upper limits on intrinsic strain h 0 as a function of GW frequency in 1 Hz bands. The method was the same as in Lindblom & Owen (2020), with the same false dismissal rate of 10% (90% confidence) integrated over a population of sources with randomly oriented spin axes: A semi-analytic estimate was confirmed with 1000 software-injected signals per upper limit band.
The left panel in Fig. 1 displays our upper limits on h 0 as a function of frequency, minus two bands starting at 222 Hz and 264 Hz which the injections indicated had slightly more than 10% false dismissal rate. The discontinuity at 125 Hz is due to the difference in integration times used above and below that frequency. The (red) horizontal line in the left panel represents the fiducial value 2.5 × 10 −25 of the indirect limit h age 0 from energy conservation. Our search places limits on GW emission from NS 1987A that are significantly better than this astrophysical upper limit over the band searched, and better than the strictest h age 0 of 2.0 × 10 −25 over part of the band.
The efficiency of our search can be expressed in terms of two quantities derived from h 0 . One common measure (Wette et al. 2008) is the factor Θ in where S h is the harmonic mean psd of all SFTs and T is the total live time of data used. Our Θ is about 36 in the low frequency band and 40 in the high frequency band, slightly worse (higher) than the youngest supernova remnant searches so far-see Lindblom & Owen (2020) for example. The sensitivity depth, defined by (Behnke et al. 2015) as is about 36 Hz −1/2 and 29 Hz −1/2 for the low and high bands respectively. This is also slightly worse (lower) than the youngest supernova remnant searches so far, as expected due to the extreme youth and wide parameter space of NS 1987A. Upper limits on h 0 can be converted to upper limits on fiducial neutron star ellipticity using (Wette et al. 2008, e.g.) 9.5 × 10 −5 h 0 1.2 × 10 −24 and to upper limits on a particular measure of r-mode amplitude, α, (Lindblom et al. 1998) using Owen (2010, The numerical values are uncertain by a factor of roughly two or three due to uncertainties in the unknown neutron star mass and equation of state. Upper limits on and α for NS 1987A from this search are shown in the right panel of Fig 1. Indirect limits on and α derived from h age 0 are not shown since, on this logarithmic scale, they are close to the direct observational limits. Our upper limits on h 0 beat the indirect limit h age 0 . The latter limit is astrophysically interesting despite the youth of the source. Equation (1) is derived under the assumption that the star has spun down significantly since birth. Without that assumption, but with constant braking index n, the limit h age 0 is multiplied by (Sun et al. 2016) where f b is the GW frequency at birth. If f f b , this factor is 1 and we recover Eq. (1). If f b ≈ f as assumed by Sun et al. (2016), essentially imposing a smaller ellipticity, this factor can be much less than 1, but it does not need to be. It is straightforward to insert our direct limits on and α and integrateḟ to find that our assumptions are consistent-i.e., that these values of and α result in significant spin-down over 30 years and hence Eq. (1) is valid despite the youth of NS 1987A. Thus our search had a chance of detecting a signal, and the lack of detection represents the first GW observational constraints on NS 1987A.

CONCLUSIONS
We have performed the first search for GWs from NS 1987A that covered a physically consistent range of spin-downs and achieved a sensitivity better than the indirect limit from energy conservation. We also showed that this limit is applicable to NS 1987A despite the youth of the source. While we detected no astrophysical signal, we set direct observational upper limits that beat the indirect limit and thus for the first time constrain the GW emission of NS 1987A if it is emitting within the frequency band searched. Our constraints on r-mode amplitude are not competitive with the standard theoretical prediction (Bondarescu et al. 2009), but our constraints on ellipticity are within the predicted range of elastic deformations of quark stars (Owen 2005). If NS 1987A is made of baryonic matter and the protons in its core are not yet superconducting, our constraints imply upper limits on the internal magnetic field of about 10 15 G for the twisted torus configuration likely formed with the neutron star (Ciolfi & Rezzolla 2013). If the protons are now superconducting, the field is likely mainly poloidal and our constraints limit the field to less than about 10 16 G (Lander 2014).
Our search achieved this with a simple coherent integration of O2 data. Better methods and data are available, and we expect this will motivate further searches and improvements of search methods for rapidly evolving continuous wave signals.