Gravitational Wave Searches for Post-Merger Remnants of GW170817 and GW190425

For a rotating solid body, the strain amplitude of radiated gravitational waves is [24]

where $I_{zz}$ is the moment of inertia of the body about its axis of rotation, $\epsilon$ is its ellipticity, $f$ the frequency of the emitted gravitational waves (assumed here to be twice that of the body’s rotational frequency), $D$ is the distance to the source from the detector, $G$ is the gravitational constant, and $c$ is the speed of light.

Continuous gravitational wave search methods typically use a matched filtering process, where a signal model is compared to data in order to calculate a detection statistic. The gravitational-wave signal $s$ at a detector at time $t$ depends on the characteristic strain of the incoming gravitational wave [Eq. (1)], and the antenna pattern of the detector, encoded in four functions $h_{i},i=1,\dots,4$ [25, 24]. The signal $s$ is then given by

The parameters $\mathcal{A}_{i}$ are known as the canonical amplitudes of the gravitational wave. They are functions of the parameters $\phi_{0},\psi,h_{0}$ and $\cos\iota$ [26]. The parameters $\phi_{0}$ and $\psi$ are, respectively, the initial phase of the gravitational wave at reference time $t=t_{\text{ref}}$ and its polarisation angle. The final parameter $\iota$ is the angle between the rotational axis of the source and its sky position vector $\vec{n}$. Finally, the vector of phase parameters $\vec{\lambda}$ is composed of the sky position, frequency, and derivatives of frequency in time of the gravitational wave.

The likelihood ratio $\Lambda$ is the ratio of the probability that a signal $s(t,\mathcal{A})$ is present within noisy data against the probability that only noise is present. It is defined by [24]

where $x$ is a continuous representation of the detector data. The scalar product $(\cdot|\cdot)$ is defined as

where $T$ is the length of data used, and $\mathcal{S}_{h}$ is the single-sided spectral density of the detector noise (here assumed, for simplicity, to be constant in time).

The $\mathcal{F}$-statistic is the maximisation of Eq. (3) over the canonical amplitudes $\mathcal{A}_{i}$, and is commonly used in continuous gravitational wave searches [24]. We observe that Eq. (3) is linear in the parameters $\mathcal{A}_{i}$. This allows for the $\mathcal{F}$-statistic to be analytically maximised in these parameters:

We then wish to find the optimal phase parameters $\vec{\lambda}$ which maximise $2\mathcal{F}$.

We compute $2\mathcal{F}$ for a set of phase parameters $\{\vec{\lambda}\}$. This set is the template bank, and each $\vec{\lambda}$ is a template. If a signal is present in the data with parameters $\vec{\lambda}_{S}$, it is unlikely that they will coincide exactly with any given $\vec{\lambda}$ in the template bank. Any recovered signal from some $\vec{\lambda}$ within the template bank will therefore have some mismatch between it and the signal parameters $\vec{\lambda}_{S}$. The mismatch between $\vec{\lambda}$ and $\vec{\lambda}_{S}$ is defined using the signal-to-noise ratio $\rho(\mathcal{A},\vec{\lambda}_{S},\vec{\lambda})$, and is given as [25]

If a particular template $\vec{\lambda}$ lies close to the signal parameters $\vec{\lambda}_{S}$ such that $\Delta\vec{\lambda}=\vec{\lambda}_{S}-\vec{\lambda}$ is small, a second order Taylor expansion of Eq. (6) leads us to the parameter space metric $\mathbf{g}$ [27, 28, 25]:

where $\cdot^{T}$ represents the matrix transpose. We can see that Eq. (8) defines $\mu$ as the magnitude of the vector $\Delta\vec{\lambda}$, as measured by $\mathbf{g}$, and one can therefore treat the mismatch as a geometric distance between $\vec{\lambda}$ and $\vec{\lambda}_{S}$ within the parameter space.

A convenient approximation to $\mathbf{g}$ is the phase metric, $\mathbf{g}_{\phi}$. The phase metric relies only on the parameters $\vec{\lambda}$, and is defined as [29, 25]

The $\partial_{i}$ are partial derivatives with respect to the phase parameters which make up $\vec{\lambda}$. The $\langle\cdot\rangle$ are time averages defined as

The function $\phi(t,\vec{\lambda})$ describes the phase evolution of the gravitational wave signal. The phase is typically given as [24]

where $\vec{r}$ is the position vector of the detector with respect to the Solar System Barycentre (SSB), $\vec{n}$ is the unit vector pointing from the SSB to the source, $f_{\text{max}}$ is the maximum frequency of the gravitational wave frequency over the search band, and the superscript $(l)$ represents a time derivative. The phase model can also be written as a time integral of a gravitational wave frequency model, $f_{\text{GW}}(t,\vec{\lambda})$:

The frequency model $f_{\text{GW}}$ is typically chosen to be a first or second-order Taylor expansion [e.g. 30, 31, 32]. Taylor expansions are good signal models for the close to monochromatic signals that continuous wave search methods typically target. Additionally, Taylor expansion signal models are linear in their parameters $\vec{\lambda}$ by construction. Signal models with this property minimise the number of templates needed to cover the parameter space, and hence the computational cost of carrying out a search. Linear signal models achieve this by facilitating the implementation of lattice-based algorithms which minimise the size of the template bank while still fully covering the parameter space [33, 34].

In this section we have given an overview of the techniques used to search for continuous gravitational waves in the case of a single detector. These methods can be generalised to cases where multiple detectors are used; specifically, the definition of the $\mathcal{F}$-statistic provided here is extended to scenarios involving more than one detector as well as different noise contributions in [35]. Searches which use the multi-detector $\mathcal{F}$-statistic have improved sensitivity as they provide more data for which signal power is accumulated. The single detector $\mathcal{F}$-statistic is useful in vetoing signal candidates, as seen in Section VI.2. In the searches presented here, we use both multi-detector and single detector $\mathcal{F}$-statistics.

Four scenarios are typically considered for the evolution of the remnant of a BNS coalescence event. The remnant may: i) immediately collapse into a black hole; ii) form a hypermassive neutron star before collapsing into a black hole; iii) form a supramassive neutron star before collapsing into a black hole; or iv) form a stable long-lived neutron star [7]. Each of these scenarios is expected to be accompanied by a different gravitational wave signal.

Hypermassive and supramassive neutron stars are born with masses above what a non-rotating neutron star is able to support without collapsing into a black hole. The threshold mass to form a stable non-rotating neutron star is unknown, due to the uncertainty surrounding the neutron star equation of state. It is generally agreed, however, that this threshold mass exceeds $2~{}\text{M}_{\odot}$, as electromagnetic observations have found pulsars with masses exceeding this limit [36, 37]. Some equations of state predict maximum neutron star masses as large as $2.8~{}\text{M}_{\odot}$ [38, 39, 40]. Hypermassive neutron stars are expected to support their weight through thermal pressure and differential rotation, and have expected lifetimes on the order of milliseconds [41, 42, 43]. Supramassive neutron stars are less massive than hypermassive neutron stars, support their weight through rotation alone, and are expected to live anywhere from seconds to hours [7], to possibly up to two years [44].

Remnant neutron stars are expected to be born with large ellipticities and spinning at high frequencies [40], making them strong candidates for gravitational wave follow-up searches. Hypermassive neutron stars may have extreme deformations, with possible ellipticities as large as 0.87 [45]. The estimated ellipticities of supramassive neutron stars are affected by their assumed lifespan, and range from $10^{-7}$ to $10^{-4}$ [40]. Remnant neutron stars which are stable have estimated initial ellipticities ranging from $10^{-7}$–$10^{-2}$ provided they possess a large ($10^{15}$ G) magnetic field [46, 47, 48]; other estimates of maximum ellipticities range between $10^{-8}$–$10^{-6}$ [48, 49, 50, 51].

The gravitational wave event GW170817 was the first detection of a BNS coalescence. This event coincided with the detection of a gamma ray burst (GRB). Electromagnetic observations following the initial gravitational wave detection suggest that a short-lived neutron star was present immediately after the merger [52]. Initial gravitational wave follow-up of GW170817 focused on short- to long-transient length signals of $<$1 s to $<$500 s [12]. The searches carried out in [12] made use of the STAMP [53] and coherent WaveBurst (cWB) [54] algorithms. These search methods make use of unmodelled techniques, which are typically used to search for signals with unknown waveforms or large parameter spaces. Across the 1–4 kHz frequency band which [12] covered, peak search upper limits for 50% search confidence of $2.1\times 10^{-22}~{}\text{Hz}^{-1/2}$ and $5.9\times 10^{-22}~{}\text{Hz}^{-1/2}$ were achieved for signals of 1 s and 500 s duration respectively.

Subsequent searches looked for longer-lived signals for a GW170817 remnant, ranging from 3 hours to 8.5 days [13]. Four independent algorithms were used: STAMP [53]; unmodelled Hidden Markov model tracking in conjunction with the Viterbi algorithm [55, 56, 57]; and adaptations of the modelled Hough transform method [58] known as FrequencyHough [59, 60] and Adaptive Transient Hough [61]. These algorithms searched for a potential gravitational wave signal covering frequency bands from 300–4000 Hz; no detection was claimed. The searches were sensitive to a BNS remnant signal at a distance of 1 Mpc, and therefore not capable of detecting a post-merger signal from GW170817.

A search carried out in [14] used BayesWAVE [62], a Bayesian inference analysis to achieve upper limits on strain amplitude. No claim of detection was made. A search using a machine learning method was performed on one week of data in [15]. Sensitivity upper limit estimates were consistent with those achieved in [13].

At the time of GW190425, only two detectors were online. This resulted in poor sky localisation of the source, and despite electromagnetic follow-up efforts, no counterpart was found [63, 64]. The total system mass for GW190425 sits 5 standard deviations above the expected galactic distribution for neutron stars [11]. The remnant mass was estimated at $3.4~{}\text{M}_{\odot}$. If a neutron star formed after the coalescence, at such a high mass it is likely to have promptly collapsed into a black hole. A search for a gravitational wave signal following GW190425 may then shed light on the remnant object of this event. The detection of a gravitational wave signal from this event would indicate that massive neutron stars are possible and give further insight into the neutron star equation of state. Alternatively, it may confirm that neutron stars of these masses are unstable and again inform the equation of state.

The common parameters of the searches are summarised in Table 1; the setup was studied in detail in [23]. Search parameters specific to each source are given in Table 2. The searches take place over a frequency band of 100–2000 Hz on the 1800 s of data following the GW170817 and GW190425 events. We choose 1800 s of data for computational cost constraints and because a signal duration of this length following GW170817 has not previously been explored. The value $k_{\text{max}}$ is determined using the values of $I_{zz}$ and $\epsilon$ substituted into Eq. (34) of [23]. We use fiducial values for $I_{zz}$ and optimistic accepted values of $\epsilon$. The value $k_{\text{min}}$ is chosen to be 10% of $k_{\text{max}}$. The value of $S=2$ has been chosen as a value of $S=3$ increases the computational cost of the search with no expected improvements to search sensitivity upper limits. A single piecewise segment is used since for the 1800 s of data multiple segments are not required, however, suggestions for multiple segment searches are made in  [23]. The value of $\mu_{\text{max}}=0.2$ is chosen as a commonly used value for continuous wave searches.

The sky position of GW170817 is taken as that of its host galaxy, NGC 4993, known precisely from its electromagnetic counterpart [4]. Despite two detectors being online at the time of GW190425, only the L1 detector was able to observe the event due to its greater sensitivity. Its sky localisation was therefore poor, and no electromagnetic counterpart was detected. For this search, we take the sky position of GW190425 to be its most likely location according to the sky map available from GWOSC [68]. This choice is justified, as follows.

Unlike searches for binary black hole/neutron star systems, where the sky position is determined by triangulation of time delays between detectors, continuous wave search methods rely on the sky-position-dependent Doppler modulation of the signal phase due to the Earth’s spin and orbit. When the time-span of the data being analysed is long ($\gg$ days), the motion of the Earth causes an appreciable Doppler modulation of the signal phase, and therefore it is important to know the sky position of the source to high accuracy. In this search we use 1800 s of data, however, which means that the phase of a potential signal will only have minimal Doppler modulation. Moreover, over such a short time-span, the sinusoidal Doppler phase modulation, when Taylor expanded in time, is well-approximated by linear and quadratic terms, and is therefore degenerate with the terms which arise from the frequency $f$ and spin-down $f^{(1)}$ phase evolution, respectively:

Any error between the chosen sky position used and the true sky position of GW190425 should therefore not affect the detectability of any signal, except that such a signal may be recovered with errors in the frequency and spin-down parameters.

To carry out the search, the frequency range of 100–2000 Hz is split into four sub-bands; 100–550 Hz, 550–1500 Hz, 1500–1990 Hz and 1990–2000 Hz. Each of these sub-bands is then further divided into smaller frequency ranges to make up individual search jobs. Using these sub-bands is necessary in order to request the appropriate computational resources and for minimising the time to carry out the search. The density of templates throughout the parameter space changes with $f_{\text{max}}$, as shown in Fig. 3. Consequently, the search setup also changes with frequency. At frequencies below 550 Hz, the parameter space is extremely narrow. This leads to a significant portion of the parameter space, at its boundaries, being unaccounted for by templates. To remedy this, template padding is applied at the boundaries of the parameter space [23]. Conversely, the density of templates progressively increases in the frequency bands 550–1500 Hz, 1500–1990 Hz and 1990–2000 Hz.

The density of templates per unit frequency is determined by finding small frequency ranges across the 100–2000 Hz band which give template banks of sizes of approximately $10^{6}$. These small frequency ranges were chosen to have their upper frequency limit at 100 Hz intervals beginning at 100 Hz. The lower frequency limit was then found by manual selection until a value which gave a template bank size of $10^{6}$ was found. The density of templates was then determined by dividing the size of each of these small template banks by the width of its associated frequency band. Linear interpolation between these results then achieved an estimate of the number of templates per unit frequency across the full 100–2000 Hz range [23]. Figure 3 presents this template density.

The rapid increase in template density with frequency seen in Figure 3 is due to the GTE. Each time derivative of the GTE introduces a factor of $f_{0}^{n-1}$. Each derivative of the GTE is then proportional to $f_{0}$, and hence $f_{\text{max}}$, by; GTE $\propto f_{\text{max}}$, GTE^′ $\propto f_{\text{max}}^{n}$ and GTE^′′ $\propto f_{\text{max}}^{2n-1}$. As these form the boundaries of the parameter space, Eqs.(17)–(20), the size of the parameter space then scales strongly with $f_{\text{max}}$. This increase in parameter space size, even for a fixed frequency band, results in the strong frequency scaling of template density seen in Fig. 3.

Using Fig. 3a, each of the frequency ranges 100–550 Hz, 550–1500 Hz, 1500–1990 Hz and 1990–2000 Hz is divided into smaller ranges for each search job, each with an equal number of templates. To get the frequency range of each search job, the curve given in Fig. 3 is integrated from the start of the frequency sub-band (e.g. 550 Hz for the 550–1500 Hz band) up to a frequency which gives the desired number of templates for the given job. This sets the frequency range of the first search job. By then repeating this process from the final frequency of the previous job, a frequency range for each job is calculated, with each job containing approximately the same number of templates. An example of how the frequency range has been divided for the 550–1500 Hz band is shown in Fig. 3b. The number of templates and estimated computational cost of each job in each frequency sub-band is shown in Table 3.

The total number of templates is calculated by integrating the curve given in Fig. 3. Once known, an acceptable computational cost for each job was chosen. The search was performed on the OzSTAR supercomputing cluster, for which it was known that approximately $2.5\times 10^{4}$ templates per second could be computed [23]. Using this number, the computational cost of a search job is estimated from its frequency range and template bank size.

One might consider why in Table 3 we have split the 655 search jobs unevenly across four frequency sub-bands rather than use 655 search jobs all of equal computational cost. The search jobs have been split this way due to the higher density of templates at large frequencies, as demonstrated in Fig. 3a. To have 655 search jobs each of equal computational cost, one search job would cover a frequency band of 100–1242 Hz, calculated using data from Fig. 3a. Such a wide band is likely to contain instrumental lines which will saturate the loudest recorded $2\mathcal{F}$s. This problem is alleviated by purposefully splitting jobs into smaller frequency bands. Furthermore, at high frequencies to split jobs equally by computational cost, the search bands will be forced to become extremely narrow. This narrowing would force that padding be required, which would unnecessarily increase the computational cost of carrying out these jobs. The total number of CPU hours used to carry out the search can then be further reduced by using fewer jobs at high frequencies which cover wider frequency bands. Furthermore, search jobs cannot be divided to have equal frequency as this would make the computational cost of jobs at high frequencies impractical.

In Table 3 the number of templates, and corresponding computational cost, per individual job differs between each frequency sub-band. The frequency band 100–550 Hz must be separated from the other sub-bands due to template bank padding being required at frequencies below 550 Hz [23]. As the density of templates is lower in this region, the search jobs in this band are considerably cheaper in computational cost than those in the other frequency bands. The 100–550 Hz band is still split across 5 search jobs, despite their lower computational cost. This is to reduce the impact of instrumental lines saturating the loudest recorded $2\mathcal{F}$s. This is discussed in Section V.3. The remaining sub-bands are separated from one another due to the changing shape of the parameter space. If they were treated as a single sub-band, such that each search job from this sub-band had equal computational cost, the frequency range of each job would decrease as higher frequencies are reached. For the high frequency jobs, their frequency bands would narrow sufficiently such that they would require padding. This would then increase the computational cost of the search. Furthermore, the additional padding templates used by each search job would overlap in frequency, leading to double counting of templates at similar frequencies. This would again unnecessarily increase the computational cost of the search.

For each job, the largest $1,000$ multi-detector $\mathcal{F}$-statistics were recorded, alongside their associated templates, and the corresponding individual detector $\mathcal{F}$-statistics. With 655 search jobs, a total of 655,000 templates and their corresponding $\mathcal{F}$-statistics were recorded.

At low frequencies, the frequency bands of each search job are wide, covering many tens of Hz. Frequency bands of these widths are likely to encompass numerous instrumental artefacts, or lines, arising from noise sources in the detectors. If a strong instrumental line occurs within a search job band, templates with frequencies coinciding with the line will likely have very large $\mathcal{F}$-statistics. Given that a list of only $1000$ $2\mathcal{F}$s are stored from each job, it is possible that the templates coincident with the instrumental line will saturate this list; no templates at other frequencies in the search band will then have a $2\mathcal{F}$ value large enough to be recorded in the top 1000 templates, and any gravitational wave signal may be missed.

To remedy this behaviour, additional search jobs were run to cover the frequency ranges where no $2\mathcal{F}$s were recorded. To cover these frequency ranges, 7 additional searches jobs for GW170817 and 42 search jobs for GW190425 were used. For brevity we will continue to state that 655 search jobs and 655,000 $2\mathcal{F}$s were recorded despite these additional jobs. By including results from these additional jobs we achieve $2\mathcal{F}$ values across the entire frequency range. These results are shown in Fig. 4.

Figure 4 presents the top 1,000 candidates from each search. Through visual inspection of Fig. 4, groups of large-$2\mathcal{F}$ candidates are identified; the central frequencies of these groups are then cross-examined with known line lists from the second (O2) and third (O3) observing runs for GW170817 and GW190425, respectively [74, 75]. Groups of large-$2\mathcal{F}$ candidates which coincide with known lines are then vetoed. In Fig. 4 all vetoed candidates are shown in red, and the remaining potential candidates in blue. The identified lines used for vetoing are listed in Tables 6 and 7 in the Appendix.

Another candidate veto commonly used in continuous wave searches is a detector $2\mathcal{F}$ veto [76, 77, 78]. For this veto, the multi-detector $\mathcal{F}$-statistic must be greater in value than the individual detector $2\mathcal{F}$s for a candidate to pass. It is applied to ensure that a signal is present in both detectors. If a candidate has a large $2\mathcal{F}$ in only one detector, it is likely to originate from a local source of noise.

All candidates presented in Fig. 4 pass the detector $2\mathcal{F}$ veto, however. This illustrates a limitation of the detector $2\mathcal{F}$ veto when applied to the searches presented in this work. First, the 1800 s of data used in these searches is much shorter than typical continuous wave searches (which span weeks to years of data). Over such long time spans, a continuous wave signal will be Doppler modulated by the motion of the Earth. Instrumental lines, however, are not affected by the movement of the Earth, and can therefore be more readily identified. The 1800 s used in these searches, however, is not long enough to observe any appreciable Doppler modulation, and thus weak instrumental lines are less well distinguished from possible signals.

Second, we observe that, in the presence of weak lines, it is possible that the multi-detector $2\mathcal{F}$ may still be greater than the contributing single detector $2\mathcal{F}$s, and thereby a weak line may pass the detector $2\mathcal{F}$ veto. This behaviour is demonstrated in Figs. 5 and 6, where the single detector $\mathcal{F}$-statistics of the largest $2\mathcal{F}$s presented in Fig. 4 are plotted against each other. We see that one detector typically has $\mathcal{F}$-statistics which are much larger than the other. This suggests that weak lines present in that detector are contributing to the multi-detector $\mathcal{F}$-statistic, while still passing the detector $2\mathcal{F}$ veto. In the bottom left of each plot a cut-off in $2\mathcal{F}$ values can be seen. The presented $2\mathcal{F}$s are only the largest 1000 $2\mathcal{F}$s from each search job; if a greater number of $2\mathcal{F}$s were recorded, this cut-off would occur at lower values of the $\mathcal{F}$-statistic.

We wish to compare the behaviour seen in Figs. 5 and 6 with that expected from an astrophysical signal. To accomplish this, we carry out searches on simulated data with an injected signal, and plot the single detector $\mathcal{F}$-statistics of the top 1000 largest candidates. Each subplot shows the results from three smaller searches, with the same setup as that presented in Table 1, and different frequency ranges given in Table 4. The simulated data contains Gaussian noise with a power spectral density (PSD) equal to the experimental data for the same source (GW170817 or GW190425) at the same frequency. Signal phase parameters are sampled at random from within the piecewise parameter space. The signal strength, $h_{0}$, is the peak strain upper limit used to determine $h_{\text{rss}}^{50\%}$ in Sec. VII. This is $1.21\times 10^{-23}~{}\text{Hz}^{-1/2}$ for GW170817 and $7.4\times 10^{-24}~{}\text{Hz}^{-1/2}$ for GW190425.

In Fig. 5 and Fig. 6, we see that the largest single detector $2\mathcal{F}$s from the simulated searches (red dots) do not overlap with the search results. For GW170817 the simulated largest $2\mathcal{F}$s are generally equal in value, and do not favour one detector over another. For GW190425 the largest $2\mathcal{F}$s from the simulated searches are generally louder in the L1 detector; this is to be expected because this detector is more sensitive than V1. If the candidates which passed the detector $2\mathcal{F}$ veto (blue dots) were from an astrophysical signal, we would expect their distributions of single detector $2\mathcal{F}$s to show similar behaviour; this is not the case, however.

In summary, the inconsistency between the single detector $2\mathcal{F}$s from the GW170817 and GW190425 searches, and from the simulated searches, suggests that the experimental data contains weak instrumental lines which are the primary cause of the candidates found by the GW170817 and GW190425 searches. In the absence of a ready-made veto capable of vetoing these lines, we instead accept them as part of the background noise in the data, with the consequence that this will likely degrade the upper limits achievable by the searches (see Sec. VII).

We now look to quantify the statistical significance of the candidates found by the searches. To do this, we use the Distromax method [79, 80]. Distromax gives the likelihood that the largest value of a detection statistic – in our case, the $\mathcal{F}$-statistic – found by a search is statistically significant. This is achieved by estimating the expected distribution of the largest value of the detection statistic, assuming that no signal is present, and comparing it to the largest value found by the search. If the largest value is consistent with the distribution, we may conclude that no statistically significant candidate was found.

Distromax operates by taking in a set of detection statistics $2\mathcal{F}$. This set of $2\mathcal{F}$ is then randomly split into $N$ batches, each containing the same number of $2\mathcal{F}$ values; this is to avoid potential correlations between $2\mathcal{F}$ values in nearby regions of parameter space. Let the maxima of each of these batches be labelled $\xi^{*}_{i}$, where $i\in[1,N]$. Let us define $\xi^{*}=\max_{i}\xi^{*}_{i}$, the loudest overall $\xi^{*}$. The probability density function for $\xi^{*}$ is given by

where $P(\xi)$ is the probability distribution function for the largest value of the detection statistic in each batch. In the case of the $\mathcal{F}$-statistic, which follows a $\chi^{2}$ distribution with 4 degrees of freedom, $P(\xi)$ converges to a Gumbel distribution:

The $\mu$ and $\sigma$ of Eq. (22) are the location and scale of the Gumbel distribution. Their values are determined by constructing a histogram of the largest $2\mathcal{F}$ from each batch, and then fitting $P(\xi;\mu,\sigma)$ to the histogram. The fitted $P(\xi;\mu,\sigma)$ is then substituted into Eq. (21) to give the probability of the largest $2\mathcal{F}$ from the entire search. A property of the Gumbel distribution (see [79] for details) is that, given $P(\xi)$ is a Gumbel distribution, $p(\xi^{*})$ is also a Gumbel distribution, with location $\mu_{*}$ and scale $\sigma_{*}$ given by

This change in variables effectively propagates the Gumbel distribution $P(\xi)$ along the $2\mathcal{F}$ axis, and for this reason we refer to it as the “propagated distribution”.

We use the distribution $p(\xi^{*})$ produced by Distromax to assess the statistical significance of candidates from the searches. If a $2\mathcal{F}$ falls within a region of $p(\xi^{*})$ which is considered likely, then this candidate is consistent with noise and can therefore be rejected as a potential gravitational wave signal. Alternatively, if a $2\mathcal{F}$ is found with a value far larger than what is predicted by $p(\xi^{*})$, then that would be considered evidence for a gravitational wave detection.

Conventionally, we would supply Distromax with the set of largest $2\mathcal{F}$s from each search job. The search setup used in this work, however, only required 655 jobs (Table 3), and that number of largest $2\mathcal{F}$s is too few samples to produce reliable results with Distromax (as discussed in [79]). To increase the number of samples, we use all of the recorded $2\mathcal{F}$ values from each search job (after applying the line veto in Sec. VI.1), instead of just the largest. Using the template density information from Fig. 3a, we re-partition the entire search frequency band into sub-bands containing approximately the same number of $2\mathcal{F}$ values. The largest $2\mathcal{F}$ from each sub-band was then selected and supplied to Distromax. In this way, the number of $2\mathcal{F}$ samples was increased from 655 to 56,192 for GW190425, and 20,596 for GW170817.

For both GW170817 and GW190425 we use a total of $N=400$ batches. To use Distromax we use the example Python notebooks provided in [80]. There is a trade off in choosing the appropriate value for $N$. A small $N$ results in fewer maxima $\xi^{*}_{i}$, however each contain a larger number of samples. These batches then are well modelled by a Gumbel distribution due to the large number of samples they contain, however with fewer $\xi^{*}_{i}$ the fitting parameters of Eq. (22) have a higher variance. On the other hand, a large value of $N$ provides a larger number of maxima $\xi^{*}_{i}$ which are not well modelled by Gumbel distributions but the fitting parameters of Eq. (22) are more precise. For a more detailed discussion of choosing the batch number $N$, see [79]. We have used a value for $N$ which has produced histograms which can be well modelled as Gumbel distributions without impacting the variance of the fitting parameters of [79] significantly.

Figure 7 plots the largest $2\mathcal{F}$s from the GW170817 and GW190425 search results presented in Fig. 4, alongside their respective probability distributions calculated by Distromax. For both searches, we see that the largest $2\mathcal{F}$ found falls within the 95% credible region predicted by Distromax. We conclude that the GW170817 and GW190425 search results are consistent with noise, and therefore no gravitational wave signal has been detected.

In the lower panel of Fig. 7, all $\mathcal{F}$-statistics used as samples by Distromax are plotted as a function of template frequency. We see that, by partitioning all $\mathcal{F}$-statistics into frequency bands of similar template bank sizes, that our sample $2\mathcal{F}$s are dominated by templates at higher frequencies. This is expected, as the density of templates per unit frequency increases with frequency (Fig. 3).