Search for neutrino counterparts to the gravitational wave sources from LIGO/Virgo O3 run with the ANTARES detector

Since 2015 the LIGO and Virgo interferometers have detected gravitational waves from almost one hundred coalescences of compact objects (black holes and neutron stars). This article presents the results of a search performed with data from the ANTARES telescope to identify neutrino counterparts to the gravitational wave sources detected during the third LIGO/Virgo observing run and reported in the catalogues GWTC-2, GWTC-2.1, and GWTC-3. This search is sensitive to all-sky neutrinos of all flavours and of energies > 100 GeV, thanks to the inclusion of both track-like events (mainly induced by νμ charged-current interactions) and shower-like events (induced by other interaction types). Neutrinos are selected if they are detected within ± 500 s from the GW merger and with a reconstructed direction compatible with its sky localisation. No significant excess is found for any of the 80 analysed GW events, and upper limits on the neutrino emission are derived. Using the information from the GW catalogues and assuming isotropic emission, upper limits on the total energy E tot,ν emitted as neutrinos of all flavours and on the ratio fν = E tot,ν /E GW between neutrino and GW emissions are also computed. Finally, a stacked analysis of all the 72 binary black hole mergers (respectively the 7 neutron star-black hole merger candidates) has been performed to constrain the typical neutrino emission within this population, leading to the limits: E tot,ν < 4.0 × 1053 erg and fν < 0.15 (respectively, E tot,ν < 3.2 × 10^53 erg and fν < 0.88) for E -2 spectrum and isotropic emission. Other assumptions including softer spectra and non-isotropic scenarios have also been tested.


Introduction
Since the first detection of gravitational waves (GWs) from compact binary mergers in 2015 [1], GW interferometers have opened a new window on the Universe, complementary to the ones already being explored with other cosmic messengers (cosmic rays, photons, neutrinos). These capabilities have already allowed the association of the GW signal GW170817 emitted by the merger of a binary neutron star system with the emission of a short gamma-ray burst (GRB), GRB 170817A, detected by Fermi and INTEGRAL in gamma rays, and the related afterglow across a wider range of the electromagnetic spectrum [2].
Neutrinos are also expected to be emitted from the relativistic outflows that characterise such mergers: see e.g., [3] for binary neutron star (BNS) mergers, [4] for neutron star-black hole (NSBH) mergers, and [5] for binary black hole (BBH) mergers. Previous searches performed with ANTARES [6], Baikal-GVD [7], IceCube [8], and Super-Kamiokande [9] have not been able to identify an excess of neutrinos and upper limits have been reported.
This article presents an updated search using the latest GW catalogues covering detections in 2019-2020 and the ANTARES data from the same period. Besides the follow-up of individual events, performed in the same way as in previous publications, first population studies are also presented by carrying out a stacking analysis for binary mergers of the same nature and taking advantage of the extensive catalogue reported by the GW community.

The GW catalogues
This paper focuses on GW sources detected during the third observing run (O3) of the LIGO and Virgo detectors and reported in the following official LIGO/Virgo catalogues: • GWTC-2 [10]: this catalogue reports detections made during the first half of O3 (April -September 2019). It contains 39 candidates, including 1 BNS, 2 NSBH, and 36 BBH events.
• GWTC-2.1 [11]: this is an update of GWTC-2, with 8 additional events not reported in the previous catalogue and that have a high probability of astrophysical origin. The catalogue also includes a larger selection of ∼ 1.2 k events with a false alarm rate of less than 2 per day but lower astrophysical probability, which are not considered in this analysis.
• GWTC-3 [12]: this catalogue covers the second half of O3 (November 2019 -March 2020) and contains 35 objects, including 4 NSBH candidates. Seven marginal candidates that do not fully satisfy the criteria for an astrophysical origin are also reported, of which only GW200105 162426 is kept as it has also been reported independently as a plausible NSBH candidate [13].
For each of these objects, the LIGO-Virgo collaboration provides a FITS file [14] with the timing of the merger t GW and the constraints on the source direction Ω as a skymap P(Ω), as well as posterior samples containing all the correlations between source direction Ω, luminosity distance estimate D L , and other source parameters such as the masses of the two merging objects m 1,2 (with the convention m 1 > m 2 ), the energy radiated in gravitational waves E GW (defined as the difference between the estimated mass of the final object and the sum of the masses of the initial objects), and the inclination between the total angular momentum and the line-of-sight θ jn . The classification among the different categories is made based on the mass estimates: BNS if m 2 < m 1 < 3 M , NSBH if m 2 < 3 M < m 1 , BBH otherwise. When considering the different catalogues, 83 objects are selected, including 1 BNS and 7 NSBH candidates.

The ANTARES telescope
The ANTARES neutrino telescope [15], located in the depths of the Mediterranean Sea, offshore from Toulon (France), has been operating in its final configuration between May 2008 and February 2022. It was composed of an array of 885 photomultiplier tubes (PMTs) enclosed in pressure-proof glass spheres, arranged in triplets over 12 vertical lines, spaced by ∼ 70 m and anchored at a depth of ∼ 2475 m.
The PMTs detect the Cherenkov light induced by relativistic charged particles originating in the interaction of a neutrino with matter surrounding the detector; the space and time pattern of PMT signals (or hits) allows the neutrino properties (direction and energy) to be inferred. Charged-current interactions of muon (anti-)neutrinos are characterised by the presence of a long muon track; the typical angular resolution for such events is < 1°for E ν > 100 GeV [16]. Other types of neutrino interactions only produce hadronic (and, in the case of ν e charged-current interactions, electromagnetic) showers, with a more localised light deposit and compact topology. Despite the shorter lever arm with respect to the muon tracks, ANTARES still achieved a median angular resolution of a few degrees for the neutrino direction [17]. In the following, the two event categories described above are referred to as the track and shower samples respectively. ANTARES data were organised in consecutive runs of at most twelve hours. The present analysis uses data from 2019 and 2020 to identify neutrino counterparts to the GW emissions described in section 1.1. A larger dataset including data from January 2018 to December 2020 is used for background estimation.
Several searches for neutrino counterparts to GW events with ANTARES have been carried out in the past. These studies were limited to event-by-event follow-ups and reported null results. See [18] for GW170104, [19] for GW170817, and [6] for other O2 events.

Analysis method
This search focuses on the selection of neutrino events in a time window of 1000 s centred on the GW emission time t GW , as motivated in [20], and in the region R 90 containing 90% of the source localisation probability, as built directly from the GW 2D skymaps P(Ω). As the reconstructed event direction does not match the true neutrino direction perfectly due to the scattering angle and the finite detector resolution, this region of interest (RoI) is extended by an angle α to account for these effects, meaning that an event with direction x would be selected if min d∈R90 (arccos (x · d)) ≤ α. This extended angle α is a free parameter of the analysis, whose choice is based on the optimisation method described below.
The ANTARES events are divided into four categories according to whether they are classified as tracks or showers and whether their reconstructed direction is upgoing or downgoing, each case based on a specific selection and optimisation procedure, as described in sections 2.1 and 2.2.
The reconstructed data are largely dominated by atmospheric muons. A selection is then applied to reduce it to an expected number of events B = 2.7 × 10 −3 , such that the detection of one event would correspond to a 3σ excess (3σ condition). This condition is separately set for each category. For a given GW, the background expectation is estimated using the dataset from 2018 to 2020. Only runs with similar data-taking conditions as the ones during the ANTARES run r GW overlapping with the GW time, characterised by the mean burst fraction 1 , are selected. As illustrated in the left panel of Figure 1, this procedure is found to allow for a proper estimate of the background while ensuring a better characterisation of the tails of the distribution as compared to a statistically-limited estimation using only the data from the run r GW .
A dedicated Monte Carlo (MC) simulation [22] has been produced to be used in this optimisation process. Simulations are done on a run-by-run basis, where each run has a specific simulation to reproduce the particular environmental and detector conditions during this run. For GW events with precise sky localisation (R 90 region smaller than 3000 deg 2 ), additional neutrinos generated solely within this region have also been produced to accurately estimate the corresponding detector acceptance for an . More details about the detector acceptance for a given neutrino spectrum are in [23].
The final cuts are optimised to ensure that the expected number of selected background events after all cuts is fulfilling the 3σ condition defined above.

Track event selection
The track selection procedure is adapted from the one presented in [18]. The upgoing and downgoing (respectively with reconstructed incoming direction below or above the horizon) event selections are different, as the latter is more likely to be contaminated by the atmospheric muon background and needs extra care.
For upgoing tracks, a cut on the track reconstruction quality parameter, Λ [24], is applied. The values of α and of the cut on Λ are optimised by ensuring the 3σ condition described above, as well as maximising the signal acceptance in the hypothesis of an E −2 spectrum. For downgoing events, the procedure is similar except that an additional cut on the number of hits employed for the reconstruction is also applied, as illustrated in the right panel of Figure 1.

Shower event selection
The selection steps for neutrino interactions yielding showering events are similar to the ones presented in [6]. Events must be contained within the detector and must not be classified as a track by the selection described in the previous section. The discrimination between neutrinos and atmospheric muons is achieved thanks to an extended likelihood ratio L µ defined by comparing the neutrino and muon hypotheses for each hit associated with the shower on the basis of its deposited charge, timing, and distance to the reconstructed shower position [17]. Another parameter is used to further reduce the background contamination: for upgoing events, this parameter is defined from a Random Decision Forest (RDF) classifier [25] while the downgoing selection exploits the number of hits used in the event fitting.
As for the track selection, the values of the cuts on these parameters and on the extension of the RoI α are optimised to ensure the 3σ condition and to maximise the acceptance, as illustrated in Figure 2.

Detector systematics
Several systematic effects may affect the detector performance, hence the obtained constraints on the neutrino emission. Three sources of uncertainty, found to be the dominant effects as already described e.g. in [6], are taken into account and evaluated independently for the four event categories (track/showers upgoing/downgoing):   The right panel shows the average acceptance, assuming an E −2 neutrino spectrum, as a function of the cut on the number of hits and the value of α. For each bin, the cut on Λ is fixed to the value that satisfies exactly the 3σ condition. The optimal working point, maximising the acceptance, is indicated by the red cross.    Figure 1. The right panel shows the average acceptance, assuming an E −2 neutrino spectrum, as a function of the cut on the RDF classifier and the value of α. For each bin, the cut on Lµ is fixed to the value that satisfies exactly the 3σ condition. The optimal working point, maximising the acceptance, is indicated by the red cross.
• The first one is related to the uncertainty on the PMT photon detection efficiency and on the water absorption length; the related uncertainties have been re-evaluated by varying these two quantities within a typical interval of ±10% in dedicated MC simulations and estimating the overall impact on the signal acceptance.
• The second source is linked to the capability of the run-by-run MC simulations to properly reproduce data conditions; the related error is estimated by comparing the variability of event rates between data and simulations.
• The last effect is the combined statistical and systematic uncertainty on the background expectation. The related uncertainty is obtained by varying the list of similar runs employed for its estimation.
The total uncertainties on the acceptance related to the first two sources are 18%, 14%, 21%, and 19% respectively for the upgoing tracks, downgoing tracks, upgoing showers, and downgoing showers. The overall uncertainty on the background is about 20% for all event categories.

Statistical analysis
For each GW event, the number of observed neutrino candidates in time and spatial coincidence in each category can be converted into a significance of the observation using Poisson statistics. In the absence of any excess of neutrino events with respect to the background expectation, upper limits on the neutrino emission are calculated. In the following section, several assumptions are made: (i) The source localisation is supposed to be within R 90 (as this is the region for which the selection has been optimised). Therefore, the GW posterior samples are restricted to those with Ω ∈ R 90 and the final constraints neglect the chances for the actual source to be localised in the rest of the sky.
(ii) There is equipartition between the neutrino flavours at Earth due to the averaging of oscillations over astronomical distances [26], starting from ν e : ν µ : ν τ = 1 : 2 : 0 at production. This allows reporting limits on the all-flavour neutrino emission.
(iii) The neutrino energy spectrum is described by a single power law dN/dE = φ · (E/GeV) −γ where φ is expressed in units of GeV −1 cm −2 and γ is the spectral index. The nominal case is γ = 2 (E −2 spectrum).
The ANTARES acceptance is estimated using the MC simulations described in section 2. It depends on the shape of the assumed neutrino spectrum (characterised by γ), on the source direction Ω, and on the event category c. It is averaged over the neutrino flavours and can be decomposed into a normalisation factor and a direction-dependent component:

Constraints on the neutrino flux
This section presents the limits on the overall flux normalisation φ obtained for an all-flavour emission with γ = 2. The cut-and-count analysis described in the previous sections corresponds to a Poisson likelihood where N (c) (resp. B (c) ) is the observed (resp. background-expected) number of events in each category, and the product is performed over the set of four event categories (c ∈ C, C = {upgoing tracks, downgoing tracks, upgoing showers, downgoing showers}). Given the selection optimisation presented in section 2, the value of B (c) is fixed to 2.7 × 10 −3 , independently for all categories. A Bayesian method is employed to obtain constraints on the neutrino flux normalisation φ starting from this likelihood. A flat prior on φ is employed, the systematic uncertainties described in section 2.3 are encoded in Gaussian priors on B (c) and on a (c) γ with the standard deviations corresponding to the uncertainties reported there, and the GW skymap P(Ω) is used as a prior on Ω.
The obtained posterior probability is then marginalised over all the nuisance parameters (background, acceptance, direction) by using MC integration techniques: toy samples (t.s.) are generated with values of B (c) and a (c) γ following the priors, and the posterior samples from GW catalogues can be used directly for the sampling of Ω. The marginalised posterior probability distribution is computed as where C is a normalisation constant that can be determined numerically by ensuring ∞ 0 P (φ)dφ = 1. The 90% upper limit φ 90 is finally obtained by solving φ90 0 P (φ)dφ = 0.90.

Constraints on the total energy
Similarly to the incoming neutrino flux on Earth, one may also constrain the total energy E tot,ν emitted in neutrinos, correcting for the source distance. This can be done under a specific assumption on the spatial distribution of the neutrino emission around the source, e.g., either isotropic or collimated into a jet. One may also consider the ratio between the total energy emitted in neutrinos and the energy E GW radiated in GW: f ν = E tot,ν /E GW .
Isotropic emission. In the case of a source emitting isotropically, the total energy emitted in neutrinos can be computed as where D L is the source luminosity distance, and E min , E max are the integration bounds. For this analysis, these bounds are fixed to E min = 5 GeV and E max = 10 8 GeV, which is the typical range where the emission is expected in most of the models (e.g., [27]). Since the total energy E iso tot,ν is proportional to the flux normalisation φ, the likelihood from Equation (3.1) can be rewritten in terms of E iso tot,ν instead of φ. The marginalised posterior probability and the upper limits are obtained similarly, where the luminosity distance is also extracted from GW posterior samples to be used in MC integration toy samples and a flat prior on E iso tot,ν is assumed. A similar rewriting is possible to obtain the limits on f iso ν = E iso tot,ν /E GW . This is a relevant parameter if it is assumed that the total energy emitted in neutrinos scales with the GW emission.
Non-isotropic scenarios. One may also consider the total energy E tot,ν for a given non-isotropic model, and the corresponding likelihoods and posteriors can be written using the relevant GW parameters for the marginalisation. For non-isotropic emission, a simple von Mises [28] jet model is considered: × exp(cos θ/ω 2 ), (3.4) where ω characterises the jet opening and θ is the angle with respect to the jet central direction.
In the following, the jet is assumed to be collinear with the total angular momentum of the merger, such that θ = θ jn from the GW data release. The total energy emitted in neutrinos for a power-law spectrum is then where the last term is the correction corresponding to the jet visibility from Earth (in the isotropic case, this term would be 2, such that the Equation (3.3) is retrieved). For a given jet opening ω, the limits on E tot,ν or the corresponding f ν may be derived as described in the isotropic case, including θ jn variable in MC integration toy samples.

Stacking analysis
The GW catalogues published by LIGO/Virgo may contain populations of sources with similar neutrino emissions, which could be more efficiently constrained by performing stacking analyses. A first version of such an analysis is presented here with two categories of GW events: the 72 BBH mergers on one side and the 7 NSBH mergers on the other. While differing in its precise implementation in this study, the method follows the ideas initially presented in [29] and already applied in [9] with Super-Kamiokande data.
The stacking approach assumes that all objects in the selected population have the same emission, either in terms of total energy E tot,ν or of f ν . In the non-isotropic scenarios, GW sources in a given population may have different jet inclinations but the shape of the jet (the parameter ω for the von Mises model) is considered to be the same for all the objects.
Assuming a flat prior on the signal parameters and considering all observations as independent, the posterior distribution for a given population S may be written as where C is a normalisation constant, X is the signal parameter to be constrained (E tot,ν or f ν ), Y may denote potential common parameters (such as ω), and Z i represents the parameters which may differ from one GW event to another within S. Upper limits on X can then be obtained with the same method as before.
The stacking approach is particularly interesting to constrain the non-isotropic models, as is detailed in the next section.

Results
A total of 80 GW events out of the 83 observed during O3, including all candidates involving at least one neutron star, can be associated with exploitable ANTARES data.
For each follow-up, the neutrino event selection is optimised according to the procedure described in section 2 and the final number of selected events in time and spatial coincidence with the GW event in each category is extracted. No event has been selected for any of the follow-ups, which is fully compatible with the expected background i∈GWs c∈C B (c) i ∼ 0.82. Therefore, only upper limits on the neutrino flux and other related quantities are reported in the following. Table 1 displays the 90% upper limits on the incoming all-flavour neutrino emission assuming an E −2 spectrum, along with GW information, for the events initially reported in the GWTC-2 catalogue. Similarly, the results for the events in GWTC-2.1 and GWTC-3 are presented in Table 2 and Table 3, respectively. The all-flavour flux limits assuming an E −2 spectrum stands mostly between 4 and 60 GeV cm −2 . Figure 3 shows the limits on the total energy emitted in all-flavour neutrinos and on f iso ν = E iso tot,ν /E GW , assuming isotropic emission, for events from the three catalogues. The limits range from 10 54 to 10 59 erg and are mostly following the expected D 2 L trend. Population studies are performed by stacking (a) all 72 BBH events, and (b) the 7 NSBH candidates. Figure 4 shows the stacked upper limits on the typical total energy emitted in neutrinos (or the ratio f ν ) within these two categories, considering jetted emission with the von Mises model (as a function of the jet opening ω) as well as for isotropic emission, and assuming E −γ spectra with γ = 2.0, 2.5, 3.0.
For individual follow-ups, as θ jn is not strongly constrained by GW detections, von Mises emission constraints are mainly dominated by geometrical effects: when the jet opening angle is narrow, there is a very small chance that the Earth is aligned with the jet, while a wide opening leads to a larger spread in total energy, which is hence more difficult to constrain. The limits on the total energy as a function of the jet opening are thus expected to show a typical parabola shape corresponding to a trade-off between these two effects, as illustrated by the envelope of individual limits shown in Figure 4. While the individual results for jetted emission are of little interest due to these limitations, Table 1. Summary of GWTC-2 follow-up results. The first four columns summarise the relevant GW information [10], including the most probable merger type, the estimated median luminosity distance DL, and the size of the region R90 containing 90% of the localisation probability, while the three last ones show the 90% upper limits on the all-flavour neutrino emission assuming an E −2 spectrum (this work), in terms of the flux normalisation E 2 dN/dE, the total isotropic energy E iso tot,ν , and the ratio f iso ν .
GW name Type D L R 90 area Upper limits on neutrino emission  Figure 3. 90% upper limits on the total energy E iso tot,ν emitted in neutrinos of all flavours (left) and on f iso ν = E iso tot,ν /EGW (right) as a function of the source luminosity distance, assuming an E −2 spectrum and isotropic emission. The horizontal bars indicate the 5 − 95% range of the luminosity distance estimate, and the markers/colours correspond to the different source categories.
Finally, it is worth noticing that the statistical power of the BBH stacking already allows the exploration of reasonable values for f ν (in the context of the assumed simple power-law spectrum), while all individual limits on this parameter are above 1, corresponding to a neutrino emission higher than the GW one.
Supplemental material with all numerical results for individual follow-ups and from the stacking analysis, for the various jet scenarios and assumed spectral indices, can be provided under request.

Discussion and conclusions
The search for neutrino counterparts to GW signals detected during the O3 run with the ANTARES detector yields no significant excess. This null result is used to extract upper limits on the neutrino emission, both in terms of the flux normalisation and of the total energy emitted in neutrinos, E iso tot,ν , in particular for the standard case with an E −2 spectrum and isotropic emission.
The ANTARES analysis presented in this article benefits from the different event topologies in the detector, leading to all-flavour limits for neutrino energies above 100 GeV. The constraints are also 10    derived as well. In the absence of the determination of any jet direction from the GW data or any electromagnetic observations, these results are only indirectly deducted from the stacking analysis.
In future observation campaigns, any such measurement, even for one or a few sources, would allow for deriving direct constraints, as well as considering more complex jet models. Though the ANTARES detector has been decommissioned in 2022, the KM3NeT experiment [32] has taken over the observation of the neutrino sky from the depths of the Mediterranean Sea. Data from the first few lines of KM3NeT/ORCA has already been used to perform an initial search for neutrino counterparts to the O3 events [33]. With more lines being deployed up to ∼ 2026, KM3NeT/ARCA and ORCA are expected to provide complementary coverage of the transient sky with respect to IceCube, increasing the chances to get significant individual detections and creating new opportunities for joint analyses.
The next GW observing period, O4, is expected to start in 2023, including the KAGRA detector in addition to LIGO and Virgo [34]. With an expected number of detections of the order of one hundred per year, population studies will become even more relevant.
Lastly, combinations with ZTF [35] or Vera C. Rubin [36] surveys and Target of Opportunity programs may help pinpoint candidate host galaxies for observed binary mergers, and allow more detailed analyses using the electromagnetic information and the environment characterisation as inputs for the prediction and constraint of the neutrino emission. More generally, there are still only a few models on the processes leading to neutrino production during binary merger events in the literature [3-5, 27, 30], and further theoretical developments in the coming years may bring new light to past, present, and future observations.