1. 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 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
t and a length
m. The second, “discriminator arm”, is equipped with a small glass-ceramic bar (
kg,
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
μ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 (quasi-static) 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
. 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 (large bar)
in
Hz bandwidth or
in
Hz bandwidth around
kHz acoustic resonance center frequency.
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 300,000 and optical pump power up to 20 W. In this case, the sensitivity will reach the level of
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 solid-state 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
and relaxation time
.
As a generalized model of the GW signal (external force impact on the detector), a short wave train (or radio pulse) with amplitude
and duration
filled with oscillations close to its resonant frequency
, 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
, where
m is the equivalent mass of the detector,
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
. 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
. 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
—the random Nyquist force, i.e., white noise with a spectral intensity
, which depends on the absolute temperature
, and the value of the losses of the detector given by the attenuation factor
, associated with the quality factor of the detector
(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
through the frequency transfer function of the oscillator
The time form of the response
is determined by the impulse response
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
of the complex envelope
of the signal wave train, which is used in further analysis
In addition to the signal
and thermal
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
of the white noise type in a limited frequency band at a spectral intensity
, the estimate of which is taken from the experiment,
. 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 , the fluctuation part is formed as the sum , including thermal noise with spectral density and measuring noise . Considering that noises and are independent, the spectral intensity of the output of the detector decomposes into the sum of the spectral densities of thermal and additive noise .
3. Optimal Filtering Algorithm
The procedure for optimal filtering of the detector output of an optoacoustic detector
lies in passing it through a linear optimal filter with a transfer function
where
means complex conjugation,
is the time of maximum response,
means some arbitrary constant. At the filter output, we get a narrow-band process
with a complex envelope
.
To explain the physical meaning of the filtering
, it is convenient to split it into two successive transformations
, thus separating the stage of the so-called “Wiener filtering”
[
15,
16] from the subsequent matched filter
[
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
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
,
, the formula for
can be simplified in accordance with the approximation
Then, collecting Formulas (
4), (
8) and (
9), we get a clearer form of the transfer function of the matched filter
. Considering that
we come to the form
which in time domain corresponds to the operation of the difference link for the variable
at the output of the Wiener filter. As a result, the procedure for optimal processing of the OGRAN detector output
is the chain of transformations
The ratio between the spectra and near the resonant frequency is determined by the condition . The role of the Wiener filter lies in the frequency cutoff of the background of the measuring white noise . The subsequent matched filtering is determined by the spectrum of thermal fluctuations and, under the condition of short GW signals , is reduced to the final procedure of the difference link for the complex envelope. Note that in practice such an envelope can be constructed in terms of the quadrature components of the output .
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
in the vicinity of the resonant frequency
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
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.,
, and a very short burst like
—an impulse for which
(the model used for Italian cryogenic antennas [
12,
13,
15]). In the first case, from (
13) we obtain
, while in the second case, we have a much smaller value
.
4. 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
constructed taking into account the phase of the quasi-harmonic variable
. 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
, which is determined by the statistics of the observed variable.
For sufficiently high threshold levels , where is the standard deviation, random outliers can be considered independent and subject to the Poisson distribution, representing the probability 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 , where ( sign of statistical averaging).
At the OGRAN output (after the optimal filter) there is a narrow-band Gaussian random process
with the correlation function
, where
) is a complex envelope with variance
and correlation coefficient
. Accordingly, for the average number of envelope outliers
above the threshold level
C during the time
the following formula takes place [
17,
18,
19]
where
at
.
The application of NPS to the “observed variable” (in this case, to outliers at the output of the difference link
—at the first step, it consists in calculating the value of the threshold level
corresponding to the selected statistical error of the 1st kind
, 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
. There are no outliers above such a threshold with an error
. To specify Formulas (
14) and (
15), it remains to give the results of calculating the parameters
and
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 , 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—) 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
, the estimate of this indicator is
, where
is the error function.
5. 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
is located in an underground laboratory protected from cosmic particle flows by a rock corresponding to 850 m 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
filled with white spirit (
). 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:
Most of the events that the BUST registers from the supernova (SN) explosion are inverse beta decay reactions:
. (interaction of an electron antineutrino with a target proton generates a neutron and a positron). At a typical average antineutrino energy
= (12–15) MeV [
21], the range of the produced positron
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
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
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
[
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
single events with a frequency of
. 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
background events is already much lower and amounts to
. It follows that clusters with
cannot be created by the background and, therefore, are candidates for registering the SN event.
Identification of suspicious 20 s intervals in experimental data can be done in two ways. In the first one, a sliding 20 s 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
) was registered. This leads to the value of the upper limit of the average frequency of gravitational collapses in the Galaxy
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].
6. 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
of the appearance of “neutrino signals” of the BUST in the channel (measuring mode) “search for collapsars” at those observational monitoring intervals
s, which fix an abnormally increased count rate “neutrino background” (see
Section 4); the designation
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
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
. Then the moment of appearance of the gravitational burst is estimated as
, 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 us 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
,
on the interval
; the symbol
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
can be represented as (see
Section 2)
, for
.
The simplified version of the model also contains the assumption that the durations of individual pulses are equal, i.e., , 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
is proportional to the square of its output signal envelope
[
18],
—for one impulse;
for a pack.
The procedure for obtaining a solution (registration of the presence of a correlation) corresponds to finding the maximum of “sufficient statistics”
D—the integral likelihood ratio with a variation in the time shift between neutrino-gravitational events
, i.e.,
here,
is noise dispersion at the output of the optimal filter,
is decision threshold
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 at obeys the distribution, for which the integral probability density has the form , where is an incomplete gamma function; using its explicit expressions, we find the following estimates , which make it possible to calculate the threshold to confirm the neutrino-gravitational correlation according to the OGRAN and BUST data.
7. 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. However, 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 MeV 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 TeV 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 MeV 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.