Constraints from LIGO O3 Data on Gravitational-wave Emission Due to R-modes in the Glitching Pulsar PSR J0537–6910

We use data from the third observing run (O3) of the two Advanced LIGO detectors (Aasi et al. 2015a). The O3 data run lasted from 2019 April 1 to 2020 March 27, with a one-month pause in data collection in October 2019 (see Figure 1). The Hanford (H1) and the Livingston (L1) detectors had duty factors of ∼76% and ∼77%, respectively during O3. For the LIGO O3 data set, the analysis uses the calibration with estimated upper limits on amplitude and phase uncertainties of ∼7% and ∼4°, respectively (Sun et al. 2020), which we use as conservative estimates of the true calibration uncertainty near the frequencies analyzed here.

Zoom In Zoom Out Reset image size

The timing analysis of J0537 using NICER data from 2017 August 17 to 2020 April 25 is described in Ho et al. (2020). During this timespan eight glitches were detected, with the last three occurring during the O3 run (see Figure 1). These three glitches divide up the O3 data set into four inter-glitch epochs, which form natural segments into which to divide the data analyses in the searches described below. Note that parameters of a revised ephemeris for epoch 4, which are presented in Abbott et al. (2020), are within the errors of the ephemeris from Ho et al. (2020) and thus do not affect our analysis here. In the first of our two search methods to be presented (${ \mathcal F }$/${ \mathcal G }$-statistic search), the data are indeed divided into four segments in this way. In the second search described (5-vector method), the long ∼170 day inter-glitch epoch is divided into two for computational reasons, giving a total of five data segments in the analysis. More detail will be given below.

The r-mode GW emission frequency f depends on the pulsar spin frequency ν and on the neutron star structure (e.g., Idrisy et al. 2015). Following the analysis of Caride et al. (2019), the relationship can be expressed as

$$\begin{eqnarray}&&f=A\nu -B\left(\displaystyle \frac{{\nu }^{2}}{{\nu }_{K}^{2}}\right)\nu ,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\dot{f}=A\dot{\nu }-3B\left(\displaystyle \frac{{\nu }^{2}}{{\nu }_{K}^{2}}\right)\dot{\nu },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\ddot{f}=A\ddot{\nu }-\left(3+\displaystyle \frac{6}{n}\right)B\displaystyle \frac{{\nu }^{2}}{{\nu }_{K}^{2}}\ddot{\nu }.\end{eqnarray} \tag{ 3 }$$

where the dots denote time derivatives of the frequencies and ν_K is the Keplerian frequency for which centrifugal forces would tear the star apart. These relations make allowance for fully relativistic gravity (via the A-parameter) and for rotation of the star (via the B-parameter) and are thus likely to be accurate as long as no other physical mechanism plays a significant role in the r-mode oscillation, e.g., a resonance between the core and crust oscillation frequencies (see Levin & Ushomirsky 2001). The precise way in which these relations are used in this paper differs slightly between the two searches we present, as will be described below.

The values of A and B for any particular neutron star depend both upon its (unknown) mass, and also upon the (unknown) equation of state (EoS) of dense matter. Caride et al. (2019) considered a range of astrophysically motivated neutron star masses, and a range of microphysically motivated equations of state, to suggest a conservative range of values of A and B to employ in a search:

$$\begin{eqnarray}&&1.39\lt A\lt 1.57,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&0\lt B\lt 0.195.\end{eqnarray} \tag{ 5 }$$

One also needs a value of the (unknown) Keplerian frequency ν_K to make use of the above equations. We again follow Caride et al. (2019) in making a conservative choice, taking the (low) value ν_K ≃ 506 Hz. This value is obtained by assuming that the Keplerian frequency coincides with the frequency of the most rapidly rotating neutron star currently known, 716 Hz for PSR J1748-2446ad (Hessels et al. 2006), and that J0537 is a low-mass star, with a consequently lower ν_K ; see Caride et al. (2019) for relevant discussion.

The GW signal from an isolated rotating neutron star is described by the model derived in Jaranowski et al. (1998). The signal depends on four extrinsic parameters and 3 + s intrinsic parameters, where s is the number of spin-down parameters included in the model. The four extrinsic parameters are the GW amplitude h₀, the star's spin axis inclination angle (with respect to the line of sight) ι, the GW polarization angle ψ, and a constant phase ϕ₀. The 3 + s intrinsic parameters are the frequency f₀, s spin-down parameters f₁,...,f_s , and two angles describing the position in the sky of the source, i.e., R.A. α and decl. δ. For GW emission from r-modes, the only difference from the model of Jaranowski et al. (1998) is that the polarization angle ψ needs to be transformed as

$$\begin{eqnarray}&&\psi \to \psi +\pi /4.\end{eqnarray} \tag{ 6 }$$

The derivation of this transformation is given and discussed in detail in Owen (2010). The position in the sky of J0537 is precisely known (J2000; Townsley et al. 2006):

$$\begin{eqnarray}&&\alpha ={05}^{{\rm{h}}}{37}^{{\rm{m}}}47.{}^{{\rm{s}}}416,\delta =-69^\circ 10^{\prime} 19\buildrel{\prime\prime}\over{.} 88.\end{eqnarray} \tag{ 7 }$$

The intrinsic parameters determine the phase modulation ϕ(t) of the signal, which in the model of Jaranowski et al. (1998) is approximated by

$$\begin{eqnarray}&&\phi (t)\cong 2\pi \sum _{k=0}^{s}{f}_{k}\displaystyle \frac{{t}^{k+1}}{(k+1)!}+\displaystyle \frac{2\pi }{c}{\boldsymbol{n}}(\delta ,\alpha )\cdot {{\boldsymbol{r}}}_{{\rm{d}}}(t)\sum _{k=0}^{s}{f}_{k}\displaystyle \frac{{t}^{k}}{k!},\end{eqnarray} \tag{ 8 }$$

where r _d(t) is the position of the detector in the solar system barycenter frame and n is the constant unit vector in the direction of the neutron star in that frame and is

$$\begin{eqnarray}&&{\boldsymbol{n}}=(\cos \delta \cos \alpha ,\cos \delta \sin \alpha ,\sin \delta ).\end{eqnarray} \tag{ 9 }$$

This formula takes into account both the intrinsic spin-down of the source and the Doppler modulations produced by the motion of the detector.

In addition to providing information on timing, electromagnetic observations also provide insight into the likely orientation of the spin axis of the pulsar. J0537 is associated with the supernova remnant N157B in the Large Magellanic Cloud. Observations of the pulsar wind nebulae (PWN) torus of the remnant have enabled an accurate measurement of the orientation of the axis of symmetry of the torus (Ng & Romani 2008). It is very likely that the symmetry axis of the torus coincides with the spin axis of the pulsar. This enables us to estimate the values of the two parameters of the GW signal: polarization angle ψ and inclination angle ι. For the former, we have ψ = 2.2864 ± 0.0384 rad. From the observations, we cannot obtain the absolute direction of the rotation of the pulsar: we have two possible values of the inclination angle ι or π − ι (Jones 2010). Thus the measurements yield ι = 1.522 ±0.016 rad or ι = 1.620 ± 0.016 rad. This information can be used to carry out versions of our searches where this prior information on these two angles is incorporated into the search, as will be described below. An inclination angle ι close to π/2 means that the GW signal is almost linearly polarized and consequently it has nearly the smallest maximum signal-to-noise ratio (S/N) achievable.

Detection statistics. We perform a coherent search of time domain segments of the data spanning between the glitches of the pulsar. The coherent search is performed using the ${ \mathcal F }$-statistic derived in Jaranowski et al. (1998). The ${ \mathcal F }$-statistic has been used extensively in searches for GWs from rotating neutron stars, including the recent directed searches for GWs from r-mode emission by Fesik & Papa (2020a, 2020b). In this search, we use an efficient implementation of the ${ \mathcal F }$-statistic presented in Astone et al. (2010) and Aasi et al. (2014). Its two main features are the use of the fast Fourier transform (FFT) to evaluate the ${ \mathcal F }$-statistic over the whole band of the data analyzed and the use of a grid of templates based on optimal coverings of the parameter space. The FFT computation is enabled by a technique called resampling which consists of interpolating the data to Solar System Barycenter time (see Section IIID of Jaranowski et al. 1998 for details). The resampling needs to be performed for all frequency derivatives only once for each sky position that we analyze. Thus in the directed search presented here, where the pulsar's position is known to a very high accuracy, the resampling needs to be performed only once.

If the angles ψ and ι are known, there is a variant of the ${ \mathcal F }$-statistic called the ${ \mathcal G }$-statistic (Jaranowski & Królak 2010) which implements the coherent search in this case. As described in Section 3.2 above, for J0537 the inclination angle ι is very close to π/2, which indicates that the GW signal is almost linearly polarized. We perform a Monte Carlo simulation where, for an observation time of 60 days, we add simulated signals to Gaussian noise with the two possible values of the inclination angle ι of J0537 consistent with observations of the PWN. We then search for the signal in both cases with a template for ι = 1.522. We perform the simulation for an array of signal-to-noise ratios (S/Ns) from 5 to 25. We find no appreciable difference in the detectability of the signals and accuracy of the estimation of the frequency parameters for the two cases. The S/N of the signals recovered using the template with a "wrong" inclination angle decreased by no more than 2% and the accuracy of the estimations of frequency and spin-down rate decreased by no more than 2%. Thus we conclude that only one value of the inclination angle is needed when using the ${ \mathcal G }$-statistic method for the search.

In this paper, we analyze the data using both the ${ \mathcal F }$- and ${ \mathcal G }$-statistic. We search the data from the two LIGO detectors and consequently use a multi-detector implementation of the statistics developed in Królak et al. (2004, 2007), Jaranowski & Królak (2009), and Cutler & Schutz (2005).

Parameter space. The relation between stellar spin frequency ν and r-mode GW frequency f is described in Equations (1)–(3), following the model of Caride et al. (2019), who used these relations to propose ranges of GW frequencies and frequency derivatives over which to search:

$$\begin{eqnarray}&&\nu \left({A}_{\min }-{B}_{\max }\displaystyle \frac{{\nu }^{2}}{{\nu }_{K}^{2}}\right)\leqslant f\leqslant \nu \,{A}_{\max },\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&-\dot{\nu }\left(\displaystyle \frac{f}{\nu }-2{B}_{\max }\displaystyle \frac{{\nu }^{2}}{{\nu }_{K}^{2}}\right)\leqslant -\dot{f}\leqslant -\dot{\nu }\displaystyle \frac{f}{\nu },\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&0\leqslant \ddot{f}\leqslant \ddot{\nu }\displaystyle \frac{f}{\nu },\end{eqnarray} \tag{ 12 }$$

where $\ddot{\nu }$ is the maximum value consistent with electromagnetic observations. The value used for ν_K and the maximum and minimum values we assume for A and B are given in Section 3.1. The parameter space for the frequency and its two derivatives defined by the above equations is illustrated in Figure 2.

Zoom In Zoom Out Reset image size

In the J0537 r-mode search of Fesik & Papa (2020a, 2020b) the prescription for the parameter space proposed by Caride et al. (2019) was also used; however, their assumed ranges of $\dot{f}$ and $\ddot{f}$ were independent of the frequency f.

Grid of templates. The problem of constructing a grid in the parameter space is equivalent to the covering problem (Conway & Slone 1999; Prix 2007). There exist optimal coverings that have the smallest possible number of grid points per unit volume of the parameter space. Our implementation of the ${ \mathcal F }$- and ${ \mathcal G }$-statistic uses an FFT which evaluates the statistic at Fourier frequencies which in general do not coincide with nodes of the optimal grid. To address this problem, we introduce constrained grids compatible with application of the FFT (Astone et al. 2010). For the case of the directed search that we present here, in which the position of the source in the sky is known, the construction of the constrained grids is described in Pisarski et al. (2011).

One important problem is to determine how many spin-down parameters f_k (see Equation (8)) we need to take into account in the search. This problem is addressed in Brady et al. (1998). In general, a longer observation time T requires the inclusion of more spin-down parameters in our templates so that we can match the signal accurately. Thus we calculate the number of grid points N_k (T) for a given number k of spin-down parameters included as a function of observation time T. The observation time at which we need to include an additional f_{k + 1} spin-down parameter is the observation time at which the curves N_k (T) and N_{k + 1}(T) intersect (see Figure 6 of Brady et al. 1998 for illustration).

The results for this search are shown in Figure 3. We see that the number of grid points N₂(T) when f₂ is included becomes larger than the number of grid points N₁(T) when only f₁ is included for observation times T greater than 117 days. In this analysis we assume a minimum spin-down age of 5000 yr, appropriate for the pulsar J0537. Thus for the expected range of frequency parameters of the GW signal from r-modes from J0537, we need to include f₂ in the templates only for coherent observation times longer than 117 days.

Zoom In Zoom Out Reset image size

To support this result, we perform a Monte Carlo simulation for an observation time of 66 days where we compare the detectability of the GW signal with one spin-down and two spin-down parameters using the ${ \mathcal F }$-statistic with templates which include only one spin-down parameter. The simulation is performed for Gaussian noise and for an array of S/Ns from 5 to 25. We find no appreciable difference in the detectability of the signals and in the accuracy of estimation of the frequency and the first spin-down parameter, whether or not the second spin-down parameter is included in the simulated signals. The S/N of the signals recovered when two spin-down parameters are included decreased by no more than 1% with respect to signals when only one spin-down parameter is included. The accuracy of the estimation of the frequency and the spin-down parameter decreased by no more than 1%.

Threshold and false alarm probability. In order to assess the significance of the candidates obtained during the search, we need to calculate the false alarm probability (FAP). We follow the approach developed in Jaranowski & Królak (2009; see Chapter 6.1.3). In that approach, the parameter space that we search is divided into a number N_c of cells. The cells are defined in such a way that the autocovariance function of the statistic between a point inside a given cell and outside that cell is less than $\tfrac{1}{2}$. This means that the statistic between cells is approximately uncorrelated and the cells can be considered as independent. Thus the probability P^T _F that a statistic ${ \mathcal L }$ exceeds a threshold ${{ \mathcal L }}_{o}$ in one or more cells is given by

$$\begin{eqnarray}&&{P}_{F}^{T}({{ \mathcal L }}_{o})=1-{\left[1-{P}_{F}({{ \mathcal L }}_{o})\right]}^{{N}_{c}},\end{eqnarray} \tag{ 13 }$$

where ${P}_{F}({{ \mathcal L }}_{o})$ is probability distribution of the statistic at a given point in the parameter space assuming that data is only noise. This is by definition the FAP. The number of cells N_c is defined as

$$\begin{eqnarray}&&{N}_{c}=\displaystyle \frac{V}{{v}_{c}},\end{eqnarray} \tag{ 14 }$$

where V is the volume of the parameter space and v_c is the volume of the cell. By inverting Equation (13), we can obtain the threshold ${{ \mathcal L }}_{o}$ for a chosen value of the FAP. To calculate the thresholds using Equation (13), we assume Gaussian noise, i.e., we assume that the $2{ \mathcal F }$- and $2{ \mathcal G }$-statistics have central χ² distribution with four and two degrees of freedom, respectively.

The search setup. During the O3 run, three glitches of J0537 are detected, as described in Section 2.2. These are glitches 6–8 presented in Tables 1 and 2 of Ho et al. (2020) (see also Abbott et al. 2020). For our analysis, we select four data segments between the times of the three glitches of the pulsar, with the first segment starting at the beginning of the O3 run and the last segment ending at the end of O3. We select data segments that have length equal to an integer multiple of sidereal days. The start times of the segments after the glitches that occurred in O3 are chosen to be 1 day after the glitch and the end times around 1 day before the glitch. This results in data segments 62, 168, 59, and 66 sidereal days long.

Table 1. Parameters of the Search of the Four Time Domain Segments Analyzed Using the ${ \mathcal F }$/${ \mathcal G }$-statistic Method

Segment	1	2	3	4
t0 (GPS)	1238163456	1243641549	1258329549	1263599949
Tspan (days)	62	168	59	66
${{ \mathcal F }}_{o}$	12.5	15	12.5	12.5
${{ \mathcal G }}_{o}$	10	12.5	10	10
δ f (Hz)	1.19 × 10−7	5.96 × 10−8	1.19 × 10−7	1.19 × 10−7
$\delta \dot{f}$ (Hz s−1)	2.21 × 10−14	1.30 × 10−14	2.28 × 10−14	1.43 × 10−14
$\delta \ddot{f}$ (Hz s−2)	⋯	3.14 × 10−21	⋯	⋯
Nt	5.02 × 109	6.09 × 1010	4.54 × 109	7.74 × 109
${ \mathcal F }$ 0.01	26.2	26.3	26.5	26.9
${{ \mathcal G }}_{0.01}$	22.9	23.0	23.2	23.6

Download table as:  ASCIITypeset image

The initial GPS times t₀ of each data segment are given in Table 1. During the longest 168 day segment we were observing for only 138 days because of the commissioning break in October 2019. In each time domain segment, we search for GW emission due to r-modes in the parameter space determined by Equations (10), (11), and (12). We adopt values of the spin parameters ν, $\dot{\nu }$, and $\ddot{\nu }$ of J0537 from Table 1 of Ho et al. (2020). For segments 62, 168, 59, and 66 days long, we select spin parameters for epochs at MJD 58600, 58723, 58836, and 58918, respectively. We include the second spin-down derivative only for 168 day segments because, as shown above, for the other segments the second spin-down parameter can be ignored. The search is performed using both the ${ \mathcal F }$- and ${ \mathcal G }$-statistic. For the ${ \mathcal G }$-statistic, we fix angles ψ and ι at

$$\begin{eqnarray}&&\psi =2.2864,\iota =1.522.\end{eqnarray} \tag{ 15 }$$

In Table 1, we give thresholds ${{ \mathcal F }}_{o}$ and ${{ \mathcal G }}_{o}$ used for the ${ \mathcal F }$- and ${ \mathcal G }$-statistic search, respectively. We also present grid spacings δ f, $\delta \dot{f}$, and $\delta \ddot{f}$ in frequency and two spin-down parameters as well as the number N_t of templates used. Moreover, in Table 1, we show thresholds ${{ \mathcal F }}_{0.01}$ and ${{ \mathcal G }}_{0.01}$ for the ${ \mathcal F }$- and ${ \mathcal G }$-statistic, respectively, corresponding to an FAP of 1% for searches in each individual segment. These thresholds are used to select candidates from the search for follow-up analysis.

The 5-vector method has its foundation in the coherent pipeline for the narrowband search of continuous waves (CWs) from known pulsars (Astone et al. 2014; Mastrogiovanni et al. 2017), which has been applied to the analysis of both LIGO and Virgo data (Aasi et al. 2015b; Abbott et al. 2017, 2019). Narrowband searches are motivated by the need to take into account a possible mismatch between the actual GW frequency and the value inferred from electromagnetic observations. The mismatch could arise from errors in the electromagnetic measurements or in the presence of differential rotation if the GW signal is mainly emitted by the inner parts of the star, rather than by the crust, or if free precession is excited. Since the EoS is uncertain (see Section 3), so is the GW frequency. As a consequence, we need to search over a wide range of frequencies (and corresponding spin-down rates).

The 5-vector r-mode search is based on an extension of the standard narrowband pipeline. The main reason why the standard pipeline is not directly applicable for this search is the presence of glitches. The 5-vector pipeline is a fully coherent algorithm: it assumes that the phase evolution is regular and therefore a Taylor expansion involving the frequency and its derivatives can be defined (Mastrogiovanni et al. 2017). In the presence of glitches, this hypothesis is only valid within each intra-glitch period. In fact, it is not possible to perform a coherent analysis along the whole GW observing run, since there are phase discontinuities at the times of each glitch. To overcome this problem, a semi-coherent search is performed: intra-glitch periods are studied coherently and independently, and the results are later combined incoherently. This is only possible thanks to NICER observations, which provided the glitch times.

Other relevant modifications are implemented in order to deal with a parameter space larger by a factor ∼10⁴ with respect to standard narrowband searches. This is due to an enlargement factor of ∼200 for the frequency and ∼50 for spin-down rate. Thus, optimization schemes to manage the computational load have been developed.

Parameter Space. The equations relating r-mode GW frequency to spin frequency are given by Equations (1)–(3), following the model of Caride et al. (2019). The ranges in f and $\dot{f}$ are obtained by inverting Equations (1) and (2) in terms of A and B:

$$\begin{eqnarray}&&A=\displaystyle \frac{1}{2}\left(3\displaystyle \frac{f}{\nu }-\displaystyle \frac{\dot{f}}{\dot{\nu }}\right),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&B=\displaystyle \frac{{\nu }_{K}^{2}}{2{\nu }^{2}}\left(\displaystyle \frac{f}{\nu }-\displaystyle \frac{\dot{f}}{\dot{\nu }}\right).\end{eqnarray} \tag{ 17 }$$

The region to investigate in the parameter space of f and $\dot{f}$ is defined by the corresponding values of A and B that satisfy the relations 1.39 ≤ A ≤ 1.57 and 0 ≤ B ≤ 0.195, as described in Section 3.1. These equations produce a parallelogram shape in parameter space, as shown in Figure 12 below, and reduce the volume to be searched over by more than a factor of 30 with respect to the choice of a rectangular grid in $(f,\,\dot{f})$. The $(f,\,\dot{f})$ region covered by this search is slightly different, at the corners, from that of the ${ \mathcal F }/{ \mathcal G }$ method, as can be seen comparing Figures 2 and 12. The central step of the 5-vector pipeline is the application of spin-down rate corrections. Then, to cover this region in the parameter space a different frequency band is selected for each value of the spin-down rate. The correct frequency range for each spin-down rate is obtained by solving Equation (2) for A and substituting it into Equation (1) in order to get

$$\begin{eqnarray}&&f=\displaystyle \frac{\dot{f}}{\dot{\nu }}\nu +2B\left(\displaystyle \frac{{\nu }^{2}}{{\nu }_{K}^{2}}\right)\nu .\end{eqnarray} \tag{ 18 }$$

Using the assumed ranges of A and B, we search the band 86–97 Hz, with a frequency step of δ f = 1/T, while the spin-down is discretized in bins of width $\delta \dot{f}=1/{T}^{2}$. Note that the second derivative $\ddot{f}$ can be considered as constant (within one frequency bin) for periods of length T < T_max, with ${T}_{\max }={\left(2/{\ddot{f}}_{\max }\right)}^{1/3}$, where ${\mathop{f}\limits^{\text{\unicode{x02025}}}}_{max}=\mathop{\nu }\limits^{\text{\unicode{x02025}}}({A}_{max}-{B}_{min}\left[\tfrac{{\nu }^{2}}{{\nu }_{K}^{2}}\right])$ is the maximum value allowed for the second derivative of the GW frequency. Assuming a braking index ≈7, we can write $\ddot{\nu }=7{\dot{\nu }}^{2}/\nu$ which gives T_max ≈ 75 days. These considerations lead to the choice of the data segments to be studied coherently and independently, as reported in Table 2. Segments 2 and 3 actually correspond to the same intra-glitch period (same rotational parameters), which we split in order not to exceed the maximum length T_max. Therefore, the second derivative is fixed for all segments for the 5-vector pipeline to the value $\ddot{f}=\bar{A}\ddot{\nu }=1.48\times {10}^{-20}\,\mathrm{Hz}\,{{\rm{s}}}^{-2}$, where $\bar{A}$ is the average value for the first-order parameter A (second-order corrections are neglected to fix this value). We fix $\ddot{\nu }=1\times {10}^{-20}\,\mathrm{Hz}\,{{\rm{s}}}^{-2}$ for all segments, which is consistent with NICER measurements reported in Ho et al. (2020). The second-order spin-down bin is $\delta \ddot{f}\simeq 2/{T}^{3}=1.30\times {10}^{-20}\,\mathrm{Hz}\,{{\rm{s}}}^{-2}$ and is large enough to include the uncertainty in $\ddot{f}$ due to the unknown exact values of parameters A and B and to the measurement of $\ddot{\nu }$. The main search parameters are shown in Table 2.

Table 2. Parameters of the Five Time Domain Periods Analyzed Using the 5-vector Pipeline with $\delta \dot{f}=\delta {f}^{2}$

Period	MJD Start	MJD End	MJD Epoch	Duration	ν [Hz]	$\dot{\nu }$ [10−10 Hz s−1]	Number of Points	δ f [10−7 Hz]
1	58574	58629	58600	55	61.9145	−1.9974	1.54 × 109	2.1898
2	58646	58701	58723	55	61.9124	−1.9973	1.63 × 109	2.1492
3	58701	58757	58723	56	61.9124	−1.9973	1.63 × 109	2.1492
4	58810	58863	58836	53	61.9104	−1.9974	1.29 × 109	2.3212
5	58873	58935	58918	62	61.9090	−1.9977	2.35 × 109	1.9026

Note. The chosen time periods are the same over which NICER provides the pulsar ephemeris.

Download table as:  ASCIITypeset image

Search pipeline description. The full search band of 86–97 Hz is divided into 1 Hz intervals for computational reasons. For each interval, the coherent part of the search pipeline, which consists of different steps (Astone et al. 2014; Mastrogiovanni et al. 2017), is applied. Schematically: application of the Doppler and Einstein delay corrections; ²⁹³ application of spin-down corrections for each bin of the spin-down grid via a heterodyning signal processing technique; computation of network detection statistic using 5-vectors. These are complex arrays containing the Fourier transform of the data ²⁹⁴ at the five frequencies at which the power of a CW signal, emitted at a given frequency, would be split by sidereal modulation. In order to reduce the effect of spectral leakage due to frequency discretization, which can produce a sensitivity loss of up to ∼36%, we use an "interbinning" procedure (Mastrogiovanni et al. 2017). This consists of estimating the data FFT values at the halfway bin $k+\tfrac{1}{2}$ in terms of the values at the kth and (k + 1)th bins as

$$\begin{eqnarray}&&{\tilde{x}}_{{FFT},k+\tfrac{1}{2}}\approx \displaystyle \frac{\pi }{4}\left({\tilde{x}}_{k}-{\tilde{x}}_{k+1}\right).\end{eqnarray} \tag{ 19 }$$

The detection statistic is computed separately on the natural grid and on the grid of shifted bins.

The incoherent part consists of tracking the evolution of the GW signal through the glitches. For each point of the parameter space of the first segment with GW parameters $(f,\dot{f})$, we first compute the corresponding (A, B) using Equations (16) and (17). Then keeping A and B constant, the GW frequency/spin-down rate pair $(f,\dot{f})$ in the other segments is constructed by taking into account the change of rotational parameters due to glitches. In practice, a different grid in the plane $(f,\dot{f})$ needs to be used for each segment, as the duration of each segment is in general different (see Table 2). When a pair $(f,\dot{f})$ is extrapolated from one segment to the following segment, the pair corresponding to the point of the grid closest to the extrapolated value is taken. At the end, for each point in the parameter space of the first segment, a weighted average of the corresponding detection statistic over the five analyzed segments is computed. The weights are computed by the inverse of the median of the data power spectrum in the 1 Hz band to which the point belongs. The fundamental assumption we make in this procedure is that the coefficients A and B do not change across the glitches. As their value is related to the neutron star EoS and thus should not be affected by glitches, this assumption is reasonable.

Once the overall weighted statistic is computed for each point of the parameter space, the loudest outliers (i.e., those having the highest value of the detection statistic) are selected every 1 mHz over the corresponding whole spin-down range. For each outlier, the significance is computed through the p-value, obtained by comparing the outlier statistic value to the statistic noise distribution, following the procedure used in Abbott et al. (2017, 2019). More specifically, the noise distribution is computed from the data itself, taking the 1 Hz band which includes the outlier, with the exclusion of the 1 mHz band to which the outlier belongs, plus the preceding and the following 1 Hz bands. Outliers with p-value smaller than 1% (after properly taking into account the trial factor) are upgraded to the level of candidates and followed up in order to increase their significance or, conversely, to demonstrate they are incompatible with an astrophysical GW signal.

Candidates. We performed the search in LIGO data for GW emission from r-modes in J0537 during the four time domain segments using both the ${ \mathcal F }$- and ${ \mathcal G }$-statistic methods. In the frequency band searched, the data contain several known periodic interferences that correlate well with our templates. To identify the origin of the interferences in each case, we also searched the data from the Hanford and Livingston detectors separately. In Figure 4, we present results of the search for the case of 66 day segments using the ${ \mathcal F }$-statistic. We see in this case several strong lines are present in the Hanford detector data. These lines are not present in Livingston detector data but contaminate the search of the network of detectors. Most of these lines are integer multiples of 1 Hz. As these single-detector lines cannot be of GW origin, we veto them before further analysis of the candidates. The lines vetoed in the case of the 66 day segment are marked by vertical lines in Figure 4. In searches of all data segments, the strong lines that led to threshold-crossing outliers only appeared in the Hanford detector data.

Zoom In Zoom Out Reset image size

In Figure 5 we present candidates obtained in our search of the four time domain segments using the ${ \mathcal F }$-statistic and in Figure 6 candidates using the ${ \mathcal G }$-statistic. These candidates are after vetoing the lines. The continuous horizontal lines correspond to the threshold for 0.1% false alarm probability. The threshold is calculated for the total number of cells in all the four time domain segments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also give thresholds corresponding to 1% FAP for the parameter space searched limited to each individual segment. These thresholds are drawn as dashed horizontal lines.

We see that none of the candidates is strong enough to cross the global threshold for the FAP of 0.1%. However there are several crossings of the 1% local threshold for the 168 day segment. These subthreshold candidates are listed in Table 3, together with values of their frequency and spin-down parameters. The S/Ns of the candidates are also given.

Table 3. Subthreshold Candidates for the ${ \mathcal F }$/${ \mathcal G }$-statistic Search

f [Hz]	$\dot{f}$ [Hzs−1]	$\ddot{f}$ [Hz s−2]	S/N
${ \mathcal G }$-statistic 168 days
94.7057	−3.04 × 10−10	1.04 × 10−20	6.8
${ \mathcal F }$-statistic 168 days
87.6870	−2.82 × 10−10	3.47 × 10−21	7.2
89.5434	− 2.88 × 10−10	6.94 × 10−21	7.2
88.9333	−2.86 × 10−10	3.47 × 10−21	7.0
88.5916	−2.85 × 10−10	3.47 × 10−21	7.0
91.3000	−2.94 × 10−10	6.94 × 10−21	7.1
94.2063	−3.04 × 10−10	3.47 × 10−21	7.1
96.5544	−3.11 × 10−10	0	7.1
96.1719	−3.09 × 10−10	3.47 × 10−21	7.1

Download table as:  ASCIITypeset image

Subthreshold candidates occur only in the 168 day segment. We searched for coincident candidates to the subthreshold candidates in the three other data segments, correcting the frequency and spin-down parameters of the candidates for glitches and frequency evolution. No coincident candidate is found, and therefore we conclude that we have no valid candidate for a GW signal in our search.

Upper limits and sensitivity of the search. As no significant candidate signal for GW emission from r-modes in J0537 is found, we set upper limits on the intrinsic GW amplitude h₀. We set upper limits for 0.5 Hz bands. In each band, we generate signals for an array of amplitudes h₀ for the position of J0537. For each amplitude, we generate 250 signals with f, $\dot{f}$, and $\ddot{f}$ (for the 168 day segments) chosen from uniform random distributions in their respective ranges. Moreover for the case of the ${ \mathcal F }$-statistic search, we also choose the values of the polarization angle ψ and cosine of the inclination angle ι from uniform random distributions. The signals are added to the real data and searches are performed with the same grids and search set-up as for the real data search in the neighborhood of simulated signal parameters. We search ∼±3 grid points for $\dot{f}$ and ∼±1 grid points for $\ddot{f}$ from the parameters of the simulated signal. A signal is considered as detected if the highest detection statistic value from the simulated signal search is higher than that obtained in the search of real data for a given 0.5 Hz band. The detection efficiency is the fraction of recovered signals. We estimate the ${h}_{0}^{90 \% }$, i.e., 90% confidence upper limit on the GW amplitude h₀, by fitting ²⁹⁵ a sigmoid function to a range of detection efficiencies E as a function of injected amplitudes h₀, $E({h}_{0})={\left(1+\exp (k({x}_{0}-{h}_{0}))\right)}^{-1}$, with k and x₀ being the parameters of the fit. Figure 7 presents an example fit to the simulated data, with the 1σ errors on the ${h}_{0}^{90 \% }$ estimate marked in red, for the evaluation of the ${ \mathcal F }$-statistic method in the 59 day segment at the mid-band frequency of 91.423 Hz.

Zoom In Zoom Out Reset image size

These upper limits mean that a GW signal with amplitude ${h}_{0}^{90 \% }$ from r-mode emission in J0537 would be detected in our search with 90% probability. Specifically, errors due to the Monte Carlo simulations with N_i = 250 injections per amplitude are less than 3% for the ${ \mathcal F }$-statistic and less than 1% for the ${ \mathcal G }$-statistic results. In Figure 8, we present our 90% upper limits for all four searches and for both ${ \mathcal F }$- and ${ \mathcal G }$-statistics, with the Monte Carlo simulation errors added in quadrature to the amplitude calibration uncertainty error of 7% (see Section 2.1); the latter dominates the uncertainty of our error estimate. The upper limits for the ${ \mathcal G }$-statistic are greater than for the ${ \mathcal F }$-statistic. This is because for the ${ \mathcal G }$-statistic the inclination angle ι is close to π/2 meaning that the GW signal is almost linearly polarized and consequently it has nearly the smallest maximum S/N achievable which degrades the detection sensitivity.

Zoom In Zoom Out Reset image size

A useful measure of the performance of a search is the sensitivity depth, first introduced in Behnke et al. (2015) and defined as

$$\begin{eqnarray}&&{{ \mathcal D }}^{90 \% }(f)=\displaystyle \frac{\sqrt{{S}_{h}(f)}}{{h}_{o}^{90 \% }}\,\,[1/\sqrt{\mathrm{Hz}}],\end{eqnarray} \tag{ 20 }$$

where $\sqrt{{S}_{h}(f)}$ is the noise level associated with the signal frequency f. For a network of detectors, we take S_h (f) as the harmonic mean of the single-detector one-sided power spectral densities of the data. Spectral densities of each detector are averaged over the 0.5 Hz frequency bands that the upper limits value refer to. The resulting ${{ \mathcal D }}^{90 \% }(f)$ values for the ${ \mathcal F }$-statistic and ${ \mathcal G }$-statistic analysis are shown in Figure 9.

Zoom In Zoom Out Reset image size

Outliers. Outliers are obtained by selecting the loudest point in each millihertz of the studied frequency range. In fact, following the discussion in Section 4.2, two sets of outliers are selected: those with frequency in the natural grid and those with frequency belonging to the grid of shifted bins. The two sets of outliers are shown in the top and bottom panels of Figure 10, respectively.

Zoom In Zoom Out Reset image size

It is clear from the figures that high values of the detection statistic correspond to several integer frequencies, as found also by the ${ \mathcal F }$/${ \mathcal G }$-statistic method. They are due to disturbances of instrumental origin which mainly affect the Hanford detector. ²⁹⁶ In addition, to ensure computational accuracy, the search was initially performed slightly enlarging the physical ranges of the parameters (i.e., 1.38 ≤ A ≤ 1.58 and −0.01 ≤ B ≤ 0.205). To obtain more significant results, three sets of outliers are excluded:

• 1.  
  outliers closer than 0.01 Hz to integer frequencies, in order to remove the impact of these disturbances on the noise distribution;
• 2.  
  outliers outside the physical ranges 1.39 ≤ A ≤ 1.57 and 0 ≤ B ≤ 0.195;
• 3.  
  outliers whose detection statistic in one segment is greater than the sum of the statistics on all the other four segments. ²⁹⁷ This veto limits the effect of temporary noise disturbances.

The outliers that remain after imposing these selection criteria are shown in Figure 11 for the natural frequency grid and for the shifted one. Outlier significance is estimated by means of the p-value. First, an overall p-value threshold, p_thr, is computed such that p_thr × N_p = 0.01, where N_p is the total number of points in the explored parameter space. Then, for each outlier, the detection statistic noise distribution has been built considering the corresponding 1 Hz band, after the exclusion of the 1 mHz sub-band to which the outlier belongs, plus the preceding and following 1 Hz bands. On this distribution the threshold on the detection statistic corresponding to p_thr has been computed and used to establish if that outlier is highly significant. The detection statistic threshold is shown in Figure 11 as a continuous red line and is not surpassed by any of the outliers. As before, two sets of outliers are plotted: those selected in the natural frequency grid (top panel) and those selected in the shifted grid (bottom panel).

Zoom In Zoom Out Reset image size

The detection statistic values distribution across the explored frequency range is not uniform due to the presence of noise artefacts (see, e.g., the triple peak in the range 89–90 Hz) and to the globally non-flat detector noise curve, which decreases as the frequency increases.

Figures 12 and 13 show the outlier distribution in the parameter space of $(f,\dot{f})$ and (A, B), respectively. These representations are helpful to visualize the regions defined from Equations (16) and (17), which are a rectangle for the physical parameters (A, B) and a parallelogram for $(f,\dot{f})$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As none of the outliers is loud enough to overcome the threshold, upper limits on the strain of the signal are computed.

Upper limits. Upper limits at 90% confidence level are computed every 0.25 Hz by means of software injections based on the procedure already used in standard narrowband searches (see, e.g., Aasi et al. 2015b; Abbott et al. 2017, 2019). Specifically, for each 0.25 Hz band, sets of 240 signals, with a fixed amplitude and randomly chosen parameters in their allowed range of variation, are injected in O3 data. The signal frequency and spin-down rate are shifted at each glitch epoch in such a way as to simulate glitches that occurred in the real signal we are searching for. These data are analyzed in exactly the same way as in the real search and the number of detected signals (i.e., those producing a value of the detection statistic above the detection threshold) are counted. By repeating the procedure for different amplitudes, the value corresponding to a detection efficiency of 90% is obtained through a linear interpolation among the two amplitudes corresponding to a detection probability just below and above 90%, respectively. An example of this procedure for the band 92–92.5 Hz is shown in Figure 14.

Zoom In Zoom Out Reset image size

The upper limits obtained using this method are shown in Figure 15. The lower blue curve is obtained assuming a uniform prior on the polarization angle ψ and on the cosine of the star's spin axis inclination angle ι. On the other hand, the upper red curve assumes a "restricted" prior on ψ and ι based on X-ray observations of the PWN, as discussed in Section 3.2. Specifically, we use for ψ a Gaussian distribution centered at 2.2864 rad with width of 0.0384 rad and a double Gaussian distribution for ι centered at 1.522 rad and at 1.620 rad, both with standard deviation 0.016 rad.

Zoom In Zoom Out Reset image size

The uncertainty on the upper limit is given by the combination (in quadrature) of the statistical error due to the finite number of injections (∼6%), the interpolation error in computing the amplitude corresponding to a detection efficiency of 90% (less than ∼1%), and the data calibration uncertainty (∼7%) (see Section 2.1).

We now discuss astrophysical consequences of our results. To understand if we are probing a realistic portion of parameter space and constraining theoretical models for r-mode driven spin-down of pulsar J0537, we compare our upper limits on the amplitude to the spin-down upper limit, i.e., the GW amplitude that would be needed to explain the entire measured spin-down of the star:

$$\begin{eqnarray}&&{h}_{{\rm{sd}}}=\displaystyle \frac{1}{d}\sqrt{\displaystyle \frac{10G}{{c}^{3}}{I}_{{zz}}\displaystyle \frac{\nu | \dot{\nu }| }{{f}^{2}}}\approx 4.0\times {10}^{-26}\,\left(\displaystyle \frac{90\,{\rm{Hz}}}{f}\right),\end{eqnarray} \tag{ 21 }$$

where d is distance to the source and we use I_zz = 10⁴⁵ g cm² and values appropriate for J0537 to obtain the second equality. The spin axis of the star is taken to be in the z direction, and I_zz is the moment of inertia. Note that in our case the spin-down limit is not only an upper limit, but it is exactly the predicted amplitude of the signal in the theoretical picture proposed by Andersson et al. (2018) to explain the observed inter-glitch braking index of n = 7. There is, however, an uncertainty associated with Equation (21), as the relation between f, ν, and I_zz depends on the unknown EoS of dense matter. In order to quantify this uncertainty, we make use of the results of Idrisy et al. (2015), who find that to a good approximation the r-mode (and thus GW) frequency f can be expressed as a function of compactness ${ \mathcal C }={GM}/{c}^{2}R$, independently of the exact details of the EoS, and neglecting terms of order ${\left(\nu /{\nu }_{K}\right)}^{2}$ (which is a good approximation for the rotation frequency ν = 62 Hz of J0537). The amplitude thus depends on only one parameter, which we can take to be the mass of the star M or I_zz , depending on which is more convenient.

In Figure 16, we show the physically plausible ranges for mass M and moment of inertia I_zz for the GW frequencies of our searches. Note that these ranges make no assumption as to whether r-modes actually do account for the observed spin-down; they simply reflect the range of possible values of M and I_zz consistent with a given compactness, given our lack of knowledge of the true EoS of dense matter. The upper limits on M and I_zz are physically well motivated and correspond to a causally limited EoS in the core of the star, attached to a crustal EoS at lower densities (Haskell et al. 2018), i.e., to a hard equation of state. The lower limits are observationally motivated and set by the softest EoSs considered by Idrisy et al. (2015), which are conservatively chosen as to be consistent with the 99.7% confidence limit of the observed mass of 1.97M_⊙ of PSR J1614−2230 (Demorest et al. 2010; note that the mass measurement has since been refined; see Arzoumanian et al. 2018).

Zoom In Zoom Out Reset image size

Having established the physically plausible range for I_zz , we calculate the range of h_sd from Equation (21) and compare it to the sensitivity upper limits of our searches. This is shown in Figure 17, where for simplicity we plot only a range between the upper limit set by the stiffest EoS (the causally limited EoS with crust) and the lower limits set by our softest EoS, WFF1. We see that our searches are probing a physically significant portion of parameter space for all frequencies. At frequencies above f ≳ 90 Hz, our searches are starting to probe below the limits of the r-mode-driven spin-down scenario.

Zoom In Zoom Out Reset image size

Another physical quantity of interest is the r-mode amplitude α that leads to a GW signal with an amplitude h₀, defined as (Owen 2010):

$$\begin{eqnarray}&&\alpha =\sqrt{\displaystyle \frac{5}{8\pi }}\displaystyle \frac{{c}^{5}}{G}\displaystyle \frac{{h}_{0}}{{\left(2\pi f\right)}^{3}}\displaystyle \frac{d}{{{MR}}^{3}\tilde{J}}\approx 0.017{\left(\displaystyle \frac{90\mathrm{Hz}}{f}\right)}^{3}\left(\displaystyle \frac{{h}_{0}}{{10}^{-26}}\right),\end{eqnarray} \tag{ 22 }$$

where for the second equality we assumed typical values of M = 1.4 M_⊙ and R = 12 km and used the distance to J0537. $\tilde{J}$ is the dimensionless canonical angular momentum of the r-mode, which depends only weakly on the EoS, and we fix to $\tilde{J}=0.0164$, the value calculated by Owen et al. (1998) for an n = 1 polytrope. We can use the fit of Idrisy et al. (2015) to express the frequency f in terms of compactness and reduce to a single free parameter that encodes the EoS dependence. In Figure 18, we plot upper limits on the amplitude α for our two limiting EoSs, the causally limited EoS with crust and the WFF1 EoS. We see that in both cases the sensitivity of the searches is close to or below the spin-down limit, especially for the stiffest EoS, and is starting to constrain the range of values predicted by theoretical models such as those of Alford & Schwenzer (2014) and Andersson et al. (2018). We may also visualize the constraints on the r-mode-driven spin-down scenario in terms of a constraint on the mass of the neutron star, by using the EoS-independent relations between moment of inertia and compactness of Breu & Rezzolla (2016), together with the fits of Idrisy et al. (2015), so that the only free parameter that encodes the EoS dependence in Equation (21) is the mass of the star. The results are plotted in Figure 19 where we see that our observational results allow for the r-mode-driven spin-down scenario for soft EoSs and lower-mass neutron stars. Note, however, that for our most sensitive search, the 5-vector search with uniform priors, most of the parameter space is actually excluded, and only the softest EoSs and neutron star masses M ≲ 0.9 M_⊙ are consistent with a GW driven spin-down.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is also useful to compare the upper limits on α presented above with theoretical estimates of the so-called saturation amplitude, i.e., the amplitude at which an unstable r-mode stops growing due to non-linear hydrodynamical couplings with other modes. The value of the saturation amplitude is highly uncertain and could be even orders of magnitude smaller that the upper limits presented above (see, e.g., Arras et al. 2003; Brink et al. 2004; Bondarescu et al. 2009). In this sense, the results presented here are consistent with such theoretical estimates.