Search for multi-messenger signals in NOvA coincident with LIGO/Virgo detections

Using the NOvA neutrino detectors, a broad search has been performed for any signal coincident with 28 gravitational wave events detected by the LIGO/Virgo Collaboration between September 2015 and July 2019. For all of these events, NOvA is sensitive to possible arrival of neutrinos and cosmic rays of GeV and higher energies. For five (seven) events in the NOvA Far (Near) Detector, timely public alerts from the LIGO/Virgo Collaboration allowed recording of MeV-scale events. No signal candidates were found.


I. INTRODUCTION
Recent years have seen an explosion in multi-messenger astronomy and the potential for discovery continues to increase at a rapid pace. For many years, the only extrasolar source detected by more than one messenger -defined as photons, neutrinos, gravitational waves (GW) and cosmic rays -was Supernova 1987a [1][2][3], seen in neutrinos as well as across the electromagnetic spectrum. With the advent of gravitational wave astronomy [4], the joint observation of GW170817 [5] with GRB 170817A [6,7] has been added to the list. More recently, the flaring blazar TXS 0506+056 was associated with a high-energy neutrino observed by the IceCube observatory [8].
A flux of high-energy neutrinos, of GeV-scale and higher, is expected from any compact object merger with a neutron star remnant or merger that occurs within a significant concentration of gas [9]. Additionally, any compact object remnant would emit MeV neutrinos as it cooled. The merger of two neutron stars would initially produce a hot neutron star that cools primarily via neutrino emission [10]; this hot neutron star may or may not subsequently collapse to form a black hole, but produces a neutrino flux regardless. Also, gravitational waves are expected to be emitted by core-collapse supernovae, which are known neutrino sources, provided an asymmetric collapse occurs with quadrupole or higher moments [11]. Finally, gravitational waves from unknown or exotic sources (e.g., cosmic strings [12]) may be associated with neutrino bursts. Despite these possibilities, no neutrinos have been observed to date coincident with any gravitational wave event [13][14][15][16][17][18].
These several possibilities motivate a broad search for any detectable activity in the NOvA detectors in coincidence with gravitational wave events. Although at extra-galactic distances, we do not expect a detectable flux at NOvA from compact object mergers (observed by LIGO/Virgo) or supernovae (not yet observed), it is valuable to check this hypothesis in case either a source of gravitational waves has been misidentified as extragalactic when it is not, or some observable flux production mechanism has been overlooked. As NOvA has significant sensitivity to supernova-like neutrinos, we report limits on the supernova-like neutrino flux. While we search in higher-energy channels as well, NOvA's sensitivity to the usual flux models is limited compared to other observatories. As our only likely sensitivity is to the unexpected, we do not set flux limits for higher energy signals.

II. DETECTORS
The NOvA experiment [19] consists of two detectors separated by 809 km. The detector design was optimized for the detection of ν e appearance in a ν µ beam, specifically for the discrimination between neutral current events containing a π 0 and ν e charged current events. The requirement that the radiation length be significantly longer than a detector element set the size of the scintillator cells and motivates the use of low-Z materials. The NOvA detectors have been collecting data from the Fermilab NuMI (Neutrinos at the Main Injector) beam since 2013 [20,21].
The Near Detector (ND) is located underground at Fermilab, with 22 meters water-equivalent overburden. It is designed to measure the unoscillated neutrino flux produced by Fermilab's NuMI beamline. The Far Detector (FD) is located in northern Minnesota, on the surface but slightly below grade, with a modest 3 meter waterequivalent overburden provided by 1.3 m of concrete and 16 cm of barite. The ND is relatively small, with dimensions 16 m by 4.1 m by 4.1 m and a mass of 300 t, while the FD has dimensions 60 m by 15.6 m by 15.6 m and a mass of 14 kt. The long axes point 28 • west of north; this direction is called +z, with the short axes x and y forming a right-handed coordinate system in which +x is west and +y is up. In the context of this search, the ND is a small, low-background detector as compared to the large, high-background FD.
The two detectors are functionally identical and consist of alternating vertical and horizontal planes of PVC cells [22] filled with liquid scintillator [23]. The cells are 4 cm by 6 cm and extend over the width or height of the detector. Each cell contains a single loop of wavelengthshifting fiber that extends from the readout end, down the entire length of the cell, and back to the readout end. This scheme allows for efficient light collection without the need to instrument both ends of each cell. Both ends of each fiber are coupled to a single pixel on a 32 pixel avalanche photodiode array.
The last 3 m of the ND is a muon range stack consisting of ten 10 cm thick steel plates with two scintillator planes, one horizontal and one vertical, between each steel plate. The FD has no muon range stack. With the exception of the steel plates, the detectors are 62% scintillator by mass.
Signals from each cell are continuously digitized by front-end electronics. Energy depositions over threshold are recorded for further processing. This threshold depends on position within the detectors and is typically a few MeV. Detector-wide trigger decisions are made in a farm of Linux computers. Triggers can be issued based either on the characteristics of the data or based on external signals. For instance, the "energy trigger" in each detector reads out candidate physics events if the total energy in a window of time exceeds a fixed value, and the "NuMI trigger" reads out data in a window around a timestamp received from Fermilab for each beam pulse. Data segments are available to be read out by any number of triggers independently; no trigger causes dead time for any other. The data buffer is about 30 minutes deep at the ND and 22 minutes deep at the FD, with variations caused by the number of buffer computers currently operating and the detectors' raw data rates [24].
For this analysis, the energy triggers at the ND and FD are used to collect candidate events above about 100 MeV and 50 GeV, respectively. These triggers read out for as long as a high energy burst continues, up to 20 ms. At the FD, a minimum bias "10 Hz trigger" collects either 550 µs (97% of triggers) or 500 µs (3%; occurs when the trigger lines up with a readout block boundary) at regular 100 ms intervals. This trigger is used to search for events below the energy trigger's threshold, albeit with only about a 0.55% livetime fraction.
Finally, the "LVC triggers" at both detectors receive alerts sent by the LIGO/Virgo Collaboration (LVC) [25] over the Gamma-ray Coordinates Network (GCN) each time a gravitational wave event candidate is detected. If such a trigger is received while the data are still available, 45 s of continuous data are read out, beginning 5.16 s before the LVC timestamp. The readout begins significantly before the LVC timestamp in order to capture a baseline for background subtraction, or conceivably correlated activity preceding the peak GW power. The precise time offset between the LVC timestamp and the beginning of NOvA readout is arbitrary and related to features of NOvA's data acquisition system (DAQ). Partial readouts can occur for triggers received when only some of the requested 45 s is still available; these are still analyzed (see Table I). The 45 s window is motivated partially by the length of time that a detectable flux of neutrinos is expected from a galactic supernova. An even longer readout would be better for this purpose, but 45 s was determined to be a period that could be stably recorded by the DAQ. For signals other than nearby supernovae, we do not know of a model which motivates any particular time window and which would only produce events below NOvA's energy trigger thresholds. A test 45 s trigger is issued each morning at 8:30 local time, both ensuring stability of the system and providing data for background estimates.
The LVC trigger has been active for the LIGO/Virgo Collaboration's "O3" run starting in 2019 [26]; prior to that only the energy and FD 10 Hz triggers are available. Each detector has recorded data from several other triggers (e.g., the "NuMI trigger"), but the total livetime of these other triggers was negligible and/or the characteristics of the trigger selections were unsuitable for this search.
Although not used in this analysis, we also run two triggers to collect data in the case of a galactic supernova. One responds to an alert from SNEWS [27], while the other is a self-trigger which responds to data collected by NOvA itself [28]. In the case that a gravitational wave event were caused by a nearby supernova, these triggers would collect neutrino interactions in NOvA even in the absence of an LVC trigger.

III. SEARCHES
So as not to miss any unanticipated signals, a variety of searches were performed, each designed to be as generic as possible. The energy and 10 Hz triggers were used to search for any burst within a 1000 s window centered on an LVC event timestamp in addition to the 45 s readout from the LVC triggers. We searched for bursts of (1) events selected by energy from several MeV to many TeV without regard to detailed event characteristics; (2) contained GeV-scale events; and (3) events with tracks, further broken down into several categories. First, a search was performed for events similar to those expected to be caused by O(10 MeV) supernova neutrinos. NOvA is primarily sensitive toν e through the inverse beta decay channel, with 75% of interactions expected through this channel in the no-oscillation case. We have some sensitivity to ν e via ν e 12 C → e − 12 N (5% of interactions) and any flavor through electron elastic scattering and excitation of carbon nuclei (20%).
This selection was optimized separately for the ND and FD and was designed to maximize the signal significance S/ √ B, where the simulated signal, S, used for the optimization procedure was generated using the Garching supernova flux [29] with the GENIE Monte Carlo generator [30]. The background, B, was determined from minimum bias data.
To select activity as supernova-like, first all tracks (which are mainly cosmic ray muons) and other GeVscale clusters of activity are removed, as well as any hit within 14 planes or 28 cells of such activity and within a time window extending from 2 µs before to 13 µs after the GeV-scale activity. This effectively removes any hits that were associated with the GeV-scale activity but not identified as such by the clustering algorithm, as well as removing all Michel electrons from muon decay. For highly energetic cosmic ray events, the time cut is extended from 13 µs to 200 µs, which removes neutron capture activity and spurious hits caused by electronics effects. Even with the high background level in the FD, which is on the surface, these simple cuts remove only 11% of the signal.
Further, hits near the top or sides of each detector are removed. Hits must be at least 50 (20) cells from the top of the detector in the FD (ND), 10 (4) cells from the east and west sides of the detector, and 2 (4) planes from the north and south ends. Hits of very low and very high energy are removed to eliminate noise and activity in excess of supernova-like energy. Once individual hits are selected in this manner, candidate events are formed from hit pairs consisting of one hit in a horizontal plane and one in an adjacent vertical plane. Given the z positions of the planes, the x position of the vertical cell, and the y position of the horizontal cell, the three-dimensional position of the cluster is determined. These hits must have times within 250 ns of each other after correcting the timing of each hit using the coordinate provided by the other and the propagation speed of light in the detector.
For both detectors, this selection has 20% efficiency for supernova-like inverse beta decay events within the accepted volume. The background rate is 450 Hz (0.5 Hz) at the FD (ND). While the FD background is large, the expected peak rate of selected events for a supernova in the galactic core is ∼4 kHz, which makes such a signal easily observable. Extra-galactic sources associated with gravitational wave events, however, would need to be substantially brighter in supernova-like neutrinos to be seen by NOvA.

Sub-supernova-like
To search for lower energy signals in the range of 1-10 MeV, we run two similar selections. Each is the same as the supernova-like selection above, except that any hit selected as part of a supernova-like event is removed from consideration (to create a statistically independent sample), and individual hits are selected instead of pairs. In the first of these selections, the low energy requirement for hits is lowered to just above the level associated with avalanche photodiode noise, equivalent to a few MeV. In the other, there is no low energy requirement.
Without the requirement that events include a hit in each view, three-dimensional locations of candidates cannot be reconstructed, making it impossible to know when a hit is near the end of a cell and therefore near an edge of the detector, increasing the cosmogenic background. Natural radioactivity is also selected -the NOvA design made no attempt at radiopurity. Both of these considerations increase the background rates dramatically. The background rate at the FD (ND) is 42 MHz (190 kHz) for the selection without a low energy cut and 550 kHz (38 kHz) with most electronics noise excluded by the low energy cut.
In the case of gravitational wave events for which we received an LVC trigger and read out a period of continuous data, usually 45 s, these three MeV-scale searches (supernova-like, sub-supernova-like with a low energy cut, and sub-supernova-like without a low energy cut) were run separately both for the continuous LVCtriggered data and for 1000 s of 10 Hz trigger data in the FD. While these data streams are not entirely disjoint, only 0.55% of the LVC-triggered data are also present in the 10 Hz trigger data, so we neglect the overlap.

High-energy Far Detector events
At the FD, the energy trigger is used to search for any excess of events depositing 50 GeV or higher. Besides examining the trigger rate, six selections are made to select higher-energy events with two general topologies. The energy trigger selects events in which the majority of the energy is deposited promptly, and also integrates the deposited energy up to 20 ms to select periods of time with a large total activity. Three selections are made of events in which the energy appears within a single 50 µs time window, with the requirements of an energy deposition of at least 200 GeV, 2 TeV and 20 TeV for the first, second and third selection, respectively. A second set of three selections allows the energy to arrive over a longer period of time -up to 20 ms -with total energy depositions of at least 400 GeV, 4 TeV and 40 TeV.

B. Contained events
FD data were examined for any contained activity. Such activity would be indicative of neutrino interactions in the GeV to tens-of-GeV range, however in this search no neutrino-like requirements were imposed on the event topology. To be considered a contained event, all hits of a cluster must be at least 130 cm from the bottom, east and west faces of the detector, at least 75 cm from the north and south faces, and at least 280 cm from the top. The cluster must have at least 10 hits. The number of planes between the northernmost and southernmost hits in horizontal planes must be at least nine, with the same requirement made of the vertical planes. These requirements eliminate most cosmic ray activity. Furthermore, for both views the occupancy within the smallest rectangular box that surrounds all the hits must be at least 2%, to prevent the selection of uncorrelated low energy activity, and no more than 10%, to eliminate classes of electronics noise which cause spurious hits in many adjacent channels.
The efficiency of this selection depends on assumptions about the origin of a potential signal. Most notably, physically larger events will be selected with lower efficiency because they are more often near the edges. Some loss of efficiency also occurs because background cosmic ray activity can overlap in time and space with signal events, causing them to appear uncontained. For few-GeV neutrino-like events, this effect reduces the efficiency by only a few percent. Similar considerations apply to the fully and partially-contained track selections below.

C. Track selections
In each detector, the time distribution of tracks is examined with tracks selected in nine ways, using all combinations of three track topologies and three pointing requirements. The three topologies are: fully contained tracks, tracks which start or stop in the detector, and an inclusive selection of any kind of track. In each case counts are made (1) without regard to pointing; (2) with the requirement of pointing to the LVC 90% allowed region convolved with a 1.3 • resolution; and (3) the same, but with a 16 • resolution. Convolution was performed with Healpix [31], which provides functions for analysis of binned data on a sphere. The background rate of events in the most inclusive category -any type of track with any pointing -is 110 kHz (36 Hz) at the FD (ND). When multiple tracks are detected in coincidence, they are counted as a single event with the least-contained track setting the category. This procedure ensures that the background rate closely follows a Poisson distribution. The 1.3 • convolution represents an estimate of NOvA's track pointing resolution and so this selection generically represents hypotheses that would cause charged particles that appear in NOvA to point directly back at the gravitational wave source, or nearly so. The most likely scenario would be detection of a part of a high energy neutrino interaction in the atmosphere or rock surrounding the detector. The 16 • convolution is meant to select secondary charged particle tracks resulting from lower energy interactions, and was set to represent the approximate range of reconstructed muon angles resulting from 10 GeV ν µ interactions in or near the detector. Despite these motivations, the selections are meant to be as generic as reasonably possible and do not assume any particular interaction model.
In addition, in the FD only, upward-going muon tracks are selected. Because of how light propagates in the long cells in NOvA's design, nanosecond-level timing is not available for individual hits and therefore track direction is difficult to determine for short tracks. At the FD, track direction can be determined for tracks over 8 m by fitting the timing distribution under the upward and downward hypotheses. This method is used to select upward-going muons, a potential signal of ν µ interactions either in the detector or the rock beneath it. As with the other track selections, this is repeated using the allowed sky regions convolved with 1.3 • and 16 • resolutions.

IV. NUMI BEAM VETO
The ND is exposed to Fermilab's NuMI neutrino beam which provides 10 µs-long pulses of ν µ orν µ with a mean energy of 2.7 GeV every 1.3 s to 1.4 s during beam operations. Several neutrino interactions are typically recorded in each beam pulse. All hits recorded from the beginning of each pulse to 3 ms afterwards are eliminated from this analysis. This time interval is sufficient to allow all significant effects of the neutrino interactions to end. The longest-lived such effect is caused by neutrons produced in the surrounding rock, which can thermalize in the rock, then travel through the air in the detector hall for several meters at ∼ 225 m/s before arriving at the detector and being captured with a characteristic time of 50 µs. While beam interactions also produce radioactive isotopes, including 12 B and 12 N, the rate of their decays is small compared to the background rate, so they do not motivate a longer beam exclusion window.
The NuMI interaction rate at the FD of O(1) interaction per day is negligible. Likewise, the ND is also exposed to Fermilab's Booster Neutrino Beam, but is far enough off axis to yield an event rate of about 1.5 per day, and this rate is also neglected. Table I shows a summary of NOvA data collected for each of the gravitational wave events reported by the LIGO/Virgo Collaboration from GW150914, the first detection, in September 2015, through S190707q, the last event analyzed in this report. As a result of receiving prompt triggers from LVC, low-energy data were recorded with good efficiency for five (seven) events in the FD (ND). Alerts must be received within 10 minutes to ensure a full readout. For S190426c and S190513bm, the trigger was received sufficiently late (25 and 28 minutes, respectively) that only partial readout was possible at the ND. For S190706ai, full ND readout was possible given the 18 minute trigger latency, but only partial FD readout. A prompt alert was sent for S190521r, but our connection to GCN was down at the time. In the remaining cases, data including only low-energy events were no longer available when the trigger arrived.

V. DATA SET
High-energy data, along with low-energy data with 0.55% livetime, were taken with the full 1000 s window around the gravitational wave timestamp for all other events with the exceptions of (1) S190408an, for which both detectors were down, (2) GW151012, for which the FD was down, and (3) GW150914, for which the FD was taking data, but suffering DAQ instability (we do not use such data here).

VI. ANALYSIS
For each selection described above, we searched for any excess in 1-second bins during the 1000 s analysis window, or the 45 s trigger window for LVC-triggered data streams. The bin width is intended to be similar to the duration of the initial pulse of neutrinos from a supernova as well as to that of a short gamma ray burst, such as that detected along with GW170817 [6,7], while not being finely tuned to any particular model. Different strategies are used to determine background level, depending on the characteristics of each sample. These are described below and summarized in Table II. For several selections, the background level is many hertz. In these cases, we measured the background directly as the mean rate in the analysis window, assuming that no burst of astrophysical activity will be both large enough and long enough -well over O(100 s) and spanning the time from before to after the gravitational wave burst -to significantly skew that mean. The rate of all tracks is an example of this class of selections, as shown in Fig. 1.
For all high-background selections with the exception of the few-MeV sub-supernova-like searches, the background events are nearly all uncorrelated, such that an excess can be quantified using Poisson statistics. For the few-MeV samples, many events are correlated -for instance, bursts of neutron captures from air showersso instead we measure the Gaussian width of the distribution of bin contents in a control region to determine the significance of excesses in a signal region. For the 1000 s windows, the control region is defined as beginning 500 s before the gravitational wave event timestamp and ending 5 s before. From there to 500 s after is the signal region. For the 45 s readouts, the control sample is defined as 10 to 40 s after the event, with the assumption that interesting activity is more likely in the first 10 s. Examples of each of these are shown in Fig. 2.
In the case of high-background samples with a restricted sky region, the background changes over time because the allowed region's zenith angle is changing and the cosmic ray flux is a function of zenith angle. Although the precise form, as a function of time, of the background rate is quite complex, for the relatively short time windows used in this analysis, we found it sufficient to fit a linear function to the observed rate (see Fig. 3).
For low-background samples without pointing dependence, the background level was determined by counting selected events in many uncorrelated time windows of the same trigger stream. For instance, the rate of FD events over 2 TeV was determined to be 1.0 × 10 −3 Hz. Likewise, this strategy is used for the supernova-like event search in the ND (see Fig. 4).
Finally, for low-background samples restricted by the LVC 90% allowed region, we use an ensemble of data with uncorrelated timestamps and select events using the same sky region, in zenith and azimuth, as the signal event (see, e.g., Fig. 5).
The significance of excesses was quantified by taking into account the trials factor given the number of bins searched for each gravitational wave event. Potentially interesting excesses in the first ten seconds after a gravitational wave event were considered special and only subject to a trials factor counting other bins within the first ten seconds. Each gravitational wave event was considered separately with no statistics computed using the ensemble of events.
We used a blind analysis, defining "significant" as being at least 3σ over background after the trials factor. No significant excesses were found. Post hoc visual inspection of the time distributions also did not reveal any features other than those expected in background.
Given NOvA's sensitivity to several-MeV neutrinos, but relatively small acceptance for higher energy events, the selections most likely to have shown a positive signal are the LVC-triggered supernova-like searches. Figure 6 shows the five gravitational wave events for which an LVC-triggered readout window is available in the FD. None reveals any evidence of a supernova-like burst.

A. Supernova-like neutrino fluence limits
Since no significant excesses were found in searches for a supernova-like signal, we set limits on the neutrino fluence. We assume the Garching models for 27 and 9.6 solar mass stars, without neutrino oscillations. Neutrino oscillations and other flavor-changing effects will modify the signal in NOvA and can either increase or decrease the observed interaction rate [29]. We assume that a potential supernova neutrino burst would occur in coincidence with the GW burst. The 27 solar mass model predicts a larger flux with a higher mean neutrino energy, and a time distribution more strongly peaked in the first second. The higher energy neutrinos are more efficiently detected, particularly by the Far Detector, leading to stronger limits as compared to the 9.6 solar mass model. The differing time distributions between the models, combined with background fluctuations, means the limits obtained for the  Raw events/s ND long readout: Supernova-like − S190701ah FIG. 4. Supernova-like event search for S190701ah, in the ND, using the LVC trigger. The background rate, shown as a solid line, is determined from a large number of uncorrelated time windows. The dashed line shows the expected signal for a 9.6 solar mass supernova at 10 kpc. The bin with four events has p = 9% taking into account the trials factor of the 45 bins in this plot alone and is not considered significant. two models are not simply proportional.
We perform a Bayesian analysis with a flat prior in neutrino fluence, profiling over the background level in the case of the FD (in the ND the background is fixed using uncorrelated time windows, see Section VI). Limits are shown in Table I. For the case that we read out both detectors in response to an LVC trigger, the mean 90% CL upper limit on neutrino fluence is 0.12 (0.24) × 10 12 cm −2 for the 27 (9.6) solar mass model. When only FD 10 Hz trigger data are available, the mean limits are 1.5 (4) × 10 12 cm −2 .
An upper limit on fluence, F 90 , can be converted to a lower limit on the distance, r 90 , to a hypothetical supernova: where N is the total number of neutrinos emitted. For the Garching 27 (9.6) solar mass model, N = 11 (6.8) × 10 57 . For all events with LVC-triggered readout in both detectors, distances closer than 22 kpc are excluded at 90% CL for the case of a 27 solar mass supernova. This limit varies by event and is as far as 29 kpc for S190602aq. These limits exclude 99% of the volume in which potential supernovae could occur in the Milky Way [32]. For the 9.6 solar mass model, distances up to 12 kpc are excluded in all cases for which we have LVC-triggered readout, with up to 15 kpc excluded for S190602aq. These exclusions cover 60-80% of potential galactic supernovae. Even for events in which only FD 10 Hz trigger data are available, 2.4-4 (4-8) kpc are excluded in the 9.6 (27) solar mass case, covering 4-12% (12-34%) of the galaxy.

VII. CONCLUSION
The NOvA detectors, which have sensitivity to signals, particularly neutrinos, in the MeV-TeV range, detected no significant excesses of events during the time around any of 28 gravitational wave events reported by the LIGO/Virgo Collaboration from September 2015 through July 2019 for which at least one NOvA detector was active. Sensitivity to MeV-scale events was best for S190521g, S190602aq, S190630ag and S190701ah, all binary black hole mergers during which we recorded 45 s of continuous data in both detectors.
The NOvA collaboration intends to continue operating both detectors and receiving LVC triggers through 2025. FIG. 6. Results of supernova-like neutrino search for events with FD LVC-triggered data. The readout of S190706ai is incomplete due to a late trigger, with progressively less available each second between 11 s and 14 s. A fit to a constant rate is shown as a solid red line in the upper and middle panes of each plot. Dashed lines show best fits for a 9.6 solar mass supernova at 10 kpc.