HAWC Search for High-mass Microquasars

We select target sources based on two criteria: (i) it is a confirmed XRB system with steady radio emission, i.e., a microquasar, within the sky coverage of HAWC, and (ii) it does not present transient X-ray outbursts like XTE J0421+560 (Frontera et al. 1998). Applying these conditions to the high-mass XRB catalog (Liu et al. 2006), we are left with four HMMQs as target sources: LS 5039, Cyg X-1, Cyg X-3, and SS 433. Although LS I +61° 303 may also seem to satisfy our criteria, it is at the edge of the HAWC field of view (FOV). Due to the poor detector sensitivity in that region, we do not include LS I +61° 303 in the target list.

Table 1 lists the relevant properties of the four HMMQs studied in this paper.

HAWC is a high duty cycle, wide-FOV particle sampling array consisting of 300 water Cerenkov detectors (WCDs) covering a combined geometrical area of ∼22,000 m² (Abeysekara et al. 2017a). It is located at a latitude of ∼19° N and at an altitude of ∼4100 m in Mexico. Each WCD contains 200,000 liters of purified water, and four upward-facing photomultiplier tubes are anchored to the bottom (Abeysekara et al. 2017a). The data set used in this analysis consists of cumulative observational data averaged over the time period of 1523 days, and the energy of the γ-ray events is estimated from the number of hit photomultiplier tubes per gamma-ray event. The expected energy and angular resolutions are ≥20% and ≥01, respectively, based on the Crab Nebula analysis (Abeysekara et al. 2017a), where more details about the HAWC setup, data, and general source analysis procedures can also be found.

Likelihood fitting with given spatial and spectral models is used to compute the γ-ray energy spectrum. In each energy bin i, a simple power-law spectral model is used to describe the γ-ray spectrum,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{i}={A}_{i}{\left(\displaystyle \frac{E}{{E}_{i,\mathrm{piv}}}\right)}^{-{\alpha }_{i}},\end{eqnarray} \tag{ 1 }$$

where Φ_i is the differential flux at the pivot energy E_i,piv, A_i is the flux normalization, E is the photon energy, and α_i is the spectral index. For this analysis, we use four quasi-differential energy bins as listed in Table 2. Within each bin, we adopt a spectral index α_i = 2.7, which is a good approximation for point-like HAWC sources (Abeysekara et al. 2017b). The systematic uncertainties due to the unknown spectral index and detector response functions will be discussed below. Since the binary systems have a typical size of 0.1 au and are located at a distance of several kiloparsecs, a point-source morphology is adopted for all target sources.

All four target sources are located in source-confused regions close to the Galactic plane with several nearby TeV sources. The HAWC significance maps of each region, with nearby sources labeled, are shown in Appendix C, Figure 4. The residual maps, as shown in Figure 1, are obtained with the following steps. We first fit background sources using their known locations and spectral indices from the 3HWC Catalog (Albert et al. 2020). In particular, we fit point-like background sources such as 3HWC J1819–150 and 3HWC J1913+048 with a point-source model. The regions of interest also contain four extended sources. We use a simple Gaussian morphology for 3HWC J2006+340 and 3HWC J1908+063. The 3HWC J1825–134 area consists of two pulsar wind nebulae, HESS J1825–137 and HESS J1826–130 (Abdalla et al. 2019), positioned above and below the location of 3HWC J1825–134. Hence, we apply an asymmetrical Gaussian morphological model to 3HWC J1825–134 with its semimajor axis positioned along the line joining the three VHE source locations. Finally, the Cygnus cocoon's gamma-ray profile is "flat" (Hona et al. 2020). Hence, we adopt a disk-like morphological model for 3HWC J2031+415.

Zoom In Zoom Out Reset image size

The obtained best-fit models for the nearby sources are then subtracted from the original HAWC 1523 transit maps to produce the residual maps as shown in Figure 1. Then, we fit for the flux normalization of each HMMQ to find their flux upper limits.

We calculate a test statistic (TS) for γ-ray detection based on the logarithm of the likelihood ratio when fitting with the residual maps with and without the target source in all energy bins,

$$\begin{eqnarray}&&\mathrm{TS}\equiv 2\left[\mathrm{ln}{ \mathcal L }(\hat{A})-\mathrm{ln}{ \mathcal L }(A=0)\right],\end{eqnarray} \tag{ 2 }$$

where ${ \mathcal L }$ is the Poisson likelihood function and $\hat{A}$ is the best-fit normalization found from the maximum-likelihood estimators. We obtain a priori statistical significance for a given location in the sky via

$$\begin{eqnarray}&&\sigma \approx \pm \sqrt{\mathrm{TS}}.\end{eqnarray} \tag{ 3 }$$

The best-fit normalization, $\hat{A}$, is used as an input to a Markov Chain Monte Carlo (MCMC). MCMC then estimates the distribution of the posterior likelihood around the maximum value of the likelihood with a positive uniform prior assumed. From the obtained MCMC distribution, we can finally compute the 95% credible upper limit on the flux normalization.

The HAWC data analysis involves forward-folding of the assumed morphological and spectral models through the detector response to obtain the expected gamma-ray counts. Following Abeysekara et al. (2019), we evaluate the detector systematic uncertainties by applying various versions of the detector response. Also, different spectral indices between 2.0 and 3.0 with an interval of 0.1 are applied to study the source spectrum. The systematic errors on the flux normalizations due to different detector responses and astrophysical spectral indices are computed for each source at each quasi-differential energy bin and for one full-energy bin containing data from all four bins. The errors are shown in Table 4 in Appendix E.

Due to their similarity in the source structure, such as the accretion disk–jet configuration, and in the radiation background, such as thermal photons from donor stars of similar star type, temperature, and size, the HMMQ population could, in principle, produce gamma rays with one same mechanism (Dubus 2013). By combining the observations of all HMMQs in the HAWC FOV, we can constrain the common factors that impact the γ-ray production in these microquasars.

Below we consider two generic models, referred to as scenarios I and II, for VHE γ-ray emission in microquasar jets. In the first scenario, we assume that γ-ray luminosity is proportional to the kinetic power of the jets,

$$\begin{eqnarray}&&{L}_{\gamma }={\epsilon }_{\gamma }\,{L}_{\mathrm{jet}}.\end{eqnarray} \tag{ 4 }$$

This is a general assumption that may be satisfied by different γ-ray production models such as neutral pion decay from hadronic interactions. The γ-ray flux in scenario I can be written as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\gamma }=\displaystyle \frac{{\epsilon }_{\gamma }\,{L}_{\mathrm{jet}}}{4\pi \,{D}^{2}}\,{K}_{p}{\left(\displaystyle \frac{E}{{E}_{\mathrm{piv}}}\right)}^{-p},\end{eqnarray} \tag{ 5 }$$

where D is the distance to the source, ${K}_{p}=(2-p){E}_{\mathrm{piv}}^{-p}/({E}_{\max }^{2-p}-{E}_{\min }^{2-p})$ is a normalization factor for spectral index p, and ${K}_{p}={E}_{\mathrm{piv}}^{-2}/\mathrm{log}({E}_{\max }/{E}_{\min })$ for p = 2. Also, ${E}_{\min }=1$ TeV and ${E}_{\max }=100$ TeV are the boundaries of the energy bin used for the stacking analysis (instead of the quasi-differential bins used in Section 2.2), with E_piv = 7 TeV as the pivot energy. L_jet and D of the target sources are listed in Table 1.

In the second scenario, we consider the model summarized in Dubus (2013), where gamma rays are produced when relativistic electrons accelerated by the jets upscatter optical photons from the donor star. See Appendix A for more details regarding the modeling of γ-ray production. In this model, the inverse Compton emission of an HMMQ is expected to peak at TeV energies, and the corresponding synchrotron emission is typically at 10 keV–10 MeV. The energy fluxes of the two components, F_syn and F_IC, are connected by

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{\mathrm{syn}}}{{F}_{\mathrm{IC}}}\approx \displaystyle \frac{{u}_{B}}{{u}_{0}\,{f}_{\mathrm{KN}}},\end{eqnarray} \tag{ 6 }$$

where u₀ is the energy density of the radiation field of the star (Equation (A1)) and u_B = B²/8π is the magnetic energy density. Since the stellar radiation field is in the optical band, the inverse Compton emission of VHE electrons is in the Klein–Nishina regime. The unitless f_KN factor, evaluated at the inverse Compton break energy, E_IC,bk (Equation (A4)), accounts for the suppression of the inverse Compton cross section.

The γ-ray flux in scenario II can be expressed as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\gamma }=\displaystyle \frac{{F}_{\mathrm{syn}}\,{u}_{0}\,{f}_{\mathrm{KN}}}{{u}_{B}}\,{K}_{p}{\left(\displaystyle \frac{E}{{E}_{\mathrm{piv}}}\right)}^{-p}.\end{eqnarray} \tag{ 7 }$$

We estimate the synchrotron flux F_syn using the measured X-ray (or MeV γ-ray) energy flux, F_obs,bk, between 0.1 E_syn,bk and 10 E_syn,bk, where E_syn,bk (Equation (A8)) is the peak energy of the synchrotron emission suggested by models fitted to the multiwavelength data. The energy density of the radiation field u₀ is derived from observed properties, including stellar temperature, radius, and separation from the compact object as listed in Table 1. The u₀ and f_KN used in the analysis are listed in Table 3 in Appendix A.

In both scenarios, the γ-ray flux of the ith source can be written as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{i}=K\,{C}_{i}{\left(\displaystyle \frac{E}{{E}_{\mathrm{piv}}}\right)}^{-p},\end{eqnarray} \tag{ 8 }$$

where C_i is the source-dependent contribution factor, ³⁴ ${C}_{i}={K}_{p}\,{L}_{\mathrm{jet},{\rm{i}}}/4\pi \,{D}_{i}^{2}$ in scenario I and C_i = K_p u_0,i f_KN,i F_syn,i in scenario II. K is the "weighting factor" shared between the sources, specifically K = _γ in scenario I and K = 1/u_B in scenario II.

We perform a likelihood fit for each target HMMQ to obtain their best-fit flux normalization. The sources are then stacked in the likelihood space weighted by their relevant contribution factors C_i ,

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }(p,K)=\displaystyle \sum _{i}\mathrm{ln}{{ \mathcal L }}_{i}(p,K,{C}_{i}).\end{eqnarray} \tag{ 9 }$$

The credible interval of the "linked" flux normalization for a given p is obtained using MCMC following the steps described in Section 2.2. By scanning the index p between 2.0 and 3.0 with 0.1 intervals, we can find the upper limit of $\hat{K}$ and the best-fit $\hat{p}$ that corresponds to the peak of the maximum log-likelihood $\mathrm{ln}{ \mathcal L }(p,\hat{K}(p))$. The largest difference in K obtained when varying p is used as an estimation of the statistical error due to the scanning of the index. The detector systematic uncertainties are evaluated by applying various versions of the detector response for the cases with the best-fit spectral index.

Finally, the stacked flux is used to derive the limits on the weighting factor K. Note that the stacked flux depends on the definition of the contribution factor in a physical model. The fluxes in different scenarios are not directly comparable.

Figure 2 shows the spectral energy distribution of our target sources ranging from X-rays to multi-TeV gamma rays. LS 5039 is currently the only source in our list that has been detected at TeV energies (Mariaud et al. 2016). Our limits below 10 TeV are consistent with the observation of this source by the IACTs. For Cyg X-1 and Cyg X-3, the upper limits from MAGIC (Zdziarski et al. 2017) and VERITAS (Archambault et al. 2013) are more constraining at 1 TeV but approach the HAWC upper limits as the energy goes up. Finally, for SS 433, our limits are slightly less constraining in the first quasi-differential bin but become comparable to the combined MAGIC-H.E.S.S. data (Ahnen et al. 2018) at higher energies. This could be due to a potential contribution from the SS 433 west lobe (Abeysekara et al. 2018), which is not included in the 3HWC Catalog (Albert et al. 2020).

Zoom In Zoom Out Reset image size

In Figure 2, containment bands are displayed to indicate the HAWC sensitivity at each location. A point-source model is fitted in the empty regions of the sky along the same decl. band as the target HMMQ to calculate the expected upper limits containing 68% and 95% in yellow and green, respectively. For the calculation of the sensitivities, regions with VHE gamma-ray sources such as the Galactic plane have been excluded. Indeed, our upper credible intervals, in red, are at most about 2σ above the expected HAWC limit if there was no emission (dashed black line). Hence, we do not have a clear detection of the HMMQs.

For the four sources discussed in this work, HAWC provides the most stringent upper limits above 10 TeV.

In neither stacking analysis does the combination of the four sources result in a significant detection. However, the stacked flux limits allow us to set limits on parameters of the scenarios.

In scenario I, we find the best-fit flux norm K ∑_i C_i = (2.4 ± 1.1_stat ± 0.4_sys) × 10⁻¹⁵ TeV⁻¹ cm⁻² s⁻¹ with TS = 4.4 and the best-fit spectral index $\hat{p}=2.2$. It corresponds to a 95% confidence interval (C.I.) limit on the stacked flux Φ_γ (E_piv) = 4.8 × 10⁻¹⁵ TeV⁻¹ cm⁻² s⁻¹. Through Equation (5), we obtain the limit on the jet emission efficiency above 1 TeV,

$$\begin{eqnarray}&&{\epsilon }_{\gamma }^{\mathrm{UL}}=5.4\times {10}^{-6}.\end{eqnarray} \tag{ 10 }$$

This TeV emission efficiency is 3–5 orders of magnitude lower than the emission efficiency of HMMQs in 0.5–10 keV X-rays, which typically reaches 10⁻³ to 10⁻¹ (Marti et al. 1998; Cadolle Bel et al. 2006). We note that ${\epsilon }_{\gamma }^{\mathrm{UL}}$ is derived using the jet power in Table 1. The L_jet values of the microquasars obtained by different works may differ by a factor of ∼2–4 (e.g., Paredes et al. 2006 and Casares et al. 2005). Note that the uncertainty in L_jet is not accounted for in the calculation of ${\epsilon }_{\gamma }^{\mathrm{UL}}$.

Our TeV γ-ray emission efficiency constrains the HE neutrino emission efficiency _ν of HMMQs. If VHE gamma rays are produced by the decay of neutral pions, the same proton–proton interaction should produce charged pions that decay into HE neutrinos with an emission efficiency _ν ≈ 3_γ /2. ³⁵ The _γ derived in Equation (10) suggests that a mean-orbital _ν ∼ 0.2 assumed by Christiansen et al. (2006) is overly optimistic. The emission efficiency also implies that neutrino detection of HMMQs is difficult with the current neutrino detectors (IceCube Collaboration et al. 2018).

In scenario II, using the model described in Section 2.3 and the INTEGRAL and COMPTEL observations of the sources, F_obs,bk, the stacking analysis yields the best-fit flux norm, K ∑_i C_i = (6.0 ± 8.8_stat ± 0.7_sys) × 10⁻¹⁶ TeV⁻¹ cm⁻² s⁻¹ with TS < 4 and the best-fit spectral index $\hat{p}=2.1$. The 95% C.I. upper limit on the stacked γ-ray flux is Φ_γ (E_piv) = 2.4 × 10⁻¹⁵ TeV⁻¹ cm⁻² s⁻¹, which corresponds to a lower limit on the magnetic field strength,

$$\begin{eqnarray}&&{B}^{\mathrm{LL}}=22{\left(\displaystyle \frac{{\epsilon }_{\mathrm{syn}}}{10 \% }\right)}^{1/2}\,{\rm{G}},\end{eqnarray} \tag{ 11 }$$

where _syn is an unknown factor denoting the ratio of the actual synchrotron emission by the electron population that emits VHE gamma rays to the total observed 10 keV–10 MeV flux, _syn ≡ F_syn/F_obs,bk.

To evaluate the dependence of our result on the γ-ray spectrum and consider that the γ-ray spectrum may not strictly follow a power-law spectrum, we also perform the analysis with a log-parabolic spectral model,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\gamma }=\displaystyle \frac{{F}_{\mathrm{syn}}\,{u}_{0}\,{f}_{\mathrm{KN}}}{{u}_{B}}\,{K}_{l}{\left(\displaystyle \frac{E}{{E}_{\mathrm{piv}}}\right)}^{-{\alpha }_{l}-{\beta }_{l}\mathrm{log}(E/{E}_{\mathrm{piv}})}.\end{eqnarray} \tag{ 12 }$$

We fix E_piv at 7 TeV and scan the indices α_l between 2.0 and 5.0, with 0.5 intervals, and β_l between 0.1 and 2.1, with 0.4 intervals, to find the best-fit flux norm K_l . The fit to the data, however, does not significantly improve with an extra parameter. The lower limit on the magnetic field strength is obtained to be ${B}_{l}^{\mathrm{LL}}=15{({\epsilon }_{\mathrm{syn}}/10 \% )}^{1/2}$ G, which is comparable to the limit from the power-law assumption.

The derived magnetic field strength agrees with the finding of Dubus et al. (2015), where B ≈ 20 G was obtained by fitting a relativistic hydrodynamics model to the multiwavelength observation of LS 5039. Dubus et al. (2015) conclude that a high B is unavoidable to explain the COMPTEL flux level of LS 5039. Our result extends the conclusion to all HMMQs accessible to HAWC and suggests that the large gap between the energy flux in 10 keV–10 MeV and that in VHE gamma rays could be a universal feature of HMMQs. Such a high magnetic field challenges the existing models of γ-ray binaries (Bosch-Ramon et al. 2008) and suggests that the synchrotron component is a small fraction, _syn ≲ 10%, of the observed flux between 10 keV and 10 MeV. A few caveats should, however, be noted when interpreting this result, as discussed in Section 4.