Search for sub-TeV Neutrino Emission from Novae with IceCube-DeepCore

The understanding of novae, the thermonuclear eruptions on the surfaces of white dwarf stars in binaries, has recently undergone a major paradigm shift. Though the bolometric luminosity of novae was long thought to arise directly from photons supplied by the thermonuclear runaway, recent gigaelectronvolt ( GeV ) gamma-ray observations have supported the notion that a signi ﬁ cant portion of the luminosity could come from radiative shocks. More recently, observations of novae have lent evidence that these shocks are acceleration sites for hadrons for at least some types of novae. In this scenario, a ﬂ ux of neutrinos may accompany the observed gamma rays. As the gamma rays from most novae have only been observed up to a few GeV, novae have previously not been considered as targets for neutrino telescopes, which are most sensitive at and above teraelectronvolt  ( TeV ) energies. Here, we present the ﬁ rst search for neutrinos from novae with energies between a few GeV and 10 TeV using IceCube-DeepCore, a densely instrumented region of the IceCube Neutrino Observatory with a reduced energy threshold. We search both for a correlation between gamma-ray and neutrino emission as well as between optical and neutrino emission from novae. We ﬁ nd no evidence for neutrino emission from the novae considered in this analysis and set upper limits for all gamma-ray detected novae.


INTRODUCTION
Novae are luminous outbursts that occur when a white dwarf in a binary system rapidly accretes matter from its companion star, leading to unstable nuclear burning on the surface of the white dwarf (Gallagher & Starrfield 1978).Surprisingly, it was recently discovered that these luminous events, typically identified in optical wavelengths, were often accompanied by gigaelectronvolt (GeV) gamma rays (Ackermann et al. 2014).Now, approximately one decade since the first discovery of gamma rays from novae with the Large Area Telescope (LAT) aboard NASA's Fermi satellite, there have been over a dozen novae identified in gamma rays, some identified in archival searches (Franckowiak et al. 2018), and more recently being announced in real time via chan-nels such as the Astronomer's Telegram (ATel)1 .For a recent review of novae, see Chomiuk et al. (2021).
The gamma-ray emission from novae was recently found to be strongly correlated in time with the optical emission, lending evidence for a common origin in shocks (Aydi et al. 2020).Novae have apparent efficiencies of nonthermal particle acceleration on the order of 0.3 − 1% (Li et al. 2017;Aydi et al. 2020).Although low compared to efficiencies inferred for adiabatic shocks in Galactic accelerators like supernova remnants (Morlino & Caprioli 2012), novae still provide us with a relatively nearby class of transients that can be used to understand more energetic and distant classes of nonrelativistic, shock-powered transients (Fang et al. 2020).The nonthermal gamma rays could, in principle, be produced from a variety of mechanisms, depending on the composition of the relativistic particles accelerated at the shocks.In the "leptonic" scenario, relativistic electrons Compton up-scatter the optical light or radiate via bremsstrahlung.In the "hadronic" scenario, accelerated relativistic ions interact with ambient gas, producing charged and neutral pions, which in turn decay into neutrinos and gamma rays, respectively.
For many recent observed novae, several studies favor the hadronic scenario.For example, large magnetic fields are required to accelerate particles to produce the GeV gamma rays detected from some novae, which is inconsistent with leptonic models (Li et al. 2017;Vurm & Metzger 2018).Moreover, recent observations of the spectrum of nova RS Ophiuchi during its 2021 outburst saw its spectrum extend beyond TeV energies, with fits to both H.E.S.S. and MAGIC data consistent with hadronic modeling (Acciari et al. 2022;Aharonian et al. 2022).
In addition to modeling based on gamma-ray observations, neutrinos could hold the key to distinguishing between the leptonic and hadronic models (Razzaque et al. 2010;Metzger et al. 2016).Furthermore, neutrinos could not only validate the hadronic scenario, but the level of any such neutrino flux would provide valuable information for understanding the environments of these shocks.
However, while there is a guaranteed flux of neutrinos that would accompany any gamma rays produced in the hadronic scenario, the energies of the gamma rays place a limit on the maximum energies of neutrino counterparts.On average, hadronic gamma rays carry approximately 10% of the initial relativistic ion energy, whereas the neutrinos carry about 5% of the initial energy.Thus far, with the exception of nova RS Ophiuchi, the gammaray spectra from novae detected by the Fermi -LAT have been modeled as a power law with an exponential cutoff, dN/dE ∝ E −Γ exp(−E/E cut ), where Γ is the photon spectral index and E cut is the cutoff energy (Ackermann et al. 2014;Franckowiak et al. 2018).Most individual novae have been fit with a cutoff of a few GeV, and a population analysis finds that the global spectrum from novae is well described with Γ = 1.71 ± 0.08 and E cut = 3.2 ± 0.6 GeV (Franckowiak et al. 2018).While some novae show no strong evidence for a cutoff with LAT data alone, observations with MAGIC of an individual nova, V339 Del, place a strong upper limit in the TeV band, necessitating a cutoff in the 100 GeV range, if not at lower energies (Ahnen et al. 2015).
It is for this reason that novae have not previously been targeted by neutrino telescopes.Optical Cherenkov neutrino telescopes, such as the IceCube Neutrino Observatory, typically have optimal sensitivities ≳ TeV, which is significantly higher than one could expect from even the most optimistic models for cosmic-ray acceleration in novae.Some previous searches for sub-TeV neutrino sources have been performed by Super-Kamiokande (Thrane et al. 2009;Abe et al. 2018) as well as by IceCube (Aartsen et al. 2016;Abbasi et al. 2022a), but these analyses have targeted other classes of astrophysical transients.
In this work, we describe the first search for neutrinos from novae using IceCube-DeepCore, a subarray of the IceCube Neutrino Observatory, which can be used to lower the energy threshold for astrophysical neutrino source searches.Although DeepCore has been used in past all-sky searches for generic transients (Aartsen et al. 2016;Abbasi et al. 2022a), this is the first time it has been used in a search for Galactic transients.
We begin in section 2 by describing the sample of novae used in this work, distinguishing those detected by the Fermi -LAT from those only detected in optical wavelengths.In section 3, we describe DeepCore in more depth, and list the details of the event selection used in our analyses.Then, in section 4, we describe the maximum likelihood approaches we implemented, including, in section 4.1, a search for correlations with individual gamma-ray light curves, and in section 4.2, a search for a cumulative signal from all optically detected novae.We summarize our results in section 5, and then discuss the implications of our results as well as of using DeepCore for neutrino astronomy in section 6.

NOVA SAMPLE
To create our list of novae, we began by using the catalog compiled for a search for gamma-ray emission from novae in Franckowiak et al. (2018).This catalog spanned the years 2008-2015, whereas the neutrino data used in our analysis does not begin until 26 April 2012 and extends until 29 May 2020.To gather information on more recent novae as well as to ensure that the catalog of all novae up to 2015 was exhaustive, we cross-referenced our sample against other compilations. 2,3For additional information on gamma-ray novae, we searched ATels, and include two more recent novae detected by the LAT at the >5σ level, V3890 Sgr and V1707 Sco.All gamma-ray detected novae used in this search are detected by the LAT at the 3σ level or above, with the exceptions of V745 Sco and V1535 Sco, which were only identified as "candidate" gamma-ray sources in Franckowiak et al. (2018), because they did not reach a detection at the 3σ level, but are included here for completeness.For the gamma ray detected novae, we use the same time window for the neutrino search as the LAT detection window, as specified in Table 1.Time windows for the majority of the novae match those used in Gordon et al. (2021).For the "candidate" novae V745 Sco and V1535 Sco, we use the times reported in Franckowiak et al. (2018).For the more recent novae V3890 Sgr and V1707 Sco, we use the detection times reported in their respective ATels4,5 .
Additionally, for all novae we require a measurement of the optical peak time.For novae through 2015, this was explicitly calculated in Franckowiak et al. (2018).For later novae, we accessed light curve information through AAVSO6 as well as through the Stony Brook/SMARTS Spectroscopic Atlas of Southern Novae7 , and, where available, we also included observations from the ASAS-SN transient server (Shappee et al. 2014).If a nova was the subject of a dedicated publication that included optical light curves, we took the information from the respective work.As some novae are believed to be discovered after peak brightness, and as some novae have light curves with only sparse sampling, we are not able to precisely determine the date of optical peak brightness; therefore, the novae V1404 Cen, V5854 Sgr, V1657 Sco, V3663 Oph, FM Cir, V3731 Oph, V3730 Oph, V2891 Cyg, and V1709 Sco are excluded from our analysis.
Our final sample contains 67 novae, 16 of which were detected by the LAT at GeV energies.The full table of all of the novae used in this work, including the dates of peak optical brightness as well as periods of gamma-ray detection, is included in Table 2 of Appendix A. Figure 1 displays the locations of these novae superimposed upon the sensitivity of the analysis, which will be described in Section 4.1.As expected, the novae are most commonly found near the Galactic plane, with a large number in the Galactic bulge region.It is worth noting that nova RS Ophiuchi is not included in this sample.As there is only low-energy neutrino data processed through 2020, there is no synchronous data available for an analysis of the emission seen from nova RS Ophiuchi in its 2021 outburst; however, this nova can be included in future analyses with more years of processed data.

ICECUBE-DEEPCORE
The IceCube Neutrino Observatory is a gigatonscale ice Cherenkov detector at the geographic South Pole (Aartsen et al. 2017).The detector consists of 86 vertical strings, each of which is instrumented with 60 digital optical modules (DOMs) that house 10-inch photomultiplier tubes.Neutrinos are detected indirectly via the Cherenkov radiation produced from relativistic charged particles created by neutrino interactions in the surrounding ice or nearby bedrock beneath IceCube.For the majority of the strings in the detector, the DOMs are evenly spaced between 1.45 km to 2.45 km below the surface of the ice, and the strings are arranged in a hexagonal grid.In the main array, the strings have a horizontal spacing of 125 m and the DOMs are spaced 17 m apart on each string, which optimizes the sensitivity of the detector at energies in the ∼TeV−PeV range, and has some sensitivity down to ∼100 GeV.In the center of the array, there is a more densely instrumented region of 8 specialized low-energy strings, with reduced spacings of both the strings (72 m apart) and the DOMs on each string (7 m apart), and the photomultiplier tubes in these DOMs have higher quantum efficiency than those in the main array.These 8 strings, as well as the 7 innermost strings of the main array, make up IceCube-DeepCore (Abbasi et al. 2012), and its reduced spacing and higher quantum efficiencies allow a reduced energy threshold down to ∼10 GeV.
In this analysis, we use a modified version of the GRECO (GeV Reconstructed Events with Containment for Oscillation) dataset that was developed for oscillation studies with DeepCore, as described in Aartsen et al. (2019), which we denote as the GRECO Astron-omy dataset for the remainder of this paper.This event selection is described in detail in Appendix B, including the selection criteria, the event selection response as a function of energy, as well as a description of the angular uncertainty of the events that are used in the sample.In short, the effective area is similar across the whole sky, and the selection is sensitive to neutrino emission in the energy range from a few GeV to about 10 TeV, depending on the assumed spectral hypothesis of a source.Depending on the energy and type of event morphology, events in the GRECO Astronomy dataset can have uncertainties in the angular resolution ranging from a few degrees to approximately 80 • , as shown in Figure 10 (Appendix B).The events used in this selection have angular uncertainties that are much larger than those present in ≳ TeV neutrino analyses, because the track lengths of outgoing muons (if present) are shorter in DeepCore, and shorter track lengths increase angular uncertainty.

ANALYSIS TECHNIQUES
The analyses performed for this work rely on an unbinned maximum likelihood approach that is common in searches for astrophysical neutrino sources (Braun et al. 2008).For this work, we implement an "extended likelihood" approach, which splits the entire dataset into a time period of interest (the "on-time window") with duration ∆T and the remainder of the dataset (the "offtime window"), with ∆T defined by the duration of the electromagnetic outburst for each nova.This approach has been used in many other searches for short-timescale neutrino transients (Aartsen et al. 2017;Aartsen et al. 2020;Aartsen et al. 2020).
For an on-time window with N total neutrino candidate events, we define the likelihood, L, which consists of the sum of a signal probability density function (PDF), S, and a background PDF, B, as where n s and n b are the signal and expected background event counts, respectively, and γ is the spectral index of the source, as defined in Equation 2 below.The index i iterates over all neutrino candidate events in the on-time window, and each of these events has a set of observables x i .The term outside of the product in the likelihood accounts for possible fluctuations in the rate of events, which is helpful in transient-style analyses where the expected number of background events is small (Barlow 1990).In order to calculate the expected number of events, we assume a neutrino number flux with a powerlaw spectrum: where γ is the spectral index, ϕ 0 the flux normalization, at a reference energy E 0 = 1 TeV.We maximize the likelihood with respect to both n s and γ, to characterize the normalization and spectrum of neutrino emission from a potential source.
Both the signal and background PDFs consist of terms describing the energy and the spatial distributions.Beginning with the signal PDF, S, the spatial component expects neutrinos from novae to be spatially associated with the nova location observed from photons.We define the PDF as it is used in Aartsen et al. (2017a), which is a function of the opening angle between the source location and the reconstructed event direction, ∆Ψ, and which can be approximated by the first-order non-elliptical component of the Kent distribution.The Kent distribution is an extension of a radially symmetric Gaussian that accounts for projection effects onto a sphere, which is necessary when handling events with large uncertainties in the reconstructed direction.The width of this distribution is governed by a per-event angular uncertainty, σ i .The energy term of the signal PDF is a function of each event's best-fit decl.and energy, and is calculated for a variety of possible source spectral indices.It is described fully in Aartsen et al. (2017b).
For the background PDF, B, both the energy and spatial components are parameterized directly from experimental data.The spatial component depends only on each reconstructed event's decl., and the energy component is a function of the event's reconstructed decl.and its reconstructed energy.The probability in R.A. is treated as a uniform distribution.The number of expected background events, n b , is calculated from offtime data to obtain an estimate on the all-sky rate of the sample.
The final test statistic (TS) is based on a log likelihood ratio between the best-fit signal hypothesis (with best-fit parameters ns , γ) and the background-only hypothesis (n s = 0, for which γ is undefined), and simplifies to

Gamma-ray correlation analysis
The first analysis of this work looks for neutrino emission coincident with gamma-ray emission.For each nova that has been detected in gamma rays (or was identified as a gamma-ray candidate in Franckowiak et al. ( 2018)), we use the time period of gamma-ray observations to determine the on-time window for the analysis.The time windows are consistent with what was reported for gamma-ray activity in Gordon et al. (2021) and with the time windows for the candidate emitters in Franckowiak et al. (2018).The time windows range from days to several weeks in duration.As the GRECO Astronomy dataset can only be used to search for transient neutrino sources due to high atmospheric backgrounds at these energies, we validated that we could search for neutrino emission on all of these timescales by injecting artificial signal events into our analysis and making sure that the analysis could recover an injected signal.One way of characterizing the analysis performance is by calculating the sensitivity, which is defined as the median one-sided Neyman upper limit (at the 90% confidence level) that one could expect to set under the background-only hypothesis.In Figure 1, we show the analysis sensitivity at all locations on the sky for an on-time window of 1 day, which is the minimum time window used in the gammaray correlation analysis.The sensitivity is better in the Northern Sky, where the atmospheric backgrounds are smaller, but the variation in analysis performance between the Northern and the Southern Sky is not as extreme as in higher-energy analyses, where the analysis performance can differ by orders of magnitude in the two celestial hemispheres.
Then, for each of the 16 gamma-ray novae with gamma-ray observations that were coincident with the GRECO Astronomy dataset, we perform the likelihood analysis as described above.In addition to calculating the best-fit parameters and test statistics, we also calculate p-values for quantifying the compatibility with the background-only hypothesis.These p-values are onesided and are calculated by comparing the observed test statistic to an expected distribution of test statistics created by performing pseudo-experiments assuming the background-only hypothesis.

Stacking analysis
The second analysis of this work stacks together the contribution from multiple novae in order to improve the sensitivity of the search.For this construction, the PDFs in the likelihood are weighted sums of the PDFs from each individual source, as in Aartsen et al. (2017a).While stacking analyses have the advantage of improving the average sensitivity of the sample, they require testing a particular hypothesis about the distribution of fluxes from the individual sources.There is no conclusive or agreed upon model for calibrating neutrino fluxes based on observations of novae across the electromagnetic spectrum.However, as the production of neutrinos is intimately linked to the production of gamma rays, the gamma rays could provide a handle on the expectation for a corresponding neutrino flux.If, on the other hand, the shocks that are producing neutrinos are optically thick for gamma-ray emission, then any potential neutrino flux could overshoot the corresponding gamma-ray flux (Bednarek & Śmiałkowski 2022).In this case, the optical emission could provide the best handle on the level of a neutrino flux, as it has been suggested that the total energy that goes into particle acceleration could be correlated with the optical energy output (Fang et al. 2020).
We therefore test two different weighting schemes: (1) gamma-ray stacking and (2) optical stacking.For gamma-ray stacking, we use the 16 novae that were detected in gamma rays, and weight them based on the normalization of the gamma-ray flux at 100 MeV that was fit by Fermi -LAT.For the optical stacking weighting scheme, we no longer require that a nova be detected in gamma rays.Figure 2 shows the distribution of peak optical magnitudes for all novae in the sample.We weight the novae based on their peak optical fluxes, which are converted from the optical magnitudes.We break this distribution up into novae that were detected in gamma rays and those that were not.In general, the novae that were detected in gamma rays were brighter optically, which might be an indication that all novae are gamma-ray emitters, but some are at a level that is too dim to be detected by Fermi -LAT.
We also require an on-time window for each of the novae in both of these weighting schemes.As mentioned above, using the GRECO Astronomy dataset to search for point sources is only feasible when searching for transient emission.When stacking together sources, this becomes especially important because, if there are k sources each with an on-time window of ∆T , the effective background is equivalent to a single source with a time window of k∆T .We find that when stacking sources, this limits us to searching for emission within a smaller time window than we could use for the gammaray correlation analysis.After injecting artificial signal into the analysis over a variety of time windows, we find that the analysis is able to precisely recover the amount of injected signal for on-time windows of 1 day or less.To maximize the amount of signal that we can integrate over, both weighting schemes used ∆T = 1 day for the duration of the stacking analyses.For many of the gamma-ray novae, the gamma-ray light curves are not resolved enough to isolate the peak day of emission.Additionally, light curves of recently detected bright gamma-ray novae, such as V906 Car and ASASSN-16ma, show correlations between the optical and gamma-ray light curves, with delays between the two light curves of significantly less than 1 day.For these reasons, we center the time window on the date of peak optical emission for both the gamma-ray stacking and optical stacking analyses.These dates, as well as the peak optical magnitudes, are included in Table 2 in Appendix A.
For the gamma-ray correlation analysis, we find that to be sensitive to an individual source with a spectrum F ν+ν (E|γ = 2.0) in the Northern Sky with an on-time window of 1 day, we require 5 neutrinos.On the other hand, an observation of a source with the same spectrum in the Southern Sky with a long on-time window would require around 50 neutrinos.When we stack sources together for the gamma-ray weighting scheme, we find that we would need around 16 (53) signal events for an F ν+ν (E|γ = 2.0) (F ν+ν (E|γ = 3.0)) spectrum, whereas for the optical stacking weighting scheme we would need around 27 (130) signal events for an F ν+ν (E|γ = 2.0) (F ν+ν (E|γ = 3.0)) spectrum.In terms of per-source contributions, this means, on average, reducing the requirement from multiple detected neutrinos from each nova to 1 (0.4) detected neutrinos from each nova for the gamma-ray (optical) weighting scheme.

Systematic Uncertainties
We estimate the effects of systematic uncertainties on our results by varying the most important systematic effects in IceCube.These systematics include ice properties, such as scattering and absorption of photons in the ice; relative DOM efficiency; and properties of the refrozen column of ice around the DOMs resulting from the drilling process, known as "hole ice."In order to calculate the effect of systematic uncertainties on our upper limits, we produced two additional Monte Carlo datasets for each systematic effect, with a discrete value of the systematic chosen to bracket the uncertainty of that systematic effect.We simulate systematics datasets with ±10% DOM efficiency, ±10% absorption coefficient, ±10% scattering coefficient, and ±1 σ variations in hole ice optical properties, which were chosen to match those used in Abbasi et al. (2022a).
We recalculate our sensitivities for each nova at 3 different spectral indices (γ =2.0, 2.5, and 3.0) using each of these datasets.We find the ratio of the sensitivity with the systematic included to the sensitivity without the systematic included, and average this value over all novae for each systematic.We then sum the components of these uncertainties in quadrature to find a total systematic uncertainty for each of the three spectral indices.We find the total systematic uncertainty on the flux upper limit, averaged over novae, to be within ±21% for a spectral index of 2.0, ±20% for an index of 2.5, and ±13% for an index of 3.0.We rescale our power-law upper limits presented in Figure 4 using the upper bound of these calculated systematic uncertainties for each spectral index, to provide a conservative upper limit.

RESULTS
After calculating test statistics for each analysis on true experimental data, no significant signal is detected, either in the gamma-ray correlation analysis or in either of the stacked analyses.
In Table 1, we include all of the results from the gamma-ray correlation analysis.We also show the observed test statistics for each of the gamma-ray detected novae, as well as the background test statistic distribution calculated by performing pseudo-experiments assuming a background-only hypothesis in Figure 3.In addition to the time windows of the analysis, we include the best-fit parameters and pre-trial p-values, before accounting for the trials correction factor accrued from analyzing multiple sources.The most significant followup comes from our investigation of the nova V1535 Sco, which had a pre-trials p-value of 2.4×10 −3 , which gives a post-trials p-value of 3.8% after accounting for the number of novae investigated, which we find to be consistent with expectations from atmospheric backgrounds.Nova V1535 Sco was not significantly detected in gamma rays, and was only identified as a candidate gamma-ray emitter in Franckowiak et al. (2018).

Power-law Upper Limits
As no significant signal was detected, we calculate upper limits on neutrino fluxes from novae during time periods of coincident gamma-ray detections.These limits are displayed in Figure 4. We also include a comparison to the gamma-ray fluxes as measured by Fermi -LAT.When comparing the levels of the observed gamma-ray fluxes and the upper limits set by this analysis, it becomes apparent that more sensitive detectors might be required before one can expect to detect a neutrino signal from novae.
Similar to the gamma-ray correlation analysis, no significant signal was detected in the stacking analysis, for either the gamma-ray weighting scheme or the optical stacking scheme.The likelihood maximization for each of the weighting schemes finds a best-fit of ns = 0, which yields T S = 0, and a pre-trial p-value of 1.0 for each of these weighting schemes.In Figure 5, we show the likelihood spaces for our observed data, and we compare this against simulations in which we injected artificial signal, to contrast how our observed data appears relative to what we might expect for a bright signal from novae.In the simulation panels, we show how well the properties of the injected signal can be recovered by the analysis.There is a tendency for the analysis to underestimate the signal, which is especially visible in the optical stacking (Figure 5, lower-right panel).This predominantly comes from the large angular uncertainties of events in the GRECO Astronomy dataset (see Appendix B, Figure 10), which causes low-energy signal events to appear background-like, and the assumption that the uncertainty contours are symmetric and Gaussian.This bias does not affect any constraints we set, as we calculate our upper limits using the one-sided Neyman construction, which handles this bias self-consistently.Note.MJDstart and MJDstop represent the beginning and end of when a nova was detected by Fermi-LAT, respectively.This is the same time window used for the neutrino search.The observed test statistic (T S), best-fit number of signal events (ns) and best-fit spectral index (γ) for the neutrino search from each nova are given here.The p-values shown are before accounting for the factor accrued from performing multiple searches.The final column gives the time-integrated flux upper limit for each nova, assuming a power-law spectrum with γ = 2.0 as defined in Equation 2 and including the systematics discussed in Section 4.3.
As we did not detect any significant signal in the stacking analysis, we calculate upper limits for a variety of power-law spectra, for both of the weighting schemes.These limits are shown in Figure 6, both in terms of the stacked energy-scaled time-integrated fluxes as well as in terms of the average number of signal events to which each of these fluxes correspond.These limits are at the level of the sensitivities we described in Section 4, which is about 1 (0.4) events from each nova, on average for the gamma-ray (optical) weighting scheme for an assumed F ν+ν (E|γ = 2.0) signal spectrum.

Physical model upper limits
In addition to calculating limits for generic unbroken power laws for the gamma-ray correlation analysis, we also inject the spectral shapes from the models developed in Bednarek & Śmiałkowski (2022), which suggests that there could be neutrino fluxes at a level larger than would be naively assumed from the levels of the gammaray fluxes.Even though the likelihood assumes a power law in the signal PDF, we can calculate an upper limit on these model fluxes by keeping the likelihood fixed but injecting the model fluxes with various normalizations until the resulting expected test-statistic distributions are distinguishable from our observed result at 90% confidence.We show these limits on the model fluxes in addition to the power-law upper limits in Figure 4.
For the cases of the fluxes from Bednarek & Śmiałkowski (2022), the spectral shapes are a function of the maximum proton energy as well as the magnetic field strength on the surface of the white dwarf at the equator.We chose the physical parameters corresponding to the tested model in Bednarek & Śmiałkowski (2022), which yields the highest number of events in the GRECO sample.We show two different maximum proton energies, E max , for the limits shown in Figure 4 and used a nominal B-field value of 10 8 G, although the interpretation does not change much when varying the parameters of the model and comparing to our analysis results.In total, depending on the exact physical parameters of the nova, the models of Bednarek & Śmiałkowski (2022) predict that the GRECO Astronomy data sample could contain O(0.01) neutrinos from the weakest novae and O(1) neutrinos from the brightest novae, in the most optimistic scenarios.Our limits, on the other hand, correspond to higher fluxes that correspond to a larger number of signal events.These limits are functions of the source position, time window, as well as the observed test statistic for each nova.For novae with shorter gamma-ray detections or with smaller observed Observed Extrapolated Figure 4. Upper limits on neutrino emission from all of the novae analyzed in the gamma-ray correlation analysis.Upper limits on power laws span the 90% energy ranges discussed in Figure 9, and include systematic effects as discussed in Section 4.3.In addition to constraining power laws, we also inject the spectral shapes from models in Bednarek & Śmiałkowski (2022), and our upper limits on those spectra are shown in pink.We compare these fluxes to the measured gamma-ray fluxes (gray).For those gamma-ray detected novae that did not show evidence for a cutoff in the gamma-ray spectra, we show the lines as dotted past the global cutoff energy found in Franckowiak et al. (2018).
test statistics, such as V1707 Sco, V392 Per, V339 Del, V745 Sco, or V407 Lup, we find that these limits correspond to an expectation of 10-30 signal events, whereas for the other novae in the sample, the fluxes must be large enough to give an expectation of 50-200 signal events.As these upper limits are above the modeled fluxes from each nova, we cannot constrain these models.

DISCUSSION AND CONCLUSION
We report on the first search for neutrino emission from novae in the ∼ GeV to ∼ 10 TeV energy range.Although there is promising evidence that novae are hadronic particle accelerators, the level of any corresponding neutrino flux must be at a level less than that which we are sensitive to with IceCube-DeepCore and current analysis techniques.Current models predict that, with current neutrino detectors, the novae that have been detected with photons in the last decade could be emitting neutrinos at levels that would result in O(0.01 − 1) detectable neutrinos in IceCube-DeepCore per nova.This is below our sensitivity, even when stacking together contributions from multiple novae.
In Section 2, we discussed that we could only incorporate novae from 2012-2020 into our sample, because newer low-energy data from IceCube-DeepCore were not yet available.As a result of this, we were not able to include nova RS Ophiuchi, the first nova detected by ground-based gamma-ray observatories, into this analysis.It is worth noting that because there was very highenergy gamma-ray emission detected from this object, a separate analysis (Pizzuto et al. 2021) was performed to search for higher-energy neutrino emission from nova RS Ophiuchi by analyzing data available from IceCube with low latency (Abbasi et al. 2021).That search found no evidence for neutrino emission from nova RS Ophiuchi.However, that search used a data sample in which 90% of detected events from a source at the location of nova RS Ophiuchi with an unbroken F ν+ν ∝ E −2 power-law spectrum would have true neutrino energies between 2 TeV and 10 PeV.This exceeds the energy range expected for neutrino emission given the detected electromagnetic emission (Acciari et al. 2022, their Figure 12).Using the neutrino flux predictions from Acciari et al. (2022), the analysis presented in this work would still not be sensitive to such a flux.Future analyses that could incorporate nova RS Ophiuchi could constrain sce- Here, ns refers to the total signal events for the entire stacked sample.The top row shows our observed data, and the bottom row shows a sample pseudo-experiment with injected signal.Contours denote 50%, 90%, and 99% containment assuming Wilk's theorem with two degrees of freedom.For our observed data (top), γ is undefined when ns = 0, so we do not include a best-fit point on these panels.Optical stacking (UL, 90%) Figure 6.Upper limits on the stacked neutrino fluxes for the gamma-ray stacking (blue) and optical stacking (orange) weighting schemes.The left plot shows the energy-scaled, time-integrated, fluxes, whereas the right plot shows the corresponding expected number of signal neutrino events.Both plots show the total contributions from all sources for each weighting scheme.
narios in which the neutrino flux exceeds the observed gamma-ray flux.Future improvements to O(GeV) neutrino detectors could drastically change the prospects for detecting neutrinos from novae.One of the chief challenges with neutrino detection from novae is that many of the predicted spectra are expected to peak at the lowest energies to which neutrino detectors are sensitive.However, planned detectors such as the IceCube-Upgrade (Ishihara 2021) or the ORCA configuration of KM3NeT (Assal et al. 2021) will represent a significant improvement in effective area between 1-10 GeV with respect to the GRECO Astronomy sample used here.The IceCube-Upgrade (Ishihara 2021) will also provide improved angular resolution for >10 GeV events.Additionally, searching for neutrino emission on longer timescales, especially when stacking, becomes difficult with lowenergy events.There are already improvements being made to low-energy reconstructions with IceCube-DeepCore (Abbasi et al. 2022b), and reconstructions with next-generation observatories have the promise of being even better due to an enhanced understanding of systematic uncertainties.Combining these improvements -enhanced effective area at the lowest energies to boost signal and better background discrimination from better event reconstructions which allow for longer integration times -have the promise of greatly enhancing sensitivity to low-energy transients such as novae.
The ability to look at longer timescales is especially important in light of the recent H.E.S.S. observations of nova RS Ophiuchi (Aharonian et al. 2022).These observations suggested that particles that are accelerated to the highest energies only attain these energies a couple of days after the peak observed in GeV gamma rays or in optical bands.While the individual catalog analysis searched for emission on timescales of multiple days for most novae, the stacked analysis was limited to a time window 1 day in duration because of the overwhelming atmospheric backgrounds.Future analyses that feature better spatially reconstructed events could reduce the effective background and have a boosted sensitivity when searching for neutrino emission correlated with this TeV peak in the gamma-ray spectra.
Even with these improvements, it may be the case that we will need a particularly nearby or bright nova to conclusively discriminate between leptonic or hadronic models of particle acceleration in novae, or to potentially discriminate between neutrino emission models that suggest neutrino emission internal to the shocks (Bednarek & Śmiałkowski 2022) or that would predict similar levels of neutrino and gamma-ray emission (Fang et al. 2020).There are several nearby recurrent novae, eg.T Coro-nae Borealis, which has an estimated distance that is approximately 3 times nearer than the best-fit distance of ∼ 2.4 kpc for nova RS Ophiuchi.If a future nova outburst from this system were to occur and have an emitted energy similar to RS Ophiuichi, the potential neutrino flux from this object could be detectable, especially with next generation detectors.For this work, we use a modified version of the GRECO (GeV Reconstructed Events with Containment for Oscillation) event selection, which was initially developed for the tau neutrino appearance analysis described in Aartsen et al. (2019).The original GRECO dataset, which was referred as Analysis A in that reference, consisted of 3 years of data from 2012 April through 2015 May covering the full sky and all neutrino flavors.This original dataset was used in an initial untriggered all-sky search for sub-TeV transient neutrino emission (Abbasi et al. 2022a).
An updated version of the GRECO dataset is used here, incorporating updates and reoptimizations developed since the original dataset was produced to identify neutrino candidate events below about 1 TeV starting in the deep, clear ice of the DeepCore detector (Aartsen et al. 2019).Simulations of the updated dataset -here referred to as the GRECO Astronomy dataset -use an updated version of the GENIE neutrino generator (Andreopoulos et al. 2010) and include new charge calibrations for both data and simulation, an updated model of the glacial ice, and newly available simulations of atmospheric muons using an implementation of MuPAGE (Carminati et al. 2008) referred to as "MuonGun."Thirteen variables based on the original GRECO dataset are combined into a newly trained boosted decision tree (BDT).These include charge and time-evolution variables, position of the DOM with the earliest local coincidence hit in the detector, the amount of causally connected charge deposition found in the larger IceCube detector, and preliminary directional reconstruction.The results of the BDT training, shown in Figure 7, show that the data is well modeled with simulation weighted by atmospheric fluxes.The GRECO Astronomy dataset provides an improved effective area relative to the dataset used in the tau neutrino appearance search, with the largest improvements above 100 GeV.While the background rates are also higher, searches for short transient phenomena like those performed in this work are not significantly affected by the additional background contamination.The effective area of this sample is shown in Figure 8, and we compare it against the IceCube "Gamma-ray Follow-up" (GFU) sample, an event selection that is often used when searching for high-energy astrophysical transients (Aartsen et al. 2017c).To highlight the complementary nature of the different event selections, GRECO and GFU, we show the energies one could expect to detect from a source with various spectral hypotheses in Figure 9.Although using lower-energy events to search for astrophysical neutrinos is not without limitations, as described below, it opens up a new energy range to search for multimessenger sources.
The final sample has an average all-sky rate of 4.6 mHz, consisting of about 60% neutrino events from cosmic-ray interactions in the atmosphere and 40% muons from atmospheric air showers.The sample contains large backgrounds relative to the total expected astrophysical signal, but these backgrounds are strongly suppressed in searches for short astrophysical transients.
Events observed in the detector can be classified into two distinct morphologies: tracks and cascades.Tracks are elongated signatures and are the result of muons traversing the detector, either from cosmic-ray interactions or from charged-current (CC) muon-neutrino interactions.Cascades, on the other hand, are more spherical, and are the result of neutral current interactions from all neutrino flavors or electron or tau neutrino CC interactions.The long lever arm of tracks allows for much better pointing resolution than for cascades.
Neutrino events starting in the detector provide a third morphology due to the physics of the interaction.Neutrinos interacting inside of the detector first produce a hadronic cascade due to deep inelastic scattering interactions off of the nuclei in the ice.Following this initial hadronic cascade, neutrinos may produce an outgoing muon track (ν µ chargedcurrent), an electromagnetic shower (ν e charged-current), a tau lepton and associated decay (ν τ charged-current), or a hidden outgoing neutrino (neutral current).The relative brightness of the initial hadronic interaction and the 10 1 10 2 10 3 10 4 10 5 E (GeV) A eff , [m 2 ] -90 < < -5 , GRECO -5 < < 90 , GRECO -90 < < -5 , GFU -5 < < 90 , GFU Figure 8. All-flavor, ν + ν averaged effective area for the low-energy GRECO event selection used for this analysis (solid lines).We include, for comparison, a higher-energy νµ + νµ-only event selection that is often used to search for astrophysical neutrino transients (dashed lines).

Median
Central 68% GRECO (this analysis) GFU (high-energy tracks) Figure 9. Central 68% of the expected true neutrino energies for sources with Fν+ν ∝ E −2 (left), E −2.5 (middle), or E −3 (right) spectra.We compare this analysis (pink) to a higher-energy event selection (gray).This distribution only shows the distribution of true energies that are present in the dataset, although higher-energy and better reconstructed events tend to have a larger effect on the analysis sensitivity than lower-energy events.subsequent emission profile of the outgoing lepton degrades and, therefore, limits our ability to clearly identify the type of neutrino involved as the energy decreases.
The GRECO Astronomy dataset contains events of relatively low energy, with low light emission, and therefore few photons reach the optical modules due to the sparseness of the array.A typical event in the GRECO Astronomy sample often has around 20-30 hits or less over a relatively small spatial extent.Any muon track contained in the GRECO Astronomy dataset also tends to be short, crossing one to two strings on average, limiting our ability to reconstruct events or differentiate low-energy tracks from cascades.
These limitations result in large directional uncertainties of the events.For tracks at low energies, the pointing resolution is worse than what is achieved in searches for high-energy neutrino sources.For comparison, in analyses that use higher-energy track events (E ν ≳ 1 TeV), the median opening angle between the true neutrino direction and the reconstructed event direction is, on average, < 1 • .For lower-energy events, such as those used here, each neutrino interaction deposits much less energy in the detector, which leads to events that are much more poorly localized, with average opening angles between the true neutrino direction and the reconstructed event direction on the order of tens of degrees.Additionally, at lower energies, the kinematic opening angle becomes nonnegligible in comparison to the resolution of the detector.
In Figure 10, we show the average angular separation between the true neutrino direction and the reconstructed event direction for all flavors of neutrino, as both track and cascade events are included in the GRECO sample.The reconstruction algorithm was designed for the original neutrino oscillation event selection and is described at length in Aartsen et al. (2019).The reconstruction takes a maximum likelihood approach where the underlying hypothesis assumes that every event begins with an electromagnetic or hadronic shower and then has an outgoing, finite muon originating at the same interaction vertex.Some of the other assumptions that go into this reconstruction, such as assuming that the muons are in the minimum ionizing regime, reflect the fact that the reconstruction was designed specifically for lower-energy events (E ν ≲ 300 GeV), and is outperformed by other reconstruction algorithms at higher energies.The mixed track and cascade hypothesis makes the likelihood an eight-dimensional parameter space, and performing the reconstruction can be extremely computationally intensive, with each event requiring about 2 minutes on average to identify a minimum.Minimizations are performed over the full space using "Nestle," a tool for posterior probabilities using a nested sampling (Barbary 2015;Feroz et al. 2009;Mukherjee et al. 2006;Shaw et al. 2007;Skilling 2006Skilling , 2004)).
Instead of using approaches that require fine scans of the likelihood space as a function of the direction (Neunhoffer 2006), a more computationally efficient approach that only relies on the final event observables was used to estimate the angular uncertainty of each event.To accomplish this, a random forest was trained using a three-fold cross validation technique with neutrino simulation that contained both track-like and cascade-like events.This model outputs a singular value, σ, that describes the angular uncertainty of each event, which is used in the likelihood-based approaches implemented to search for pointlike astrophysical neutrino emission.It was found that this approach of

Figure 1 .
Figure 1.Locations of novae investigated in this work.Those novae that were identified in GeV gamma rays are shown in orange, while those only detected at other wavelengths are shown in gray.These novae are shown on top of the sensitivity (described fully in Section 4) our unstacked analysis (described in Section 4.1) would have to a single nova, assuming an analysis time window of 1 day.

Figure 3 .
Figure 3.Comparison of observed test-statistics (gray dashed lines) to expectations from 10,000 pseudo-experiments performed under the assumption of the null hypothesis (blue).Each panel represents the gamma-ray correlation analysis for an individual nova.

Figure 5 .
Figure 5. Stacking likelihood spaces for the gamma-ray stacking analysis (left) and for the optical stacking analysis (right).Here, ns refers to the total signal events for the entire stacked sample.The top row shows our observed data, and the bottom row shows a sample pseudo-experiment with injected signal.Contours denote 50%, 90%, and 99% containment assuming Wilk's theorem with two degrees of freedom.For our observed data (top), γ is undefined when ns = 0, so we do not include a best-fit point on these panels.

Figure 7 .
Figure7.The expected contributions from each atmospheric background component, including neutrinos and simulations of atmospheric muons.Neutrino charged-current (CC) interactions for each neutrino flavor, neutral-current (NC) interactions for neutrinos of all flavors, and atmospheric muons are labeled in the legend, with the total of all contributions shown in gray.Error bars shown are statistical and do not include any relevant systematic uncertainties.Events with a BDT score below 0.13 are removed, yielding a sample with 60% atmospheric neutrinos and 40% atmospheric muons.

Figure 10 .
Figure10.Distribution of the angular separation between reconstructed event direction and true neutrino direction as a function of energy, at 3 confidence levels (20%, 50%, 80%).All neutrino flavors, weighted according to the per-flavor effective area, are included.
Franckowiak et al. (2018)f the novae used in this analysis.Those novae that were detected in gamma rays or were identified as candidate gamma-ray emitters inFranckowiak et al. (2018)are shown in the orange histogram, whereas the blue histogram shows novae that were not detected in gamma rays.The gray histogram shows the total of all novae used in the analysis.

Table 1 .
Results from the Gamma-ray Correlation Analysis.