Search for High-Energy Neutrino Emission from Galactic X-ray Binaries with IceCube

We present the first comprehensive search for high-energy neutrino emission from high- and low-mass X-ray binaries conducted by IceCube. Galactic X-ray binaries are long-standing candidates for the source of Galactic hadronic cosmic rays and neutrinos. The compact object in these systems can be the site of cosmic-ray acceleration, and neutrinos can be produced by interactions of cosmic rays with radiation or gas, in the jet of a microquasar, in the stellar wind, or in the atmosphere of the companion star. We study X-ray binaries using 7.5 years of IceCube data with three separate analyses. In the first, we search for periodic neutrino emission from 55 binaries in the Northern Sky with known orbital periods. In the second, the X-ray light curves of 102 binaries across the entire sky are used as templates to search for time-dependent neutrino emission. Finally, we search for time-integrated emission of neutrinos for a list of 4 notable binaries identified as microquasars. In the absence of a significant excess, we place upper limits on the neutrino flux for each hypothesis and compare our results with theoretical predictions for several binaries. In addition, we evaluate the sensitivity of the next generation neutrino telescope at the South Pole, IceCube-Gen2, and demonstrate its power to identify potential neutrino emission from these binary sources in the Galaxy.


INTRODUCTION
Cosmic rays (CRs) up to several PeV, the "knee" in the CR spectrum, are believed to be of Galactic origin. However, where and how these CRs are accelerated remains an open question. Interactions of hadronic CRs in the Galaxy will lead to the production of pions, which subsequently decay into gamma rays and neutrinos with energies potentially reaching hundreds of TeV. Unlike gamma rays which could also be produced by accelerated electrons/positrons, high-energy neutrinos would be the smoking-gun for CR interactions as they are the only way to produce neutrinos. In 2013, IceCube discovered TeV-PeV neutrinos of astrophysical origin (Aartsen et al. 2013a(Aartsen et al. , 2014. Nevertheless, the sources of those high-energy neutrinos are yet to be identified. Since the discovery of high-energy cosmic neutrinos by IceCube, many studies have been conducted searching for neutrino emission from point-like sources, extended regions, and the diffuse emission from CRs interacting with the Galactic interstellar medium. The isotropic distribution of neutrino events in IceCube suggests a dominant contribution from extragalactic sources and constrains the Galactic contribution to the diffuse neutrino flux to less than 14% above 1 TeV (Aartsen et al. 2017a).
X-ray binaries (XRBs) are binary systems that emit X-rays and consist of a compact object (neutron star (NS) or black hole (BH)) and a non-compact companion star. The mass from the companion is accreted onto the compact object due to the strong gravitational attraction forming an accretion disk. Depending on the mass of the companion star, XRBs can be divided into highmass XRBs (HMXBs) and low-mass XRBs (LMXBs). The companion star of an HMXB is massive, usually an O or B-type star, and the accretion is mainly a result of stellar wind capture (Liu et al. 2006). LMXBs have companion stars of A-type or later, and the mass transfer is mostly caused by Roche lobe overflows (Liu et al. 2007). These systems are bright in X-rays and sometimes in gamma rays. XRBs have been proposed as sites of CR acceleration and hadronic interactions since the 1980s (e.g., Gaisser & Stanev 1985;Berezinsky et al. 1985). Microquasars, which are XRBs with relativistic jets, are miniature analogs of quasars, a category of extremely luminous active galactic nuclei. It is expected that the similarity may not only be morphological, but also in the physics of the formation of the accretion disk and the behavior of the jet, at different luminosity scales. Hadronic processes in jets of microquasars have been widely discussed. Protons can be accelerated in the jet, and pions are generated through interactions with the external radiation field of the accretion disk and/or internal synchrotron photons (pγ) (Levinson & Waxman 2001;Distefano et al. 2002;Romero & Vila 2008). Other scenarios focus on hadronuclear interactions (pp), e.g., jet-cloud/wind interactions when the jet traversing the matter field of the ejected clouds or stellar wind from the companion star (Aharonian & Atoyan 1991; Romero et al. 2003;Bednarek 2005). For other XRBs without collimated beams, hadronic interactions happen in a wider shocked region. CR can be accelerated in the magnetosphere of a spinning NS and then interact with matter from either the accretion disk or the companion star (Gaisser & Stanev 1985;Kolb et al. 1985;Berezinsky et al. 1985;Cheng & Ruderman 1989;Anchordoqui et al. 2003).
Some XRBs have been observed in TeV gamma rays, which may suggest acceleration of hadrons besides an origin of accelerated leptons. Neutrino production in those sources and the potential of IceCube detection have been discussed, e.g., LS 5039 (Aharonian et al. 2006), LS I +61 303 (Christiansen et al. 2006;Torres & Halzen 2007), and SS 433 (Reynoso & Carulli 2019;Kimura et al. 2020).
XRBs exhibit both outbursts and periodic emissions. It is reasonable to hypothesize that the possible neutrino emission is related to either the periodicity or the X-ray outbursts. Conducting time-dependent studies incorporating these hypotheses provides the advantage of suppressing background. Both time-integrated and time-dependent analyses searching for neutrino emission from the entire sky have been performed by Ice-Cube (Abbasi et al. 2012;Aartsen et al. 2013bAartsen et al. , 2014Aartsen et al. , 2015 and ANTARES (Adrian-Martinez et al. 2014;Albert et al. 2017) previously, without significant signal detection. Here, we present a comprehensive high-energy neutrino search focusing on XRBs using 7.5 years of IceCube data. This study explores the periodic, X-ray outburst-correlated, and persistent emission scenarios of the possible neutrino flux from XRBs, while covering a broader list of sources compared to previous studies.

DETECTOR & DATA SET
The IceCube Neutrino Observatory, built under the ice surface at the South Pole, is a 1 km 3 Cherenkov neutrino telescope. The detector is composed of an array of digital optical modules (DOMs), each equipped with a photomultiplier tube (PMT) and on-board read-out electronics. The PMTs in the DOMs collect the Cherenkov photons emitted by the relativistic charged particles produced in neutrino interactions when traversing the ice. The number of detected photons and their arrival time information are used to reconstruct the energy and direction of each event. Different neutrino flavors and interactions can result in different event signatures inside IceCube detection volume. Among them, muon tracks from ν µ /ν µ charged-current interactions can be reconstructed with a good angular resolution, typically 1 • (Aartsen et al. 2017b), which makes them ideal for neutrino point source searches.
In this search, we use 7.5 years of all-sky muon track data collected between 2011-05-13 and 2018-10-14, with an effective livetime of 2711 days. The data sample being used has an event selection focusing on high-quality through-going neutrino track events from the entire sky with an energy threshold near 100 GeV. The final event rate is ∼6 mHz, yielding a total of 1,502,612 events. A better sensitivity in the Northern Hemisphere is expected due to the suppression of the atmospheric muon background by Earth for up-going events. Details of the data sample are described in Aartsen et al. (2017c).

ANALYSIS
The three analyses use an unbinned maximum likelihood-ratio method (Braun et al. 2008(Braun et al. , 2010 to search for an excess of neutrino events above the background of atmospheric muon, atmospheric and isotropic astrophysical neutrino flux in the direction of selected XRBs. In all analyses, the likelihood function describing the signal includes both a directional correlation and an energy spectrum assumed to follow a power law, i.e., ∝ E −γ . For the time-dependent analyses, each incorporates a unique temporal term in the likelihood to model the signal. As the data is expected to be backgrounddominated, which is uniform in time, the background probability density function (PDF) is constructed from time-randomized data. The test statistic (TS) is obtained by maximizing the likelihood ratio with respect to a set of parameters, which include the number of signal events (n s ) and a power-law spectral index (γ). Other parameters are analysis-specific, e.g., time of the neutrino emission and the duration. A more detailed description of the method can be found in the appendix Sec. A. The sources studied are from a catalog of 114 Note-The most significant source in each analysis with its best-fitted TS,ns,γ. Both post-trial and pre-trial (bracketed) p-values are shown. The trial correction is calculated by considering the number of sources studied in each analysis respectively and has a relation ppost = 1 − (1 − ppre) N , where N is the source number.
Galactic HMXBs (Liu et al. 2006) and 187 LMXBs (Liu et al. 2007). Furthermore, we added 7 TeV XRBs from TeVCat 1 very-high-energy gamma-ray catalog (Wakely & Horan 2008), which are not in the HMXB or LMXB catalog. Additional selection criteria are then applied to this initial source list in the two time-dependent analyses (as explained in Sec. 3.1 and Sec. 3.2), leading to two different source lists with an overlap. We also include a time-integrated search on four individual sources. The rest of this section includes the description of each analysis and its results. Table 1 summarizes the most statistically significant source in each analysis.

Periodic Analysis
Periodic searches have been previously performed by IceCube using data of 2007-2009 with a partially-built configuration and 2008-2012, including the first year of the completed detector (Abbasi et al. 2012;Aartsen et al. 2015), assuming that the neutrino emission from XRBs is modulated by orbital periods. In the periodfolded phase space, these neutrinos should appear as at the same phase. As an improvement on prior work, the temporal signal PDF in the phase space is modeled by a von Mises distribution instead of a Gaussian distribution (see Sec. A.1 for details). The choice of the von Mises distribution is to satisfy the periodic boundary condition imposed by the wrapped event phase. Modeling the neutrino emission with an aperiodic Gaussian can cause a loss of statistical power, especially if the emission phase profile is wide and/or if the peak is close to the boundary. The parameters to be fitted in the temporal term are the phase peak, Φ 0 , the concentration, κ which represents the spread of the events in phase, and the orbital period, P . A Gaussian prior centered at the measured period, P exp is added to the likelihood ratio to facilitate the optimization of P (see Sec. A.1 for details).
In addition to the initial source list discussed at the beginning of this section, binary sources from the 8 yr Fermi-LAT 4FGL gamma-ray catalog (Abdollahi et al. 2020) are also included in this analysis. Due to the limited sensitivity of IceCube to neutrino-induced track events in the Southern Sky, only sources with a declination greater than −5 • are selected. Since the focus is on sources with known periodicity, the next step excludes sources without a measured period, leaving 55 sources.
Results -No evidence for periodic neutrino emission is found. The most significant source is V635 Cas, which is an HMXB consisting of a Be star and a NS. The pretrial significance is 2.6 σ, which results in a post-trial p-value of 0.25. The distribution of signal-like events is shown in the top panel of Fig. 1 as a function of the phase. The full results are tabulated in Table 2.

Flare Analysis
This analysis focuses on searching for a correlation between the neutrino emission and the X-ray activity of a source. Hard X-ray light curves are used to construct the time PDF. Light curves are obtained from data reported by Swift/BAT 15-50 keV band 2 (Krimm et al. 2013) and MAXI 10-20 keV band 3 (Matsuoka et al. 2009). The Xray light curves are binned in days, and a Bayesian block algorithm is applied to find the optimal segmentation of the data and identify flares (Scargle et al. 2013). After the light curves are divided into blocks, the flux data points inside each segment can be fitted as a constant value taking into account the uncertainty of the X-ray data. The normalized blocked light curves over the total livetime then act as the time PDF, mimicking the assumed neutrino flaring. Time-related parameters introduced here are the threshold of the X-ray flux f th which removes irrelevant variation thus picks flares, and the potential time lag L t between the X-ray and the neutrino emission. A more detailed description can be found in Sec. A.2. Starting from our initial source list, sources without available Swift/BAT or MAXI hard X-ray light curves are removed. Furthermore, the variability and excess variance (see Sec. A.2) of the light curves are evaluated such that sources with weak or steady emission are also removed. This step is applied only to the X-ray data in the time frame overlapped by the neutrino data sample. As neutrino emission is more likely to be correlated with harder X-ray emission, which is more probable to be nonthermal, the Swift/BAT light curves are selected over MAXI light curves cases where both pass the selection criteria. After this selection, there are 102 sources from the initial source list left to be analyzed.
Results -No significant signal events are found in the flare analysis, and results of all sources in this analysis can be found in Table 3. The most significant source in the flare analysis is the microquasar V404 Cyg, a low-mass BH XRB, with a post-trial p-value of 0.75. There are five sub-TeV neutrino events within 1.5 • of the source at the time of the major X-ray flare in 2015, and the best-fit threshold indicates a time duration of 11 days, as seen in the bottom panel of Fig. 1. This giant X-ray flare was observed with a duration of approximately 13 days by Swift/BAT 4 .

Time-integrated Analysis
To complement the time-dependent studies, we conduct a time-integrated search for neutrinos on four notable sources: Cyg X-3, LS 5039, LS I +61 303, and SS 433. These are sources widely discussed as potential CR accelerators (see e.g., Gaisser & Stanev 1985;Aharonian et al. 2006;Christiansen et al. 2006;Torres & Halzen 2007;Reynoso & Carulli 2019;Kimura et al. 2020). Time-integrated tests on these four sources use the method described in Braun et al. (2008).
Results -In the time-integrated analysis, we did not find any significant signal. The most significant excess is found for Cyg X-3, which exhibits a post-trial p-value 0.036 after considering the 4 trials in the time-integrated analysis. Within 1 • around the source location, there are 44 events, and the most energetic among them has a deposited energy about 5 TeV, leading to a soft best-fit power-law spectrum (γ = 3.25). In the absence of any significant signal, 90% confidence level (C.L.) upper limits (ULs) are computed for the sources with assumed spectra, which are shown in Table 2, Table 3, and Table 4. There is an overlap of 30 sources between the periodic and flare analyses. V635 Cas, the source with the smallest p-value in the periodic analysis, shows clear flaring episodes in the Swift/BAT light curves, which yields a pre-trial p-value of 0.16 in the flare analysis. The most significant source in the flare analysis V404 Cyg has a pre-trial p-value of 0.92 in the periodic analysis. In both periodic and flaring analyses, Cyg X-3 is one of the top 10 sources. On the whole, all results are consistent with the null hypothesis.
For non-microquasar XRBs with an NS as the compact object such as V635 Cas, the most significant source in the periodic study, it is possible that the CRs are accelerated in the magnetosphere, and neutrinos can be produced via CR interactions with the accretion disk through pp interactions. A similar scenario has been studied in detail for A 0535+26, an HMXB with an accreting neutron star, suggesting possible periodic neutrino emission by Anchordoqui et al. (2003). The flux prediction and the UL derived from the periodic analysis are shown in Fig. 2. The UL is calculated for a powerlaw spectrum with a cutoff at 1 TeV. The current UL is comparable to the predicted flux, yet not enough to impose a constraint on the predicted neutrino flux from the model. The range of the estimated sensitivity with 15 years of IceCube data due to varying the duty cycle is also shown in Fig. 2  could likely be constrained by a future IceCube analysis with an extended livetime. The most significant source in the flare search, V404 Cyg, is a microquasar. For microquasars, relativistic jets are expected to be CR acceleration sites. Possible neutrino emission is expected from the beam dump on either radiation from the compact object itself or gas from the companion star. Here, we employ the parameters for neutrino flux prediction in Distefano et al. (2002), based on the photohadronic model of Levinson & Waxman (2001). The flare of V404 Cyg in June 2015 was observed in a multi-wavelength campaign, and the jet activity during that outburst was studied, e.g., in Miller-Jones et al. (2019); Tetarenko et al. (2019). A simple estimation of the neutrino flux using the jet model can be performed with the radio jet information when the source is in an outburst state. The comparison between the ULs and the predicted fluence is shown in Fig. 3. The CR interaction region in the source is estimated from the flaring duration, and the spectrum is assumed to follow a power-law distribution (γ = 2) with an exponential cutoff at 100 TeV, motivated by spectrum of Fermi acceleration and maximal expected   Figure 4. Red and purple lines indicate a comparison between current UL and estimated 10 yr sensitivity (dotted) & discovery potential (solid) in IceCube-Gen2. As the nearby events cut at several TeV, an exponential cutoff at 5 TeV is also applied for computing ULs. The shaded regions show predictions of pp and pγ scenarios. The inclusion of a cutoff is also to be compared to the shaded pink region, which includes a cutoff of CR energy at 100 TeV with the spectral index ranging from 2.4-2.7. The gray shaded region shows the uncertainty from the collision radius.
Galactic CR energies. Regarding jet models, the magnetic energy and the kinetic energy in the jets are usually assumed to be in equipartition, leading to large magnetic fields, which may imply a strong attenuation of pions and muons at high energies, decreasing the neutrino flux compared to the simple calculations (Reynoso & Romero 2009). An alternative scenario considers neutrino emission through jet-wind interactions further away from the jet base (Reynoso & Romero 2009). Cyg X-3, the most significant source in the timeintegrated search, is one of the microquasars identified as a gamma-ray source early in observations. Many predictions have been calculated in the past decades assuming different models for high-energy emission from microquasars (see e.g., Gaisser & Stanev 1985;Kolb et al. 1985;Berezinskii et al. 1986;Bednarek 2005). For comparison, we take Baerwald & Guetta (2013) and Sahakyan et al. (2014), which discussed the general pγ and pp scenarios. It is important to mention that Cyg X-3 lies in the direction of the Cygnus X region and is close to the Cygnus OB2 association, which is potentially one of the PeV point sources detected by LHAASO (Cao et al. 2021). Source contamination from the Cygnus X complex cannot be excluded.
For TeV XRBs, in a recent ANTARES time-integrated point-source search (Illuminati 2021), HESS J0632+057 has a pre-trial p-value 0.02 while this periodic analysis finds a pre-trial p-value 9.0 × 10 −3 , the second most significant in the analysis. However, the flare search gives a pre-trial p-value of 1 while it is not included in the time-integrated search. Nevertheless, the hadronic component of the TeV gamma-ray observation cannot be constrained, and the significance is not large enough for any conclusion. For SS 433, more years of observation are needed to constrain the hadronic fraction of the observed TeV emission by HAWC Abeysekara et al. (2018) and flux prediction in e.g., Kimura et al. (2020) while prediction in Distefano et al. (2002) is strongly constrained. More data are also needed to constrain LS 5039, which lies in the Southern Sky, and LS I +61 303, where the neutrino flux calculation from e.g., Torres & Halzen (2007) is not constrained by the upper limit.
The next generation of the IceCube Observatory, IceCube-Gen2, will provide a factor of eight extension in volume (Aartsen et al. 2020b), leading to an expected five-fold increase in the effective area compared to Ice-Cube, corresponding to an improvement in sensitivity by the same order. Here, we extend the study to IceCube-Gen2 and estimate the sensitivity and discovery potential for V404 Cyg, as an example of a flaring source and Cyg X-3, for persistently emitting sources. The estimated improvement can be seen in Fig. 3 and Fig. 4. Here, the effective areas of muon tracks are computed from the proposed IceCube-Gen2 configuration, and the projection is evaluated with the same method as in Aartsen et al. (2020b) without considering contribution from the existing IceCube detector. In comparison with theoretical calculations, it demonstrates the potential to either identify those sources or better constrain models in the future. Therefore, providing a better estimation of the hadronic component of TeV gamma-ray sources.

CONCLUSION & OUTLOOK
In this paper, we presented a comprehensive study on neutrino emission from XRBs, a long-standing candidate for the Galactic sources of CRs and neutrinos. With no significant signal events found, we set ULs on the neutrino emission in the scenarios presented with a discussion. Our estimation with the improved sensitivity of IceCube-Gen2 demonstrates the potential of future detection, presenting a promising outlook of identifying XRBs as Galactic cosmic-ray accelerators in the upcoming years. The general form of the likelihood in this work can be written as where i indicates the i-th event, n s is the number of signal events, N is the total number of events, and S i , B i are the signal and background PDFs evaluated for event i, respectively. The signal PDF is the product of three terms, where P, E and T are the PDFs of space, energy and time, respectively. x i is the reconstructed incoming direction, σ i is the estimated angular uncertainty, E i is the reconstructed muon energy, and t i is the arrival time of the event. The spatial clustering of signal events, represented by P, can be modeled as a 2D Gaussian distribution; the energy term E is modeled by a simple power-law spectrum with spectral index γ, and the detector response is determined by the Monte Carlo simulation. In this work, the fitting range of γ is set to 1-4.
constructed by binning the experimental data in reconstructed declination and energy, based on the assumption that background events dominate the data. It is independent of the right ascension due to the Earth's rotation. The time distribution of the background is assumed to be uniform, considering that the effect of seasonal variation is negligible.

A.1. Periodic Search
The signal time PDF is expressed as a phase PDF modeled by the von Mises distribution where φ i is the phase of the event which is a function of the arrival time, t i and the orbital period, P of the source using the equation φ i = (t i −T 0 )/P , with T 0 set to MJD 55690. The value T 0 is chosen to reduce the error accumulated by the possible uncertainty in P . κ measures the concentration of events within a period, and Φ 0 is the mean phase of the distribution. κ, Φ 0 , P , along with n s and γ, are free parameters. The TS in this analysis is Similar to the previous time-dependent analysis by Aartsen et al. (2015), the likelihood ratio is modified by a marginalization term 2π/σ Ψ (κ) to prevent the minimizer from showing biases towards short flares, where σ Ψ is the standard deviation of the von Mises distribution given by The second term is the prior probability of the orbital period of the source. It is approximated by a normal distribution centered at the measured orbital period P exp , and the variance σ 2 P is chosen as the reported uncertainty of P exp . The role of the prior is to facilitate the minimization of the orbital period with the information provided by the electromagnetic observations of the sources.

A.2. Flare Search
The method is similar to the IceCube triggered flare search in Aartsen et al. (2015). The Bayesian block algorithm is applied to the daily binned X-ray data to characterize statistically significant variations, thus remove noise and identify flares. The signal time PDF is constructed from the Bayesian blocked X-ray light curve and can be written as where f is the count rate of the light curve, which is a function of time and f th is the count rate threshold. Only the light curve above the threshold is kept, which picks out the flares we care about while the irrelevant variation below the threshold is removed. t lag is the time lag, which accounts for the time difference between X-ray flares and neutrino flares. This time lag comes from two factors: the first is the time difference between X-ray and neutrino emissions, which is not expected to be a large value; the second is the uncertainty of binning, as X-ray light curves used in this analysis are daily binned and the real emission time is unknown. We allow a uniform time lag prior of ±0.5 days for fitting. The TS is defined as the maximal log-likelihood ratio TS = −2 log[(L(n s = 0)/L(n s ,γ,f th ,t lag )].
The source selection is based on the variability which describes how variable the data is, and the excess variance which identifies flares of X-ray light curves (Krimm et al. 2013) where F X, i is the individual measurement of a source flux, F X, avg is the weighted average and N X is the total number of measurements. σ X, i includes both statistical and systematic errors.

B. SOURCE TABLES
The initial source list includes 301 X-ray binary sources from the HMXB and LMXB catalogs. On top of those, 7 binary sources from the TeVCat, which are not listed in the previous two, are added to the initial list. In the periodic search, 31 sources from the HMXB, 22 from the LMXB, and 1 from TeVCat HESS J0632+057, which have resolved orbital periods and declinations above −5 • are selected. An extra source, 3FGL J0212.1+5320 from the 4FGL catalog, is added to form the analysis source list of 55 sources. In the flare search, there are 43 HMXB and 58 LMXB from the catalogs covering the whole sky with one extra HMXB HESS J0632+057 from TeVCat, which compose the total 102 sources after the selection. There are 30 sources overlapping in the two analyses.   Note-Sources for the periodic point-source analysis with their best-fitted TS,ns,γ,κ,Φ0, P (in days), and pre-trial p-values. The 90% C.L. flux ULs are parameterized as dNν µ +νµ /dEν = Φ 90%, E −γ νµ+ν (Eν /TeV) −γ where γ indicates the assumed spectral index of a power-law spectrum, and their units are 10 −12 TeV −1 cm −1 s −1 . The equatorial coordinates are provided in the J2000 epoch, which also applies to other tables in this paper.   Note-Sources for the flare analysis with their best-fitted TS,ns,γ, thresholdf th , time lagT lag and pre-trial p-values. The threshold has units 10 −2 count cm −2 s −1 and the superscript m indicates the use of MAXI light curves while others use Swift/BAT light curves. The time lag is represented in days. The 90% C.L. fluence ULs are parameterized as E 2 dNν µ+νµ /dEν = F 90%, E −γ νµ+νµ (Eν /TeV) −γ · 10 −3 TeV cm −2 , where γ indicates the assumed spectral index of a power-law spectrum. Note-Sources for the time-integrated analysis with their best-fitted TS,ns,γ and pre-trial p-values.The 90% C.L. flux UL are parameterized as dNν µ+νµ /dEν = φ 90%, E −γ νµ+νµ (Eν /TeV) −γ · 10 −12 TeV −1 cm −1 s −1 , where γ indicates the assumed spectral index of a power-law spectrum.