SEARCHES FOR PERIODIC NEUTRINO EMISSION FROM BINARY SYSTEMS WITH 22 AND 40 STRINGS OF ICECUBE

Recently, the observation of Very High Energy (VHE) γ-ray emission from several high-mass X-ray binaries has established a new subclass of VHE or "γ-ray-loud" binaries. While much of the evidence from multi-wavelength observations favors leptonic emission, it is likely that a hadronic component is also accelerated in the jets of these binary systems. These binary systems are powerful accelerators with jets capable of accelerating cosmic rays and can have a role as PeVatron accelerators of galactic cosmic rays. The observation of neutrino emission would be clear evidence for the presence of a hadronic component in the outflow of these sources. Four such binary systems, PSR B1259-63, LS 5039, HESS J0632+057, and LS I +61° 303 have been identified as persistent TeV γ-ray emitters (Aharonian et al. 2005a, 2005b, 2007; Albert et al. 2006), while Cygnus X-1 is a possible candidate. The binary pulsar system PSR B1259-63 was also recently discovered to have periodic emission in GeV photons (Abdo et al. 2011). PSR B1259-63 is formed by a B2Ve star orbited by a young 48 ms pulsar (Tavani & Arons 1997) both exhibiting a strong wind. As observed in Aharonian et al. (2005a), its VHE emission could come from inverse Compton scattering on shock-accelerated leptons from the interaction zone between the pulsar and wind from the star though a hadronic interpretation cannot be excluded. On the other hand, the driving factor of the VHE emission in Cygnus X-1, most probably a black hole orbiting a supergiant O9.7 star (Ziółkowski 2005), could be the interaction of the black hole with the strong stellar wind of the star. However, there has been no other evidence for steady VHE emission of Cygnus X-1, though a VHE flare of about 1 hr was observed (Albert et al. 2007). HESS J0632+057  was recently seen to have a periodic modulation in X-rays by Swift, with heightened TeV emission coincident with the X-ray maximum (Bongiorno et al. 2011).

LS I +61° 303 remains a mystery even after four decades of observations over a wide range of wavelengths, from radio (Gregory 2002), soft and hard X-ray (Torres et al. 2010; Li et al. 2011; Zhang et al. 2010; Sidoli et al. 2006; Harrison et al. 2000), GeV (Abdo et al. 2009a), and TeV photons (Albert et al. 2009; Acciari et al. 2009). The best measurement of its period, P₁ = 26.4960 ± 0.0028 days, comes from radio data (Gregory 2002 and references therein) with the orbital zero phase taken by convention at JD 2,443,366.775 (Gregory & Taylor 1978), but the same modulation has also been detected in other wavelengths, notably in the GeV/TeV band emission (Abdo et al. 2009a; Albert et al. 2009). Together with LS 5039, discovered in the TeV γ band (Aharonian et al. 2005a), LS I +61° 303 lacks strong evidence supporting the black hole or neutron star nature of the compact object. This prevents their clear classification as microquasars or pulsar systems (Paredes 2011; Paredes et al. 2000, and references therein). A discussion of the different theoretical models for these systems is presented in Bosch-Ramon & Khangulyan (2009). As with all of the binary systems above, the detection of multi-TeV neutrinos would complement the VHE photon observations and unequivocally prove the existence of hadronic acceleration.

The analysis in this paper has been performed using a likelihood method in which the underlying hypothesis is that the neutrino emission is periodically modulated due to the geometry of the X-ray binary systems. Neutrinos would be produced by a beam of hadrons accelerated by the compact object and interacting with the matter of the massive star and its atmosphere. The periodic modulation would be connected to the orbital motion of the system. This modulation is observed in photons from radio to X-ray and in the VHE band. The analysis is designed to incorporate only minimal assumptions regarding the neutrino emission. The period is fixed to that observed in an electromagnetic band, while the phase and duration of neutrino emission are free parameters and not constrained to match the photon emission. This is to account for the fact that photons can be absorbed when the accelerator is behind the large star of the binary system, while neutrino production can be enhanced if enough matter is crossed. The neutrino energy spectrum is fit with a simple power law with the index also a free parameter.

This paper is organized as follows. In Section 2, we describe the IceCube observatory and the data taken with two detector configurations (Achterberg et al. 2006). The analysis method is described in Section 3, and the expected sensitivity and discovery potential are shown. To avoid bias, the search has been performed in a blind fashion by defining cuts before looking at the true times (equivalently, the right ascension values) of the final event sample. In Section 4, we present the results of the search performed on a catalog of seven galactic binary stars in the northern sky. The selected objects are considered as microquasars in Distefano et al. (2002), where their expected emission of neutrinos is calculated. While that paper is not specifically about the periodic emission from these sources, the objects considered there are promising neutrino emitters for which radio observations allow identification of jet parameters such as the Lorentz factor and the luminosity of the jet. All of the sources considered are located in the Northern Hemisphere where IceCube is most sensitive (Abbasi et al. 2011).

Construction of the IceCube Neutrino Observatory started with a first string installed in the austral summer of 2005/2006 (Achterberg et al. 2006) and was completed in 2010 December. It is composed of an array of 86 strings with a total of 5160 Digital Optical Modules (DOMs) instrumented between a depth of 1450 and 2450 m. The deep ice detector is complemented at the surface by the extensive air shower array IceTop. Each DOM consists of a 25 cm diameter Hamamatsu photomultiplier (Abbasi et al. 2010) and electronics to digitize the photomultiplier tube (PMT) output voltage (Abbasi et al. 2009b), all in a spherical, pressure-resistant glass housing. IceCube observes the Cherenkov photons emitted by relativistic charged particles produced in high-energy neutrino interactions. The PMT signals are digitized and the charge and arrival time of photons are measured. The data taking and performances of the detector during the two seasons analyzed here are described in Abbasi et al. (2009a, 2011). The main background comes from down-going muons due to cosmic-ray interactions in the atmosphere above the detector: in the 40 string configuration these trigger the detector at a rate of about 950 Hz and in the 22 string configuration at about 350 Hz. During the austral summer the atmosphere above the South Pole gets warmer and thinner and the probability that pions generated in cosmic-ray air showers decay rather than interact increases (Tilav et al. 2009), causing the trigger rate to vary by about ±10%. A series of event selections and higher-level event reconstructions are applied to remove these downward-going events, while retaining upward-going tracks from muons induced by neutrinos which crossed the Earth. At the final level of analysis, this remaining background of upward-going atmospheric neutrinos comes from many different directions on the other side of the Earth. Temperature effects average over a wide terrestrial region and the seasonal modulation is only a few percent. This variation has a period of one year, much longer than any period considered in this search.

The searches presented here used two data samples. The data taken with the 22 string configuration have a livetime of 275.7 days, 89% of the operation period from 2007 May 31 to 2008 April 5, or Modified Julian Date (MJD) 54251–54561. The sample is described in Abbasi et al. (2009a) and consists of 5114 events, which are mostly neutrino-induced upward-going muons with declinations from −5° to +85°. The deadtime is mainly due to test and calibration runs during and after the construction season. The livetime of the 40 string data used in the analysis is 375.7 days, which is 92% of the nominal operation period from 2008 April 5 to 2009 May 20 (MJD 54561–54971). The handling and processing of the data to obtain the final neutrino candidate event sample are fully described in Abbasi et al. (2011). The final 40 string sample contains 36,900 atmospheric neutrino and muon events distributed over the whole sky of which 14,121 events are upward-going (the rest are downward-going events from the southern sky, used for neutrino source searches strictly above PeV energies). The median angular resolution for the final sample for energies greater than 10 TeV is <1°. The energy of each event is estimated using the density of photons along the muon track due to stochastic energy losses of pair production, bremsstrahlung, and photonuclear interactions which dominate over ionization losses for muons above 1 TeV. The energy resolution is about 0.3 in log₁₀ of the muon energy in the detector between 10 TeV and 10⁵ TeV. The estimated muon energy is a lower bound on the primary neutrino energy, since for interactions that occur outside the detector the muon loses energy over an unknown distance before reaching the detector. (Energy distributions used internally within the analysis therefore refer to the observable muon energies.) The muon neutrino flux upper limits at 90% CL for time-integrated searches (depending on declination) are between E²dN/dE ∼ 3–20 × 10⁻¹² TeV cm⁻² s⁻¹ in the northern sky where the sources considered in this paper are located.

The 22 and 40 string data samples used in this paper were also used to look for bursting neutrino sources in Abbasi et al. (2012) where the stability of the data taking is discussed in detail. Azimuthal geometry effects of the 22 and 40 string IceCube detectors (due to the fact that they are more elongated in one direction than in others) and the rotation of the Earth interfere constructively for source periods that match to multiples of a half a sidereal day, which is not the case for any of the source periods tested.

The limits in this paper were produced assuming a flux of only muon neutrinos and antineutrinos at the Earth with simulated energies from 10⁸ to 10¹⁹ eV. For standard neutrino oscillations over astronomical distances (Athar et al. 2000), equal fluxes of all neutrino flavors at the Earth are expected from a source producing neutrinos via pion decay with a ratio of ν_e: ν_μ: ν_τ = 1: 2: 0. For the assumption of equal fluxes of muon and tau neutrinos at the Earth, the resulting upper limits on the sum of both fluxes are about 1.7 times higher than if only muon neutrinos are considered (Abbasi et al. 2011). This sets better limits than the expected factor of two due to oscillations if no tau neutrinos were detectable. This is due to the tau decay channel into muons with a branching ratio of 17.7%. In addition to this, tau leptons with energy greater than several PeV that may travel far enough to be reconstructed as tracks in IceCube before decaying. For an E⁻² neutrino spectrum, the contribution due to the detectable tau neutrino flux for sources at the horizon is 10% and rising to 15% for sources in the Northern Hemisphere.

The main systematic uncertainties on the flux upper limits come from photon propagation in ice, absolute DOM efficiency, and uncertainties in Earth's density profile and muon energy loss. For an E⁻² spectrum, the estimated total uncertainty is about 16% (Abbasi et al. 2011). They are included in the upper limit calculations following the method of Conrad et al. (2003) with the modification described in Hill (2003).

The likelihood method used in this analysis was described in full detail and demonstrated in Braun et al. (2008, 2010), and applied to the 40 string data in Abbasi et al. (2012). In the likelihood ratio method, the data are modeled as a combination of signal and background populations. The probability density functions (PDFs) for signal and background consist of three terms: a space term, an energy term, and a time term. The first two are implemented in the same way as in Abbasi et al. (2011). For signal, the space term characterizes the clustering of event directions around the hypothesized source location (effectively, the point-spread function for the reconstructed muon, since the interaction angle between the incoming neutrino and outgoing muon is subdominant at the energies of these data samples). For background, the space term is simply estimated by time scrambling of the real data. The energy term for background is similarly a PDF built from the energy estimates of events in the real data (selected from a declination band similar to the declination of the source being searched for). For signal, distinct energy PDFs are constructed for simulated events arising from a range of neutrino source spectra from E⁻¹ to E⁻⁴. The chief purpose of the energy term in this search is not to determine the spectrum of the source (if one were detected). Rather, it is to enhance the detectability of a source if its spectrum is relatively hard (e.g., E⁻²) by leveraging the difference in the energy distribution of the signal events compared to nearby background events, which are primarily atmospheric neutrinos with a soft (∼E^−3.7) spectrum.

The third term in the PDF incorporates timing information. For signal, a periodic emission with Gaussian time profile is assumed. The period is fixed to that determined by photon observations, while the phase and duration of the neutrino emission are left as free parameters. A Gaussian shape is used for the profile to provide a smooth function with the fewest assumptions about the exact time profile of the neutrino emission. The time PDF for the ith event can thus be expressed as

$$\begin{equation} S^{\mathrm{time}}_{i} = \frac{1}{\sqrt{2\pi }\sigma _T} e^{ -\frac{ \left| \varphi _i-\varphi _0 \right|^2}{2\sigma _T^2} }, \end{equation} \tag{ 1 }$$

where σ_T is the width of the Gaussian, φ_i is the phase of the event, and φ₀ is the phase of the peak of the neutrino emission. The fit parameters are σ_T and φ₀. For background, the time term is a flat function because in the absence of detector biases the background events are randomly distributed in time.

For each candidate source, the likelihood ratio analysis finds best-fit values for four parameters: the number of signal events, the spectral index of the signal, the peak phase of the signal, and its duration. An initial estimate of the significance is made by assuming that the likelihood ratio follows a χ² distribution and converting to a (pre-trial) p-value. To ensure a robust estimate of the final significance, however, this assumption is not used, and a correction for the number of trials is also included. For the final significance, the analysis is performed on time-scrambled data and the same catalog of sources. The final post-trial p-value is given by the fraction of analyses which yield a smaller (pre-trial) p-value for any of the sources in the catalog.

We calculate the sensitivity and median upper limit at the 90% confidence level using the method in Feldman & Cousins (1998). The discovery potential is the flux required to achieve a p-value less than 2.87 × 10⁻⁷ (5σ of the upper tail of a one-sided Gaussian) in 50% of trials. It should be noted that the threshold significance to claim a discovery in IceCube is set to 5σ. Figure 1 shows the sensitivity and the discovery potential for the analysis, together with the corresponding values from the time-integrated search (Abbasi et al. 2009a).

Zoom In Zoom Out Reset image size

Compared to the time-integrated analysis, searching for periodicity in neutrino emission results in a better discovery potential if the duration of the emission σ_T is less than about 20% of the total period (see Figure 1). As the time-dependent search adds two additional degrees of freedom to the analysis, the discovery potential is, on the other hand, roughly 10%–15% better using the time-integrated search if neutrinos are actually emitted at a steady rate or over a large fraction of the period. For both the 22 string and 40 string analyses, if the emission has a σ_T of 1/50, the method requires about half as many events for discovery as the time-integrated search.

The seven predefined sources, listed in Table 1, were used for the initial search with the 22 string data from 2007–2008. The most significant outcome in this sample was for the source SS 433, with a pre-trial estimated p-value of 6%. In identical analyses of time-scrambled data, we find at least one of the seven tested sources to be more significant in 35% of the analyses.

Table 1. Results of the Analysis with the IceCube 40 String Data Sample from 2008–2009

Source	Period	p-value	T0	$\hat{\varphi }_{o}$	$\hat{\sigma }_{T}$		Time-dependent UL	Time-integrated UL
	(days)	(pre-trial)	(MJD)	(phase)	(period−1)	Reference	(TeV−1 cm−2 s−1)	(TeV−1 cm−2 s−1)
Cygnus X-3	0.199679	0.0019	54896.693	0.819	0.02	Abdo et al. (2009b)	3.01 × 10−11	6.64 × 10−12
Cygnus X-1	5.5929	0.080	41874.707	0.031	0.02	Brocksopp et al. (1998)	4.08 × 10−12	7.41 × 10−12
LS I +61° 303	26.498	0.23	43366.775	0.916	0.02	Acciari et al. (2008)	1.82 × 10−11	9.78 × 10−12
GRS 1915+105	30.8	0.43	53945.7	0.502	0.045	Neil et al. (2007)	2.57 × 10−12	3.23 × 10−12
SS 433	13.0821	0.35	50023.62	0.779	0.02	Hillwig & Gies (2008)	3.15 × 10−12	3.03 × 10−12
XTE J1118+480	0.169934	0.28	52287.9929	0.985	0.13	McClintock et al. (2003)	7.29 × 10−12	8.18 × 10−12
GRO J0422+32	0.21214	0.037	50274.4156	0.831	0.02	Webb et al. (2000)	1.46 × 10−11	6.89 × 10−12

Notes. T₀ is the time of zero phase for the binary systems tested. $\hat{\sigma }_{T}$ is the standard deviation of the best-fit Gaussian, as a fraction of the period of the binary system. The value for σ_T is constrained to be larger than 0.02 to prevent the method from isolating a single event. We also include the reference used for the orbital information. The last columns give the 90% CL upper limits on the normalization of the flux for an E⁻² neutrino spectrum, for the time-dependent and time-integrated hypotheses. The upper limits also incorporate a 16% systematic uncertainty.

Download table as:  ASCIITypeset image

The analysis was subsequently performed on the 40 string data from 2008–2009, which provided twice the sensitivity of the previous sample (see Figure 1). The most significant outcome in this sample was for Cygnus X-3. The pre-trial estimated p-value of this source is 0.0019. To account for trials and for the fact that the likelihood ratio is not perfectly χ² distributed, the analysis is performed again on time-scrambled data. An equivalent or more significant outcome from any of the sources is found in 1.8% of scrambled samples (expressed in Gaussian standard deviations, a 2.1σ excess), so the result is compatible with random fluctuations of the background. The best-fit peak emission is found to be at phase $\hat{\varphi }_0 = 0.82$ and $\hat{\sigma }_T = 0.02$. The best-fit number of source events is $\hat{n}_s = 4.28$ and spectral index is $\hat{\gamma }_s = 3.75$. The full results of the analysis on each source (with time-dependent and time-integrated flux upper limits) are given in Table 1.

Table 2 compares the 40 string time-integrated limits to the model predictions in Distefano et al. (2002) for each source. The model predicts the neutrino flux based on the radiative luminosity associated with the jet from radio observations in quiescent states and during flares, the durations of which are specified in Table 4 in that paper. The figure shows limits for both the persistent and time-dependent cases for a time window similar to the observed flare but not coincident to it (since IceCube was not active at the time of radio observations noted in the paper). For the persistent case of SS 433 the model predicts more than 100 events during the 40 string data taking period, a flux level which is excluded by previous searches by the AMANDA detector (Achterberg et al. 2007). Distefano et al. (2002) noted that for the specific case of SS 433, the model may be biased because the source is surrounded by the diffuse nebula W50, which can affect the estimate of the radio emission used in the model for this source. For Cygnus X-3, the IceCube limits are near the prediction with the 40 string data.

The main parameters on which the neutrino flux depends in this model are the fraction of jet kinetic energy converted to internal energy of electrons and magnetic field, η_e; the fraction of the jet luminosity carried by accelerated protons, η_p; and the fraction of proton energy in pions, f_π, which strongly depends on the maximum energy to which protons can be accelerated. We show as an example for the case of a three-day burst of Cygnus X-3 how the parameters are constrained by our result. We assume equipartition between the magnetic fields and the electrons and the proton component (η_p = η_e) for setting a constraint on f_π < 0.11. If equipartition does not apply, we assume f_π = f_{π, peak} as given in Table 2 in Distefano et al. (2002; for Cygnus X-3 f_{π, peak} = 0.12) and constrain η_p to be less than 92% of η_e. In deriving these limits we have assumed that the Lorentz factor of the jet is well known from radio measurements, but in many cases there is a large uncertainty on this parameter. Hence, our limits for the parameters of this model may have different implications that we cannot disentangle: protons may not be dominant in the jet, they may lose smaller energies into pion decay than the values considered in Distefano et al. (2002), or the Lorentz factor is lower than the value indicated in Table 1 in that paper.

The exploration of the GeV and TeV photon sky with the instruments on board the Fermi spacecraft and the ground-based Cherenkov telescopes has heralded the golden age of γ-ray astronomy. The connection to neutrino astronomy is clear: high-energy processes which cause the observed VHE emission can be responsible for the observed high-energy cosmic rays. This implies hadronic acceleration mechanisms in astrophysical sources which can result in an observable neutrino flux with giant neutrino telescopes like IceCube.

The available photon observations have made it possible to enhance the sensitivity of searches for neutrino fluxes by incorporating assumptions derived from the γ-ray data. One crucial development has been the formulation of time-dependent neutrino flux searches, postulating a connection between the time modulation of the high-energy emission and the possible neutrino flux. This assumption has increased the sensitivity of these searches in comparison to their time-averaged counterparts.

This paper has presented a search for neutrinos from objects with periodic photon/broadband emission. Seven X-ray binaries in the Northern Hemisphere were selected as candidate sources in analyses of IceCube 22 string and 40 string data. The most significant source in the catalog is Cygnus X-3 with a 1.8% probability after trials (2.1σ excess). Comparing the time-integrated limits for each source to model predictions from Distefano et al. (2002), we show that our limits can constrain the fraction of jet luminosity which is converted into pions and the ratio of jet energy into relativistic leptons versus relativistic hadrons, under some assumptions. For instance, for Cygnus X-3 and equipartition between electrons and protons, the fraction of proton energy in pions is limited to about 11%. All of the results in this paper are compatible with a fluctuation of the background.

We thank D. Guetta and E. Waxman for helpful discussions on neutrino flux prediction models. We acknowledge the support from the following agencies: U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, University of Wisconsin Alumni Research Foundation, the Grid Laboratory Of Wisconsin (GLOW) grid infrastructure at the University of Wisconsin-Madison, the Open Science Grid (OSG) grid infrastructure; U.S. Department of Energy, and National Energy Research Scientific Computing Center, the Louisiana Optical Network Initiative (LONI) grid computing resources; National Science and Engineering Research Council of Canada; Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation, Sweden; German Ministry for Education and Research (BMBF), Deutsche Forschungsgemeinschaft (DFG), Research Department of Plasmas with Complex Interactions (Bochum), Germany; Fund for Scientific Research (FNRS-FWO), FWO Odysseus programme, Flanders Institute to encourage scientific and technological research in industry (IWT), Belgian Federal Science Policy Office (Belspo); University of Oxford, United Kingdom; Marsden Fund, New Zealand; Japan Society for Promotion of Science (JSPS); the Swiss National Science Foundation (SNSF), Switzerland; A. Groß acknowledges support by the EU Marie Curie OIF Program; J. P. Rodrigues acknowledges support by the Capes Foundation, Ministry of Education of Brazil.