Search for gravitational-neutrino correlations on ground-based detectors

The problem of joint data processing from ground-based gravitational and neutrino detectors is considered in order to increase the detection efficiency of collapsing objects in the Galaxy. The development of the"neutrino - gravitational correlation"algorithm is carried out within the framework of the theory of optimal filtration as applied to the well-known OGRAN and BUST facilities located at the BNO INR RAS. The experience of analyzing neutrino and gravitational data obtained during the outburst of supernova SN1987A is used. Sequential steps of the algorithm are presented, formulas for estimating the statistical efficiency of a two-channel recorder are obtained.


Introduction
After the registration by ground-based gravitational-wave (GW) interferometers LIGO and Virgo of bursts of gravitational radiation from the merger of relativistic binaries: -black holes [1,2] and neutron stars [3], as well as after the discovery of a relativistic object simultaneously generating gravitational and electromagnetic (gamma) signals, with subsequent optical afterglow [3], one can confidently speak about the emergence of a new GW channel of astrophysical information. This new channel, effectively supplementing the previously developed methods for registering cosmic rays and particle fluxes, contributes to the development of a new heuristic approach to astronomical observations, called "multi-messenger astronomy" or in Russian transcription -"multi-channel astronomy".
To date, there are no operating large GW interferometers in Russia, although a project to create such an installation in Siberia is known [4,5]. The only gravitational detector in the kilohertz frequency range is the combined optoacoustic antenna OGRAN, developed jointly by the Russian Academy of Sciences and Moscow State University. Details of the structure, principle of operation, and technical parameters of this setup can be found in [6,7]. The latest upgraded version, along with the current experimental characteristics, is presented in [8]. The proposed article is devoted to the analysis of the data processing technique of this antenna, for which it is first useful to briefly recall its principles and design scheme.
OGRAN is an installation located in the underground premises of the Baksan Neutrino Observatory of the Institute of Nuclear Research, Russian Academy of Sciences, designed to search for collapsing stars in the Galaxy together with the BUST neutrino telescope. Both instruments are sensitive enough to register collapses in our Galaxy as rare events with an average rate of 0.03 events per year.
Structurally, OGRAN ( Figure 1) consists of two arms (or "feedback loops"). The first one, the "detector arm", contains a large aluminum bar -an acoustic resonator of longitudinal vibrations with a mass M ≈ 2 t and a length L ≈ 2 m. The second, "discriminator arm", is equipped with a small glass-ceramic bar (M ≈ 30 kg, L ≈ 40 cm). Both bars have a cylindrical shape with drilled central axial channels for a light guide. Pairs of mirrors are mechanically attached to the ends of the bars, forming Fabry-Perot (FP) standards. When illuminated from a common (tunable) laser with a wavelength of λ ≈ 1 µm, these FP standards form the optical degree of freedom of the OGRAN antenna, with which GW can interact directly, i.e. acting on the light field in the FP. Antenna tuning in the operating mode binds the laser frequency to the selected optical FP resonance in the detector arm. In this case, information about the acoustic oscillations of the detector (large bar) is encoded in the frequency of the optical pump. In the arm of the discriminator (small bar), tuning to optical resonance is carried out by moving one of the FP mirrors due to feedback circuits at low (quasistatic) frequencies. As a result, frequency demodulation is performed and the output signal reproduces the acoustic vibrations of the gravitational detector. In general, the design of the OGRAN antenna is a differential circuit known as the "comparator of optical standards" [9,10].
The upgraded OGRAN version has the following optical parameters. Optical pump power in each of the arms is ∼ 200 mW , FP sharpness (finesse)in the detector arm is 30,000 and 70,000 in the discriminator arm. Interference contrast 0.3 ÷ 0.7. The sensitivity curve of the antenna is shown in Figure 2. It can be seen that for the current OGRAN version those metric variations h are registered which cause deformation perturbations of the length of the detector It should be noted that this sensitivity is not enough to detect signals from collapsing stars at a distance of about 8 kpc (the center of the Galaxy). An increase in sensitivity by two orders of magnitude without changing the design of the antenna is, in principle, possible by forcing its parameters: sharpness (finesse) of FP resonators in the arms up to 300000 and optical pump power up to 20 W . In this case, the sensitivity will reach the level of 10 −22 Hz −1/2 in a band of about 100 Hz [6,7,8].
We also note that OGRAN is located in the underground premises of the BNO, not far from the BUST neutrino detector (BUST) which plays an important role in the search for and detection of collapsing objects within the Galaxy. The joint coordinated use of these tools increases the likelihood of success in this task. The goal of this study is to develop an efficient algorithm for joint analysis (processing) of data from both detectors -neutrino and gravitational.

2
Mathematical model of the gravitational detector and GW signal The acoustic part of the OGRAN detector is the same as that of resonant solidstate detectors [12,13], i.e. appears to be the main antisymmetric mode of longitudinal oscillations of a cylindrical bar (with length L much greater than its radius r), which makes it possible to model the detector as a one-dimensional oscillator with two parameters: normal mode period τ 0 and relaxation time τ r . As a generalized model of the GW signal (external force impact on the de- tector), a short wave train (or radio pulse) with amplitude F s and duration τ s ≪ τ r filled with oscillations close to its resonant frequency ω 0 , is considered. In accordance with the GR formulas, in the weak field approximation, the action of a monochromatic gravitational wave (GW) with frequency ω on a one-dimensional detector is equivalent to a quasi-resonant swinging force with an amplitude F s = 0.5mω 2 Lh 0 , where m is the equivalent mass of the detector, h 0 is the dimensionless amplitude of the metric variations of the transferred GW associated with the gravitational energy flux density in the wave [14,15]. Since due to the equivalence principle, the action of the GW on the detector is determined by the transferred acceleration, regardless of its mass, it is natural to introduce the notation f s = F s /m. In the mathematical formulation of the problem of calculating the effect of a GW train on a gravitational detector, we use the apparatus of complex envelopes [16]. Thus we represent f s (t) = Re[f s (t)e jω0t ]. For a wave train, the form of the complex envelope depending on time is where the A is a wave train amplitude andτ -wave train duration.
To describe the thermal (Brownian) oscillations of the detector, we introduce f T (t) -the random Nyquist force, i.e. white noise with a spectral intensity N 0 = 4k B Θmδ, which depends on the absolute temperature Θ, and the value of the losses of the detector given by the attenuation factor δ = τ −1 r , associated with the quality factor of the detector Q = ω 0 /(2δ) ≫ 1 (the so-called "Langevin approach" to description of Brownian motion [17] ). As a result, the response spectrum of the acoustic detector (the resonant unit of the OGRAN setup) is obtained by passing the GW train f s (t) through the frequency transfer function of the oscillator The time form of the response s(t) is determined by the impulse response g(t) through the convolution transformation where the sign ↔ means the transition from the temporal form of the response to the spectral form, which reduces to the product of the spectra of the considered functions. We immediately present the form of the spectrum f s (jω) of the complex envelopef s (t) of the signal wave train, which is used in further analysis In addition to the signal f s and thermal f T effects on the detector, one should take into account the "measuring noise" (or sensor noise), which for OGRAN is associated with the optoelectronic conversion of the acoustic oscillations of the bar into a measured electrical signal. Usually fluctuations of measuring instruments are represented as an additive interference ξ(t) of the white noise type in a limited frequency band at a spectral intensity N ξ , the estimate of which is taken from the experiment, < ξ(t)ξ(t + τ ) >= N ξ δ(τ ). Thus, roughly (without fine details) the output of an opto-acoustic gravitational detector can be written as As indicated above, the signal part is represented by the convolution s(t) = f s (t) * g(t), the fluctuation part is formed as the sum n(t) = η(t)+ξ(t), including thermal noise η(t) = f T (t) * g(t) with spectral density N T (ω) = N 0 | g (jω) | 2 and measuring noise ξ(t) . Considering that noises η(t) and ξ(t) are independent, the spectral intensity of the output of the detector decomposes into the sum of the spectral densities of thermal and additive noise N (ω) = N T (ω) + N ξ .

Optimal filtering algorithm
The procedure for optimal filtering of the detector output of an optoacoustic detector y(t) = s(t) + n(t) lies in passing it through a linear optimal filter with a transfer function where * means complex conjugation, t 0 is the time of maximum response, const means some arbitrary constant. At the filter output, we get a narrow-band process z(t) = Re[z(t)e jω0t ] with a complex envelopez(t).
To explain the physical meaning of the filtering K opt (jω), it is convenient to split it into two successive transformations K opt (jω) = K 0 (jω)K 1 (jω), thus separating the stage of the so-called "Wiener filtering" K 0 (jω) [15,16] from the subsequent matched filter K 1 (jω) [16]. This is a consequence of the fact that in the problem under consideration there are two types of noise: white additive measurement noise and thermal noise, which has a finite spectrum in a given limited frequency band. As a result, we arrive at a representation in which the transfer function of the Wiener filter for the total noise with an effective bandwidth Ω ∼ ω 0 N 0 /4m 2 N ξ 1/2 reads as and the transfer function of the matched filter against the background of thermal noise reads as With small losses in the vicinity of the resonant frequency ω 0 ± ω, ω ≪ ω 0 , the formula for K 1 can be simplified in accordance with the approximation Then, collecting formulas (4,8,9), we get a clearer form of the transfer function of the matched filter K 1 (ω); Considering thatf s (jω) ≈ (1/2)f s (j (ω 0 + ω)) we come to the form which in time domain corresponds to the operation of the difference link for the variable x (t) ←→ x (jω) at the output of the Wiener filter. As a result, the procedure for optimal processing of the OGRAN detector output y(t) is the chain of transformations The ratio between the spectra x (jω) and z (jω) near the resonant frequency ω 0 is determined by the condition 2x (ω + ω 0 ) The role of the Wiener filter lies in the frequency cutoff of the background of the measuring white noise ξ(t). The subsequent matched filtering K 1 (jω) is determined by the spectrum of thermal fluctuations N 0 |g (jω)| 2 and, under the condition of short GW signals τ ≪ τ r , is reduced to the final procedure of the difference link x(t) − x(t − t) = z(t) for the complex envelope. Note that in practice such an envelope can be constructed in terms of the quadrature components of the output z(t).
Summarizing the procedure for optimal processing of the OGRAN output, we write out in explicit form the expression for the transfer function of the optimal filter K opt (jω) in the vicinity of the resonant frequency ω 0 and the formula for calculating the signal-to-noise ratio (SNR) at its output, namely here Ω is the previously introduced denotation of the Wiener filter effective bandwidth.
It is interesting to compare SN R = q for two signal models: the wave train model adopted by us, for which the GW burst spectrum width is less than the Wiener filter bandwidth , i.e. Ω τ > 1, and a very short burst like δ (τ )an impulse for which Ω τ ≪ 1 (the model used for Italian cryogenic antennas [12,13,15]). In the first case, from (13) we obtain q 1 ≈ A 2 τ /2N 0 , while in the second case we have a much smaller value q 2 = q 1 Ω τ /2.

Envelope outliers detections
The optimal filtering procedure of the OGRAN output signal described above, in fact, leads the operator to the need to observe bursts or outliers at the output of the difference link △ x(t, τ ) = z(t) constructed taking into account the phase of the quasi-harmonic variable z(t). In the absence of a priori information about the moment of arrival of the GW signal and any information from parallel recording channels (multi-messenger astronomy), the only method for detecting external influences is the so-called Neyman-Pearson strategy (NPS) [18]. In this methodology, an indicator of the presence of an impact is the excess of the number of emissions above the average for the selected threshold level C α , which is determined by the statistics of the observed variable.
For sufficiently high threshold levels (C/σ) > 1, where σ is the standard deviation, random outliers can be considered independent and subject to the Poisson distribution, representing the probability P {k} of the number of outliers k that occurred during the observation time T . In the simplest case, the average number of outliers λ is assumed to be constant thus a flow of uniform intensity is considered. Then the "Poisson law" is written as P {k} = λ k (1/k!)e −λ , where λ =< n (C, T ) > (<> sign of statistical averaging).
At the OGRAN output (after the optimal filter) there is a narrow-band Gaussian random process z(t) = Re[z(t)e jω0t ] with the correlation function k(t) =< z(t)z(t + τ ) >= Re k (τ )e jω0τ , wherek(τ ) = σ 2 ρ(τ )) is a complex envelope with variance σ 2 =< z 2 (t) > and correlation coefficient ρ (τ ) . Accordingly, for the average number of envelope outliers | z (t)| above the threshold level C during the time (0, T ) the following formula takes place [17,18,19] < n(c, where ρ ′′ 0 = d 2 ρ(τ ) dτ 2 at τ = 0. The application of NPS to the "observed variable" (in this case, to outliers at the output of the difference link △ x(t, τ ) -at the first step, it consists in calculating the value of the threshold level C α corresponding to the selected statistical error of the 1st kind α < 1, -the so-called "probability of false alarm" in the theory of detection, or "probability of chance" in nuclear physics (the mathematical term "significance level" is also used). For a Poisson distribution with an average value (14) it is not difficult to find that the threshold corresponds to the condition for the absence of outliers P {k = 0} = 1−α. There are no outliers above such a threshold with an error α. To specify formulas (14), (15), it remains to give the results of calculating the parameters ρ ′′ 0 and σ 2 Within the framework of NPS, two approaches are possible to make a decision about an external influence (presence of a signal).
The first one (the standard approach) is tracking a random variable, the number of outliers during the observation time T , over a certain threshold level (C/σ), the height of which is limited only by the independence condition of the outliers (the absence of their correlation). Then the excess of the observed number of releases over the average value, corresponding to high reliability of the event (standard -0.95) will serve as a statistical criterion for registering an external impact. This corresponds to the usual rule of accepting the "non-null hypothesis" i.e. hypotheses of "the presence of an external perturbation" that changes its own statistics of the observed variable.
The second approach ("absolute maximum criterion") is the observation at a high threshold level given by the condition (15) of a single "excess event", i.e. the main (largest) maximum crosses the threshold level (15) at least once.
Obviously, the first approach should in principle be more sensitive to weak influences, since it allows operation at relatively low levels. However, the degree of registration reliability (confidence limit) is expected to be higher in the second case. In general, the choice of approach depends on the physics of a particular problem. To conclude this section, we present a formula for the "probability of correct detection" D when using the "absolute maximum criterion" (second approach). It is not difficult to find that with a known ratio SN R = q >> 1, the estimate of this indicator is D ≈ Φ √ q − Cα σ , where Φ (x) is the error function.

BUST and neutrino detection technique
In the BNO, the registration of neutrino signals accompanying the appearance of collapsing stars in our Galaxy is included in the program of observations of the Baksan Underground Scintillation Telescope (BUST). This large-sized detector with dimensions (17 × 17 × 11) m 3 is located in an underground laboratory protected from cosmic particle flows by a rock corresponding to 850 meters of water equivalent [20]. Structurally, the BUST consists of four horizontal and four vertical planes filled with cells (counters) of liquid scintillators. The total number of counters is 3184 with a total scintillator mass of 330 tons. Each counter is an aluminum container (0.7 × 0.7 × 0.3) m 3 filled with white spirit (C n H 2n+2 ; n ≈ 9). However, in the "collapse program" only a part of the most protected (internal) counters is used, -the number of 1200 with a mass 130 t targets. Thanks to improved protection, they have a relatively low background event count rate: f = 0.02 s −1 Most of the events that the BUST registers from the supernova (SN) explosion are inverse beta decay reactions:ν e +p −→ n+e + . (interaction of an electron antineutrino with a target proton generates a neutron and a positron). At a typical average antineutrino energy E νe = (12 − 15) M eV [21], the range of the produced positron e + will, as a rule, be contained in the volume of an individual counter. Therefore, the signal from SN will appear as a series of events, when only one counter is activated from the total number of cells on the installation ("single event") Thus, the strategy for searching for a neutrino burst at the BUST is to register a cluster (group) of single events during the time interval τ = 20 s (estimation of the duration of a neutrino burst from SN ). For SN at a distance of 10 kpc, the total energy emitted in a neutrino is 3 · 10 53 erg. With a target mass of 130 tons (three lower BUST planes), the estimate of the expected number of events from a single collapse is N ∼ 35 [20,21]. Of course, there is a random background of such events, created by the radioactivity of the environment, cosmic ray muons, false alarms of counters, etc. The background is such that it creates a cluster of k = 8 single events with a frequency of 0.138 year −1 . At the same time, no more than 2 such clusters can be expected to appear in 10 years. The rate of formation of clusters from k = 9 background events is already much lower and amounts to 7 · 10 −3 year −1 . It follows that clusters with k > 9 cannot be created by the background and, therefore, are candidates for registering the SN event.
Identification of suspicious 20-second intervals in experimental data can be done in two ways. In the first one, a sliding 20-second time interval moves in discrete steps from one single event to the next so that there is always at least one event in the cluster (at the beginning of the interval). In this case, the Poisson distribution law of events may be violated. Another processing option is possible when the event clusters do not overlap. The beginning of each time interval coincides with the end of the previous one. The first interval is chosen arbitrarily. In this case, the distribution of clusters according to the number of events is strictly Poisson, but a cluster with a higher multiplicity can be lost due to some of the events falling into a neighboring cluster. The work [20] refers to the use of both methods for mutual control. It also states that the "radius of sensitivity" of the BUST is approximately 20 kpc. This region includes about 95 % of the stars in our Galaxy. For more distant SNe, the number of single events in a cluster will be less than nine, in which case it is necessary to investigate correlations with other installations. For the period from 1980 to 2021, the net observation time for collapses at the BUST was ∼ 36 years which is the longest time of observation of the Galaxy on the same setup. During this time, not a single event of a candidate for the collapse of a stellar object (a cluster with k ≥ 9) was registered. This leads to the value of the upper limit of the average frequency of gravitational collapses in the Galaxy f col < 0.064 year −1 at the 90 % confidence level [22]. Theoretical estimates of the frequency of occurrence of galactic SNe with core collapse give a value of approximately 2-5 events per century [23].

Algorithm for searching for neutrino-gravitational correlations
Let us refine the form of the expected neutrino signal that accompanies the process of collapse of a relativistic star whose mass exceeds the critical one. According to [24,25], the most active phase of the process develops in a short time of the order of ∼ 1 s. The neutrino luminosity curve usually has a large initial peak followed by weaker peaks reflecting radiation during bounces in the process of monotonic compression [24]. There are also other multistage collapse scenarios [26,27] that predict the presence of a neutrino pulse flux at longer times, ∼ 20 s. In this regard, it is reasonable to study the algorithms for the joint registration of events in which bursts of neutrino pulses significantly correlate in time with the excitations of the gravitational detector. In fact, a similar algorithm was used in the processing of neutrino-gravity data during the SN 1987A flare [28,29,30,31,32].
The initial reference points will be considered the times t k (E) of the appearance of "neutrino signals" of the BUST in the channel (measuring mode) "search for collapsars" at those observational monitoring intervals (∆T = T 2 − T 1 ) ≈ 20 s, which fix an abnormally increased count rate " neutrino background" (see Section 4); the designation t k (E) is adopted to emphasize also the possibility of estimating the energy release of the event E). Further, at the same time intervals, the signals (emissions) of the OGRAN gravitational detector are studied.
Let us consider a separate (single) neutrino event (count) registered at the BUST at the time t k (E) with energy release E. We assume that it can be accompanied by a delayed (or advanced) gravitational burst arising at the output of the gravitational detector optimal filter matched with the GW radiation model in the form of a train of duration τ s . Then the moment of appearance of the gravitational burst is estimated as τ k = (t k (E) + τ s + δτ ), where δτ is some shift between neutrino and gravitational events, limited by |δτ | ≤ ∆τ . The range of used shifts ∆τ is set from physical considerations using astrophysical models of collapse [24,25,26,27] Let's move on to the problem of detecting the correlation between the registered packets of neutrino -gravitational signals. We use a model in which the gravitational detector output signal (after the optimal filter) is written as z(t) = ϑS(t) + n(t), S(t) = k s k (t − τ k ) on the interval T 1 ≤ t ≤ T 2 ; the symbol ϑ = (1, 0) is a formal parameter of the presence or absence of a signal. This model assumes that the signal disturbance is given by the sum of individual non-overlapping pulses. An individual impulse in a burst s k (t) can be represented as (see Section 2) s k (t) = a k Re exp (−jϕ k ) f k exp (−jw 0 t) * g(t), for 0 ≤ t ≤ τ s .
The simplified version of the model also contains the assumption that the durations of individual pulses are equal, i.e. τ sk = τ s , and the universality of delay shifts between gravitational and neutrino signals. In addition, the considered neutrino-gravitational events are considered independent. Then the logarithms of the likelihood ratios of the individual events are added. It is natural to apply the maximum likelihood ratio criterion to our problem since its mathematical apparatus is well developed for detection problems against the background of Gaussian noise.
In particular, it is known that for narrowband normal noise (the OGRAN detector), the logarithm of the likelihood ratio Λ(z) is proportional to the square of its output signal envelope R(t) [18], lnΛ k (z) ∝ R 2 [t k (E) + τ s ] -for one impulse; lnΛ(z) = k lnΛ k (z) for a pack.
The procedure for obtaining a solution (registration of the presence of a correlation) corresponds to finding the maximum of "sufficient statistics" Dthe integral likelihood ratio with a variation in the time shift between neutrino -gravitational events δτ , i.e.
here Q is the number of pulses in the packet, σ 2 is noise dispersion at the output of the optimal filter, C α is decision threshold ϑ = 1 with significance level α, the calculation of which requires knowledge of the distribution function of sufficient statistics D.
Taking into account that with a Gaussian noise background the statistics of the variable k R 2 (τ k ) /σ 2 at ϑ = 0 obeys the χ 2 distribution, for which the integral probability density has the form F 1 (x) = Γ (Q, x/2) Γ (Q), where Γ (m, x) is an incomplete gamma function; using its explicit expressions, we find the following estimates F 1 (C α ) ≈ (1 − α) ν , ν ≈ ∆τ (−ρ ′′ 0 ) 1/2 , which make it possible to calculate the threshold C α to confirm the neutrino-gravitational correlation according to the OGRAN and BUST data.

Conclusions
A priori, the development of an algorithm for estimating the neutrino-gravitational correlation seems to be an incorrect task due to completely different principles (mechanisms of action) of the corresponding detectors. In this work, this can be done by reducing the output signals of the detectors to one type of random process -the Poisson pulse stream. The question of the presence of a statistical connection between two stochastic impulse flows is solved based on the "correlation criterion". The combined observable variable ("optimal" or "sufficient statistics") is chosen as the sum of gravitational signals arising in a close time neighborhood of neutrino events. In fact, such a variable is proportional to the cross-correlation function of the two considered Poisson processes. The correlation criterion corresponds to the maximum value of the observed variable in the space of mutual temporal shift of impulse processes of different physical nature. In this approach, the solution of the problem to the maximum, in addition to calculating the probability of the presence of a correlation, also gives an estimate of the delay of the neutrino signal relative to the gravitational one (flux shift corresponding to the maximum) and, thereby, an estimate of the rest mass of the neutrino (at the GW light speed). Of course, the algorithm presented in this paper is a simplified scheme of the first approximation. It is limited by the assumptions that the considered Poisson fluxes are homogeneous (with a constant average pulse repetition rate) as well as the delay between neutrino and gravitational signals. The model in which the flow of neutrino signals has a time-varying average pulse repetition rate will apparently be more adequate to the physics of the BUST operation in the "search for collapsars" mode. However, this significantly complicates the calculation of the statistical characteristics of detection in the calculations of such a scheme. The authors are planning such a study in the next work.
The last remark is related to the degree of generality of the proposed search algorithm for neutrino gravitational correlation. In other words, can it be used for other more sensitive gravitational detectors such as LIGO/Virgo and neutrino telescopes such as Super-Kamiokande, DUNE, Hyper-Kamiokande, and JUNO. The answer, by and large, is no. Our development is directly related to a specific pair of OGRAN -BUST, since the optimal filtering of the output signals of the detectors depends on their principle of operation and design. But most importantly, the experiment under discussion -the search for collapses in the Galaxy -belongs to the category of "registration of rare events" -at best, one event in 30 years. Such observations are usually not included in the priority programs of the aforementioned detectors.
Preference is given to significantly more probable events with a higher probability of occurrence. This is achieved by taking into account a large number of distant galaxies while increasing the radius of the observation sphere and, consequently, the distances to possible sources of gravitational and neutrino signals. In this case, registration of neutrino bursts with energies of the order of 10 M eV is impossible due to the insufficient sensitivity of scintillation detectors to cover very distant objects. They are being replaced by Cherenkov-type telescopes which include the installations mentioned above. The energies of recorded bursts from cosmic catastrophes (including special types of collapses) are in the range from GeV to T eV and higher. Examples of such registrations can be found in [33] for the Ice Cube facility in Antarctica and [34] for the Antares facility in the Mediterranean Sea. Registration of low-energy events with energies on the order of tens and hundreds of M eV by the Ice Cube facility was discussed in [35].
Attempts to compare these events with data from LIGO gravitational detectors are contained in [36]. Algorithms for searching for such correlations will apparently be refined. However, a common feature of the refined versions will be the optimization of the time setting of the shift interval between the gravitational and neutrino signals.